Hello,
I am trying to fit my model (initial model – it may not actually produce a good fit) to a set of experimental data but the fminsearch seems to produce a rather odd result.
The model aims to capture cell proliferation when subject to a drug of a particular concentration.
In non-dimensional form, my equations become:
Alpha 1-4 are unknown, but alpha 5 is known from my parameter estimation of the no drug model (i.e. logistic growth) which gives me g (growth rate, which I use for my ND timescale) and Csbar – which are given in the code.
Below is the main function that calls the ode's (note K in the code is Csbar in the model equations above). The fminsearch, and both ode solvers are presented as attachments.
function erri=ExpData(params)format longglobal g K alp1 alp2 alp3 alp4 alp5%These values (g and K) were determined from the no drug model.
%g will determin the non-dimensional timescale
%K is the carrying capacity - used to scale back to dimension level of S
%These are both estimated from the no drug case.
g = 0.54007;K = 85132.4387;%Parameters estimated in this model.
alp1=params(1);alp2=params(2);alp3=params(3);alp4=params(4);% alp5=params(5);
%Dimensional time in days
texp=[0 3 6 9 12]; %Experimental time in days
%Non-dimensional time-using time in days - multiplying by g - the timescale
texpND = g*texp;%Med Drug
%Divided both by K to get into ND form (CND = C/K) - where ND means
%non-dimensional
CN_DLL = [10640.3 22684.2 37263.6 60832.2 68035]/K;errLL = [1556.5 1613.7 2765.4 1555.4 3573.9]/K;% Define Species IC - first entry is species, and second entry is for bound
% drug level
x0 = [10640.3/K 0];%----------------------
%Splitting up times into pre and post wash
%Drug only in contact for 60mins (1/24).
tpw1 = linspace(0,1/24,1000);tpw2 = linspace(1/24,12,1000);%Converting time to non-dimensional
tpw1ND = tpw1*g;tpw2ND = tpw2*g;%Ode Solver - Main work - Prewash.
[t,xpw1] = ode45(@MedDrugODEPreWash,tpw1ND,x0);%Extracting final value which will then be inputted as an IC into the post
%wash case.
zpw0 = [xpw1(length(tpw1ND),1) xpw1(length(tpw1ND),2)];[t,xpw2] = ode45(@MedDrugODEPostWash,tpw2ND,zpw0);%Extracting the SMC data:
zi = xpw1(:,1);zi=transpose(zi);zl = xpw2(:,1);zl = transpose(zl);%Extracting the bound drug data:
bi = xpw1(:,2);bi=transpose(bi);bl = xpw2(:,2);bl = transpose(bl);erri=0;%errl=0;
for j=1:length(texpND) erri= erri +((zi(j)+zl(j))- CN_DLL(j))^2; % errl=errl+(zl(j)-CN_DLL(j))^2;
%errT=erri+errl;
enderrifigure(1)yyaxis leftscatter(texpND, CN_DLL, 'filled', 'g', 'Linewidth',2.5)xlabel('Time (ND)'), ylabel('Normalized cell count'), title('Med Drug (1mM/L) Estimation');set(gca,'fontsize',24)hold on% scatter(t,CN_D1, 'filled', 'r', 'Linewidth', 2.5)
% hold on
plot(tpw1ND,zi, 'b--','LineWidth',2.5)plot(tpw2ND,zl, 'g--','LineWidth',2.5)% errorbar(t, CN_D1, err, 'r', 'LineStyle','none','Linewidth',1.5);
% hold onerrorbar(texpND, CN_DLL, errLL, 'g', 'LineStyle','none','Linewidth',1.5);% yyaxis right
% plot(tpw1ND,bi, 'r--','LineWidth',2.5)
% plot(tpw2ND,bl, 'k--','LineWidth',2.5)
legend ('Experimetnal Data - Medium Drug','Ode Solver - SMC Pre wash','Ode Solver - SMC Post wash','Error bars','Ode Solver - b Pre wash','Ode Solver - b Post wash')%title('Release Profile', 'FontSize',28,'FontName','Times New Roman');
%title([ 'da0=', num2str(da0),', T=', num2str(T), ', s=', num2str(s), ', gam=', num2str(gam), ', diff=', num2str(diff),'err=', num2str(err)], 'fontsize',24)
drawnow
Interestingly about this code, is that I have 2 ode solvers for the above equations. For the first 60 minutes, drug is present and the equations are as shown above. However, after this time drug is removed and thus alp3 = 0.
As mentioned above, maybe the model just won't capture the data due to the growth rate attained from the no drug model being too steep.
Any help, or clarity would be greatly appreciated.
Note: I also tried an additional set of equations, where the Eq. 1 was multiplied by some constant alp5 that allowed for the model simulation result to pass through the data points quite nicely (i.e. I know it fits, but may not be optimal). So, then when I run the fminsearch code again, the curve moves away from the data points and towards the y – axis, yet the error value decreases.
This is why I think I may have messed up somewhere in my code/error calculation.
If anything is unclear, I would be happy to try solve these issues in the comments.
Thank you.
Best Answer