Clear Filters
Clear Filters

Fitting a series of coupled ODEs to a dataset with fminsearch

5 views (last 30 days)
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 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.
% 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
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);
zl = xpw2(:,1);
zl = transpose(zl);
%Extracting the bound drug data:
bi = xpw1(:,2);
bl = xpw2(:,2);
bl = transpose(bl);
for j=1:length(texpND)
erri= erri +((zi(j)+zl(j))- CN_DLL(j))^2;
% errl=errl+(zl(j)-CN_DLL(j))^2;
yyaxis left
scatter(texpND, CN_DLL, 'filled', 'g', 'Linewidth',2.5)
xlabel('Time (ND)'), ylabel('Normalized cell count'), title('Med Drug (1mM/L) Estimation');
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)
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.

Accepted Answer

Star Strider
Star Strider on 12 Jul 2019
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.
Star Strider
Star Strider on 6 Feb 2020
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.)

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!