simultaneouslty fIt ode soln to multiple datasets with lsqcurvefit

I have a data matrix, A, whose columns I want to simultaneously fit to the solution of the ode below using lsqcurvefit. Each column of A corresponds to a different value of L. For example, if

L=1:25
time=[0:0.001:3.8]'

The data matrix has the dimensions time x L. I have tried using meshgrid to represent the independent variables, but my attempts return the error:

"The entries in tspan must strictly increase or decrease."

So it seems that the function is reading matrix generated by meshgrid as a continuous time series.

Any guidance would be greatly appreciated.

function [C]=conformselect(K,time)
Y0=[1;1;0];
[tSol,YSol]=ode45(@diffeq,time,Y0);
function dYdt = diffeq(time,Y)
a=Y(1);
b=Y(2);
c=Y(3);
k1=K(1);
km1=K(2);
k2=K(3);
km2=K(4);
dadt=-k1*a+km1*b;
dbdt=k1*a-(km1+k2*L)*b+km2*c;
dcdt=k2*L*b-km2*c;
dYdt=[dadt;dbdt;dcdt];
end
A=YSol(:,1);
C=YSol(:,3);
end

Answer 1

[Star Strider]

Your function appears to be essentially correct.

The lsqcurvefit function can fit matrix dependent variables. However you need to return all the columns you want to fit in your ‘C’ output. You are currently returning only one column, ‘YSol(:,3)’.

For various approaches to this sort of problem, see:

• Monod kinetics and curve fitting (link)
• Solving Coupled Differential Equations (link)
• Parameter Estimation for a System of Differential Equations (link)

[John Hackett]

Thanks for taking the time to look at my problem. If A is just a single column vector, this works like the posted examples. In this problem, the columns of A do not represent different species, rather the time dependence of one species, c, at different values of L. This function generates some data with random noise:

function [time,ligand,conc_sim]=conformselect_sim
time=[0:0.01:4]';
ligand=[1:10]';
K=[1 2 10 2]
Y0=[1 1 0]
for i=1:length(ligand)
    L=ligand(i)
    conc_sim(:,i)=conformselect(K,time)+0.05*rand(length(time),1)
end
function [C]=conformselect(K,time)
[tSol,YSol]=ode113(@diffeq,time,Y0);
function dYdt = diffeq(time,Y)
a=Y(1);
b=Y(2);
c=Y(3);
k1=K(1);
km1=K(2);
k2=K(3);
km2=K(4);
dadt=-k1*a+km1*b;
dbdt=k1*a-(km1+k2*L)*b+km2*c;
dcdt=k2*L*b-km2*c;
dYdt=[dadt;dbdt;dcdt];
end
C=YSol(:,3);
end
end

I may be on the wrong track, but I think I need to be passing the independent variables generated from

[L T]=meshgrid(ligand,time)

to the conformselect function above. I know how to pass column vectors to conformselect, though I don't know how to modify it to accept matrices of simulated data and the independent variable matrices L and T. I have seen a number of errors, but it always seems to get hung up in the ode solver. What am I missing?

[Star Strider]

I can’t run your code since I don’t have the arguments.

If it’s getting ‘hung up’ on the ODE solver call, see if changing to ode15s improves your result. Most kinetics problems are ‘stiff’, creating serious problems for the non-stiff solvers.

I don’t believe meshgrid has any current application to your problem.

Since you’re fitting matrix dependent variables, you simply need to return a matching matrix of values from your objective function (that contains the ODE integration code) so that lsqcurvefit can use them to tune the parameters. In your current code, you’re returning one column in ‘C’. That will throw an error if you use it with a matrix of dependent variables in lsqcurvefit.

[John Hackett]

I wish I had made this much progress. The conformselect function will not take matrix variable and I don't know how to modify it to do so. For example, if I pass a time matrix to conformselect:

ligand=[1:10]';
time=[0:0.1:4]'*ones(1,length(ligand))
K=[1 2 10 2]
Csim=conformselect(K,time,ligand)

it returns:

Error using odearguments (line 87) The entries in tspan must strictly increase or decrease.

I know that I am only returning one column of C, but the code doesn't make it that far. The ode solver reads the time matrix as a continuous series.

[Star Strider]

I have no idea what ‘conformselect’ even does. This limits my ability to help you with it.

This will create a (41x10) matrix of identical columns of ‘time’:

time=[0:0.1:4]'*ones(1,length(ligand));

If you’re using ‘time’ as the ‘tspan’ variable, define it as:

time=[0:0.1:4]';

If you want to multiply it by the length of the ‘ligand’ variable, just do:

time = [0:0.1:4]*length(ligand);

or to multiply it by the last value of ‘ligand’ (that here will produce the same result):

time = [0:0.1:4]*ligand(end);

If you want ‘time’ to be the same length as ‘ligand’, use the linspace function:

time = linspace(0, 4, length(ligand));

[John Hackett]

I think we are going backwards-I'll repost code to avoid confusion. Let me start with what I want to build on. The following code calls conformselect with a single value for L, a column vector with time, and a four-element rate constant vector K. The output is the concentration of species C:

function [time,conc_sim]=simulate_data
%%Input parameters
L=1
time=[0:0.1:4]
K=[1 2 10 2]

%%Simulate Data
conc_sim=conformselect(K,time)

%%Function conformselect
function [C]=conformselect(K,time)
Y0=[1 1 0]
[tSol,YSol]=ode15s(@diffeq,time,Y0);
function dYdt = diffeq(time,Y)
a=Y(1);
b=Y(2);
c=Y(3);
k1=K(1);
km1=K(2);
k2=K(3);
km2=K(4);
dadt=-k1*a+km1*b;
dbdt=k1*a-(km1+k2*L)*b+km2*c;
dcdt=k2*L*b-km2*c;
dYdt=[dadt;dbdt;dcdt];
end
C=YSol(:,3);
end
end

This works fine and I can easily extract the nested function conformselect to fit data using lsqcurvefit. Indeed, ode15s appears to be the most efficient.

Heres the first problem. I want to modify conformselect this so that it takes a vector L, i.e

L=1:10

to produce a matrix conc_sim that has length(time) rows and length(L) columns. The columns of the matrix conc_sim will be the concentration of species C for each value of L.

I understand that this could be achieved by adding a loop over L, but I don't see how this will be conducive to later integration into lsqcurvefit.

[John Hackett]

Sorry about all of the confusion and headache-the answer was in front of me (the loop over L). Not elegant but it works:

%%Define Ligand Concentrations,Rate Constants,and Time Span
ligand=[0.5 1 2 4 6 10 15 10];
K0=[1 2 10 2];
dt=0.001
time=[0:dt:3.8]'
error=0.05
a0=1
global ligand a0
%%Simulate Data With Error
conc_sim=conformselect2(K0,time,ligand)+error.*rand(length(time),1);

%%Fit Simulated Data
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt','Display','iter-detailed');
[Kfit,resnorm,residual,exitflag]=lsqcurvefit(@conformselect2,rand(1,4),time,conc_sim,zeros(4),Inf*ones(4),options);
%%Comparison of initial and fitted rate constants
K0
Kfit
%%Plot the data
plot(time,conc_sim,'.','MarkerSize',1);hold on
plot(time,conformselect2(Kfit,time),'k-','Linewidth',1)
xlim([0 max(time)]);xlabel('time(s)');ylabel('Concentration')

function [conc_sim]=conformselect2(K,time,ligand)
global ligand
for i=1:length(ligand)
    L=ligand(i);
conc_sim(:,i)=conformselect(K,time,L);
end

function [C]=conformselect(K,time,L)
Y0=[1 K(1)/K(2) 0];
[tSol,YSol]=ode15s(@diffeq,time,Y0);
function dYdt = diffeq(time,Y)
a=Y(1);
b=Y(2);
c=Y(3);
k1=K(1);
km1=K(2);
k2=K(3);
km2=K(4);
dadt=-k1*a+km1*b;
dbdt=k1*a-(km1+k2*L)*b+km2*c;
dcdt=k2*L*b-km2*c;
dYdt=[dadt;dbdt;dcdt];
end
C=YSol(:,3);
end
end

[Star Strider]

I’m lost. I still don’t understand.

So problem solved?

[John Hackett]

Yes, thanks for your help. As you mentioned, the function needed to return all of the columns of C. I was distracted by looking for an alternative way to achieve that without a loop, but it turned out to work fine.

[Star Strider]

My pleasure.

If my Answer helped you solve your problem, please Accept it!