Fitting system of ODEs to data - multistart?

4 次查看(过去 30 天)
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

采纳的回答

Mary Fenelon
Mary Fenelon 2019-6-19
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 个评论
Matt Bilyard
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!).

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Global or Multiple Starting Point Search 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by