Fitting system of ODEs to data - multistart?

Question

Matt Bilyard 2019-6-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/467710-fitting-system-of-odes-to-data-multistart

编辑： Jan 2019-6-25

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:

give me an idea if this is correct, and perhaps which of multistart and globalsearch I might be best starting with?
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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mary Fenelon 2019-6-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/467710-fitting-system-of-odes-to-data-multistart#answer_379841

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Matt Bilyard 2019-6-21

Thanks very much for pointing me in this direction - it took me a while to work out how to adapt this for my specific problem but I think I'm getting somewhere now. (This is very much a side project alongside an unrelated main project for me, hence progress is slow!).

请先登录，再进行评论。