MATLAB: Fitting system of ODEs to data – multistart

Question

Hi,

I am a researcher in chemical biology who has very recently started using Matlab – as such, I am very much an amateur at it! This is also my first post on here so I apologise if any formatting etc. fails and so forth – please let me know if so.

I have a system of rate equations, as ODEs, that model conversion of a metabolite to its various oxidised derivatives. I'd like to fit experimental data to these ODEs to try and get values for the 7 unknown rate constants, "k(1) to k(7)", in each equation (realistically, semi-quantitative, i.e. approximate relative magnitudes is fine).

So far I've been using lsqcurvefit to do this, but the values and fit I get out depend very much in what initial values I choose for the rate constants k(1) to k(7), presumably since multiple local minima exist. I've copied the code for this below. The fit also isn't that good for the lower-magnitude data (though this could be the model itself, of course) – see attached images of an example.

I've heard that I might use multistart or globalsearch in order to improve this, i.e. find a global minimum. However, despite extensive reading around, including the various tutorials, I am really stuck as to how practically to do this for this fairly complex system!

I wonder if anyone could:

1. give me an idea if this is correct, and perhaps which of multistart and globalsearch I might be best starting with?
2. point me in the right direction as to how I'd begin setting up multistart or globalsearch for this system, based on the code for the lsqcurvefit fitting below? I'm more than happy to figure out things independently, but a push along the right lines to start me off would be really appreciated given my total inexperience in this area.

Thanks a lot,

Matt

Code for the lsqcurvefit fitting process (includes data to fit):

function cytosinefitcell
 
 
    function C=kinetics(k,t)
        c0=[95.8989;0;0;0;0];
        [T,Cv]=ode45(@DifEq,t,c0);
        
        function dC=DifEq(t,c)
            dcdt=zeros(5,1);
            
            dcdt(1)= -k(1).*c(1)+k(5).*c(4)+k(6).*c(5)+k(7).*(c(2)+c(3)+c(4)+c(5));
            dcdt(2)= k(1).*c(1)-k(2).*c(2)-k(7).*c(2);
            dcdt(3)= k(2).*c(2)-k(3).*c(3)-k(7).*c(3);
            dcdt(4)= k(3).*c(3)-k(4).*c(4)-k(5).*c(4)-k(7).*c(4);
            dcdt(5)= k(4).*c(4)-k(6).*c(5)-k(7).*c(5);
            dcdt(1) = 0;
            dC=dcdt;
            %dcdt(1) = 0 as by definition this variable is constant, i.e. steady state.
            
            
        end
        C=Cv;
    end
 
%data for t
t=[0
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    24
    48
    96
    144
    192
    240
    288
    336
    384
    432];
 
%data for c(1) to c(5)
c=[95.8989  0   0   0   0
    95.8989 0.116654171 0.000169089 0   0
    95.8989 0.221311714 0.000215784 0   0
    95.8989 0.361815956 0.001337332 0   0
    95.8989 0.443043182 0.001897635 0   0
    95.8989 0.511783075 0.002904908 2.59348E-06 0
    95.8989 0.596086415 0.003847949 2.08165E-05 0
    95.8989 0.70422927  0.004991988 3.23739E-05 0
    95.8989 0.736165548 0.005402772 4.12232E-05 0
    95.8989 0.863634039 0.006882534 4.66973E-05 0
    95.8989 0.961691531 0.007922809 7.91253E-05 0
    95.8989 1.694849679 0.02488041  0.000229454 3.88541E-05
    95.8989 2.156329216 0.046290117 0.000455848 2.44297E-05
    95.8989 2.28066375  0.058965709 0.000529312 5.23961E-05
    95.8989 2.263872217 0.060281286 0.000604404 5.68658E-05
    95.8989 2.280749095 0.059985812 0.000559687 6.81934E-05
    95.8989 2.250775532 0.059775064 0.00057279  6.40691E-05
    95.8989 2.248860841 0.059972783 0.000589628 5.50697E-05
    95.8989 2.28310359  0.059096809 0.000519199 5.7434E-05
    95.8989 2.324286746 0.058745232 0.000550764 5.6446E-05
    95.8989 2.298780141 0.059541016 0.000556106 6.05984E-05];
 
%initial guesses for k - arbitrarily set to 0 except for k(7) which has
%known range. Each "k" is really a "k_obs".
k0=[0;0;0;0;0;0;0.04];
 
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,t,c,[0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.04],[100;100;100;100;100;100;0.07]);
%upper and lower bounds are arbitrary to avoid 0 (or negative) values or
%unrealistically high values. Exception is k(7), whose range is known.
 
%plotting graph
fprintf(1,'\tRate Constants:\n')
for k1 = 1:7
    fprintf(1, '\t\tk(%d) = %8.5f\n', k1, k(k1))
end
 
tv = linspace(min(t), max(t));
Cfit = kinetics(k, tv);
 
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'mC', 'hmC', 'fC','caC', 'Location','NE')
 
end

Answer 1

This example should get you started. Information on how the solvers work can be found here. The setup is similar for both so you can easily try both.