Fitting system of ODEs to data using lsqcurvefit and ode45

I am fitting a system of 3 differential equations that model change in concentrations of molecules over time to experimental data. I'd like to optimize four parameters in the third equation (theta(1), theta(2), theta(3), theta(4)). I am receiving an error message that the initial point is a local minimum. I've tried changing the initial theta values with no luck, as well as changing optimality and function tolerance. Additionally, the last 3 rows of output vector C are NaN and I'm not sure why (this is why I included the line C(find(isnan(C))) = 0;). Any help would be greatly appreciated!

function C = kinetics(theta,t)
    c0 = [0 0 0];     %[AA_conc BB_conc MKATE2]
    [T,Cv]=ode45(@DifEq,t,c0);

    function dC = DifEq(t,c)
        dcdt=zeros(3,1);
        
        AA_conc = c0(1); 
        BB_conc = c0(2); 
        MKATE2 = c0(3); 
        
        nH_Est = 1.2;              % Hill for Est known
        nH_Ald = 1.2;              % Hill for Ald known 
        nH_het = 1.2;              %  

        kact = 1;                  % Transcriptional activation rate for all equations

        kdeg1 = 0.01;              % Degredation rate of AA
        kdeg2 = 0.01; 
        
        ec501 = 64; %nM, known (Est)
        ec502 = 200; %nM, known (Ald)
        
        tmin = 0;                    % start time
        tmax = 900;                 % end time (min.)
        inc  = 17;                  % # of time steps
        time = linspace(tmin,tmax,inc); 
        
        ON_time = 900;              % Units = min.
        OFF_time = 0;             % Units = min.
        Est  = 100*[ones(1,(ON_time/tmax)*inc) zeros(1,(OFF_time/tmax)*inc)];
        
        ON_A = 900; %ON for aldosterone, units = min 
        OFF_A = 0; %OFF for aldosterone, units = min 
        Ald = 100*[ones(1,(ON_A/tmax)*inc) zeros(1,(OFF_A/tmax)*inc)]; %units = nM
               
        EstI = interp1(time,Est,t)';
        AldI = interp1(time,Ald,t)';
        
        dcdt(1) = kact.*[EstI].^nH_Est./([EstI].^nH_Est + ec501.^nH_Est) - kdeg1.*c(1);
        dcdt(2) = kact.*[AldI].^nH_Ald./([AldI].^nH_Ald + ec502.^nH_Ald) - kdeg2.*c(2);
        dcdt(3) = kact.*[c(1).*c(2)./(theta(3)).^0.5]^theta(4)./([c(1).*c(2)./(theta(3)).^0.5].^theta(4) + theta(2).^nH_het) - theta(1).*c(3);
        dC=dcdt;
    end
C=Cv(:,3);
C(find(isnan(C))) = 0;

end

To read in the experimental data and call the function:

%% read in experimental data  
simVal = readtable('08Dec21_Table.xlsx'); 
sim330 = simVal.yAS330;
sim330 = sim330'; 
simTime = [15,18:32,45];
simTime = simTime.*60;

%% fitting  

sim330 = 100*sim330./max(sim330);
t=simTime';
c=sim330';
theta0=[kdeg3;ec503;kobs;nH_het];

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

The function currently returns the initial theta parameters without optimizing. Below is a plot of the solution to the third differential equation (using initial theta values) overlaid with the experimental data. I can provide additional information as needed. Thanks!!

[Star Strider]

Don’t do this:

c0 = [0 0 0]; %[AA_conc BB_conc MKATE2]

Estimate the initial conditions with the other parameters, appending them as the last elements in the parameter vector. Change the code in ‘Kinetics’ accordingly to accommodate that.

[Hannah Collins]

Okay, I appended to the parameter vector and updated the 'kinetics' code. I now receive the warning message "Local minimum possible. lsqcurvefit stopped because the size of the current step is less than the value of the step size tolerance." I've also tried changing maxfunctionevaluations, steptolerance, functiontolerance, optimalitytolerance, and maxiterations. Thanks!

Parameter vector:

theta0=[kdeg3;ec503;kobs;nH_het;AA_conc;BB_conc;MKATE2];

Kinetics function:

function C = kinetics(theta,t)
    %c0 = [0 0 0];     %[AA_conc BB_conc MKATE2]
    dcdt=zeros(3,1);
    
    [T,Cv]=ode45(@DifEq,t,dcdt);

    function dC = DifEq(t,dcdt)
        
        %AA_conc = c0(1); 
        %BB_conc = c0(2); 
        %MKATE2 = c0(3); 
        
        nH_Est = 1.2;              % Hill for Est known
        nH_Ald = 1.2;              % Hill for Ald known
        nH_het = 1.2;              % Optimize this value

        kact = 1;                  % Transcriptional activation rate for all equations

        kdeg1 = 0.01;              % Degredation rate of AA
        kdeg2 = 0.01; 
        
        ec501 = 64; %nM, known (Est)
        ec502 = 200; %nM, known (Ald)
        
        tmin = 0;                    % start time
        tmax = 900;                 % end time (min.)
        inc  = 17;                  % # of time steps
        time = linspace(tmin,tmax,inc); 
        
        ON_time = 900;              % Units = min.
        OFF_time = 0;             % Units = min.
        Est  = 100*[ones(1,(ON_time/tmax)*inc) zeros(1,(OFF_time/tmax)*inc)];
        
        ON_A = 900; %ON for aldosterone, units = min 
        OFF_A = 0; %OFF for aldosterone, units = min 
        Ald = 100*[ones(1,(ON_A/tmax)*inc) zeros(1,(OFF_A/tmax)*inc)]; %units = nM
               
        EstI = interp1(time,Est,t)';
        AldI = interp1(time,Ald,t)';
        
        dcdt(1) = kact.*[EstI].^nH_Est./([EstI].^nH_Est + ec501.^nH_Est) - kdeg1.*theta(5);
        dcdt(2) = kact.*[AldI].^nH_Ald./([AldI].^nH_Ald + ec502.^nH_Ald) - kdeg2.*theta(6);
        dcdt(3) = kact.*[theta(5).*theta(6)./(theta(3)).^0.5]^theta(4)./([theta(5).*theta(6)./(theta(3)).^0.5].^theta(4) + theta(2).^theta(4)) - theta(1).*theta(7);
        dC=dcdt;
    end
C=Cv(:,3);
end

[Star Strider]

In ‘Kinetics’ the ‘c0’ vector needs to be:

c0 = theta(5:7);

Append three non-zero values (they can be near zero) to the existing ‘theta0’ vector and go from there.

You coded this correctly the first time:

[T,Cv]=ode45(@DifEq,t,c0);

then changed it in the Comment.

There are more sources of relevant code in this search.

[Hannah Collins]

Thanks for the clarification. The code seems to be now trying to optimize the initial conditions c0, but not really changing the constant values for theta(1:4). lsqcurvefit is still stopping prematurely, but returns the following output values:

theta =

    0.0099
  131.5735
    0.0000
    1.1992
    0.1604
    0.1593
   10.0000

I set bounds for theta and it's unclear to me why theta(7), representing c0 for a molecule that is known to have initial concentration of 0, has a value of 10 here (upper bound limit).

ub = [1;500;1;2;10;10;10];
lb = [0;0;0;1;0;0;0];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,lb,ub,[]);

For reference, these are the three differential equations that I'm using. I only want to optimize the kobs, nH_het, ec503, and kdeg3 constants in the third equation.

% d[AA]/dt
dxdt(1) = kact*[EstI]^nH_Est/([EstI]^nH_Est + ec501^nH_Est) - kdeg1*[AA_conc]; 

% d[BB]/dt
dxdt(2) = kact*[AldI]^nH_Ald/([AldI]^nH_Ald + ec502^nH_Ald) - kdeg2*[BB_conc];

%d[mKate]/dt                    (AB_conc = (AA_conc.*BB_conc./kobs).^0.5)
dxdt(3) = kact*[(AA_conc.*BB_conc./kobs).^0.5]^nH_het/([(AA_conc.*BB_conc./kobs).^0.5]^nH_het + ec503^nH_het) - kdeg3*[MKATE2];

The current output plot is attached. Thanks!

[Star Strider]

Without the data, I can¹t determine what the problem is. Sometimes it’s just necessary to keep experimenting with different initial parameter estimates to see what works best, the reason I generally use random numbers for the initial parameter estimates (scaled appropriately). If you have the Global Optimization Toolbox, I’ll post a ga (genetic algorithm) approach to this problem that generally works. Otherwise, keep trying different initial parameter estimates until one provides a good fit. That’s simply the way nonlilnear parameter estimation works.