Fitting system of ODEs to data - multistart?

9 visualizaciones (últimos 30 días)
Matt Bilyard
Matt Bilyard el 18 de Jun. de 2019
Editada: Jan el 25 de Jun. de 2019
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

Respuesta aceptada

Mary Fenelon
Mary Fenelon el 19 de Jun. de 2019
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 comentario
Matt Bilyard
Matt Bilyard el 21 de Jun. de 2019
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!).

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Global or Multiple Starting Point Search en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by