Fitting a series of coupled ODEs to a dataset with fminsearch
3 次查看(过去 30 天)
显示 更早的评论
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/229522/image.png)
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 long
global 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;
end
erri
figure(1)
yyaxis left
scatter(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 on
errorbar(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.
0 个评论
采纳的回答
Star Strider
2019-7-12
First, please do not use global variables! Pass your data as parameters instead. That way, you know what you’re giving your functions. See Passing Extra Parameters for a full discussion.
Second, the ‘MedDrug_fitting’ function (that may be important) is not available (at least I can’t find it).
Third, let me introduce you to Monod kinetics and curve fitting. It is a prototype for fitting fitting nonnlinear differential equations to data. Give that approach a try with your differential equations and data.
10 个评论
Star Strider
2020-2-6
As always, my pleasure!
I don’t use fmincon often, and there are others here with more recent experience with it. (I will look for your new Question and see if I can solve the constrained problem.)
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Least Squares 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!