Fitting a series of coupled ODEs to a dataset with fminsearch

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 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.

Answer 1

[Star Strider]

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.

[Alistair McQueen]

Thanks for such a swift response!

Yeah you're right about the global parameters, just a nasty habbit from my undergrad courses...

The second point (my fault for lack of transparency) is the code presented in the question (i.e. not the file) is called 'MedDrug_fitting' But, to prevent any confusion I have attached it here regardless.

For the third point, (which is why I'm super grateful you responded) have been scanning through the forums on related problems the past couple of days and hence came across you several times helping other people.

I have tried applying these methods to my problem, but cannot seem to get them to work for my case. Maybe I am making a super trivial error... I am use to coding differential equations etc. but I often know the parameters, in this case I don't know the parameters and hence, trying to find the best estimate - so all the curve fitting stuff is rather new to me. I have looked at the link above and a couple others on solving for coupled differential equations that you presented in another answer, somwhere.

[Star Strider]

My pleasure.

O.K. I thought that since the function is ‘ExpData’, so is the file.

‘I have tried applying these methods to my problem, but cannot seem to get them to work for my case.’

It is not easy to follow your code.

If your data are in the ‘CN_DLL’ and ‘errLL’ vectors, you have 0nly 5 data pairs. Since is given (not to be estimated), you are then estimating 4 parameters. If your data are at all noisy, that is going to be extremely difficult for any optimisation routine to fit. (I plotted them, and it seems that they are essentially linear.)

‘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.’

It’s not obvious what your data represent in that context, since ‘texpND’ only goes from 1 to 7.

It is very difficult to follow your code. See if you can get the ‘Monod kinetics’ approach working with it with the differential equations system and data you want to fit. I can help if you have problems with it, however ar present I have no idea what your code does or what data you are fitting and to what set of differential equations you are fitting it.

[Alistair McQueen]

I understand.

I will try the approach mentioned again and see if I can get anything to work, and at the very least create a more sensible code documenting what I am trying to achieve.

I am away in Milan over the next week, but I will try my best when there to do this task also.

Hope to get something back to you as fast as I can. Once again, thank you.

[Star Strider]

My pleasure.

I can be of more help when your problem is more easily understandable.

[Alistair McQueen]

I managed to resolve the issue I was having, just took a bit of thought. Had issues with the time vectors due the washing out of the drug which occurs early doors. I was trying to calculate an error for time-points which never actually existed, which is rather naive of me... I got there though (with a few assumptions),and the fit is okay (attached)!

Thanks for you replies StarStrider, I appreciate your swift responses. I have accepted your answer now.

I am still having a few problems with my equations, but thats a more mathematical problem than a fitting one that will require more thought from myself.

[Star Strider]

As always, my pleasure.

Thank you.

[Alistair McQueen]

Hi Star Strider, I was curious if you were able to help me with another problem:

Previously I have been using fminsearch to fit my simulated data to experimental data for my PhD project. The results have been good up until now.

I have found a new dataset which requires a constraint to be added. In other words, to ensure that my model behaves as necessary, the parameters must not fall below some threshold value:

I have N = Nmax (1 - ((E1 + C)/(C + E2))).

Due to an initial condition, the above equation must not drop below that value. Therefore, the second term (fraction one) in the equation above must not fall below 0.9 (random value).

Thus, ((E1 + C)/(C + E2)) <= 0.9

Rearranging this, I can derive an expression for E1 in terms of E2.

Due to the constraints on the parameters here, fminsearch surely would not be applicable? I therefore tried to use fmincon, where I had A = [C -0.9] and b = 0.9*C (obtained from rearranging the above inequality - E1*C - 0.9E2 <= 0.9*C).

However, the results of this were rather odd.

I have attached an image (jpg format) if that may help. The first line (lowest one) is my initial guess actually fits the data fairly well. However, the fitting process continues, away from my data and then concludes at a local minimum, as the output suggests. What I don't understand is why, since the initial guess seems fairly good?

If need be, I can attach my code if it makes understanding the problem easier? I am new to using fmincon and before has only used fminsearch, but now I need the former since my parameters are bound by one another to be below some threshold.

[Star Strider]

I am not entirely sure that I understand what you want to do.

It would be best to post this as a new Question, and include your code (and data if possible).

[Alistair McQueen]

No problem, will do, thank you!

[Star Strider]

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.)