MATLAB Answers


How can I parameter estimate many unknown parameters in a system of nonlinear ODE?

Asked by george
on 8 Aug 2018 at 20:07
Latest activity Commented on by george
on 16 Aug 2018 at 17:53


I have been following the links below:


I find them helpful, but when I implement the method to my own system of ODE. I am generating wrong answers with the initial guesses I am using from a research manuscript I am following (the parameters used to refit the ode model does not match well with the exp data graphs). So I decided to randomly guess the initial parameter guess using rand. But the code takes over 2 hours to run, so I usually terminate it before the code actually finishes. I tried to use MultiStart to get better parameter initial guesses, but I am not sure if it is working. Am I doing something wrong with the code? Would it be best to use something else other than lsqcurvefit?

I have attached the excel file for the data and the code itself. Many thanks if this problem can be fixed, I have been stuck for weeks now!


Sign in to comment.

1 Answer

Answer by Arthur Goldsipe on 14 Aug 2018 at 12:48

I see that you mentioned SimBiology in a comment on the previous answer. I'm a developer on that product, so I wanted to share that it is designed specifically for this sort of problem. You can build the model graphically or programmatically with building blocks like reactions that more closely map to the underlying biology. (Not to mention the fact that this can save you a lot of bookkeeping.) You can also give variables more meaningful names, instead of having to figure out what y(6)/c(38) means.

There are also a lot of performance benefits. The way we do the ODE simulations in SimBiology is also typically quite a bit faster than doing the simulations in core MATLAB yourself. My rule of thumb is that we usually give at least a 20X speedup.

I personally have also spent a lot of time trying to make it easier and faster to do "large" fitting problems exactly like this. When using gradient-based optimization methods like lsqnonlin, we can often provide more accurate estimates of gradients than the default. This can help with the convergence of the algorithm. I have also worked on a couple of global optimization methods that make it easier to scan a large parameter space. Most recently, I implemented an enhanced scatter search algorithm that combines both global and local optimization techniques.

If this sounds interesting to you, I can take a look to see how much work I think it would be to convert your problem over to SimBiology.


I guess I will give SimBiology a try since the other posted answer was deleted. I don't know too much about SimBiology but I tried to look up references in the past to recreate the model reactions.

Here is a little snippet of what I tried to do so far:

if true
  % code
% Creating a CimBiology model
m1 = sbiomodel('m1');
% Add compartment, 13 species, 22 reaction 
c1 = addcompartment(m1, 'cell');
s1 = addspecies(m1, 'glcext');
s2 = addspecies(m1, 'glnext');
s3 = addspecies(m1, 'Xv');
s4 = addspecies(m1, 'mAb');
s5 = addspecies(m1, 'UDPGlcNAc');
s6 = addspecies(m1, 'UDPGalNAc');
s7 = addspecies(m1, 'GDPMan');
s8 = addspecies(m1, 'GDPFuc');
s9 = addspecies(m1, 'UDPGlc');
s10 = addspecies(m1, 'UDPGal');
s11 = addspecies(m1, 'CMPSa');
s12 = addspecies(m1, 'glcint');
s13 = addspecies(m1, 'glnint');
*r1 = addreaction(m1, 'glcext + glnext -> Xv'); % Can account for growth and death of it* 
*r2 = addreaction(m1, 'Xv -> glcext');
*r3 = addreaction(m1, 'Xv -> glnext'); 
*r4 = addreaction(m1, 'Xv -> mAb');* 
r5 = addreaction(m1, 'UDPGlc -> UDPGal');
r6 = addreaction(m1, 'GDPMan -> GDPFuc');
r7 = addreaction(m1, 'glcint + glnint -> UDPGlcNAc');
r8 = addreaction(m1, 'glcint -> UDPGlc');
r9 = addreaction(m1, 'glcint -> GDPMan');
r10 = addreaction(m1, 'UDPGlcNAc -> UDPGalNAc'); 
r11 = addreaction(m1, 'UDPGlcNAc -> CMPSA');
r12 = addreaction(m1, 'UDPGlcNAc -> null', 'Name', 'UDPGlcNAc degradation');
r13 = addreaction(m1, 'CMPSA -> null', 'Name', 'CMPSA degradation');
r14 = addreaction(m1, 'UDPGal -> null', 'Name', 'UDPGal degradation');
r15 = addreaction(m1, 'GDPFuc -> null', 'Name', 'GDPFuc degradation');
r16 = addreaction(m1, 'UDPGlc -> null', 'Name', 'UDPGlc degradation');
r17 = addreaction(m1, 'UDPGalNAc -> null', 'Name', 'UDPGalNAc degradation');
r18 = addreaction(m1, 'GDPMan -> null', 'Name', 'GDPMan degradation');
r19 = addreaction(m1, 'glcext -> glcint'); % 'MassAction' 
r20 = addreaction(m1, 'glnext -> glnint'); % 'MassAction'
*r21 = addreaction(m1, 'glcint -> null', 'Name', 'glcint metabolite consumption'); % 'Henri-Michaelis-Menten'*
*r22 = addreaction(m1, 'glnint -> null', 'Name', 'glnint metabolite consumption'); % 'Henri-Michaelis-Menten'*
% Undefined Reaction rate since competitive inhibition reaction 3
% r3.ReactionRate = 'Vm/(Kh/a)^n + 1';
% k3.KineticLaw = 'Unknown';
% k3.Expression 
% predefined a michaelis menten for reaction 2 ('Henri-Michaelis-Menten')
k2 = addkineticlaw(r2, 'Henri-Michaelis-Menten');
k2.Expression % Display the built in kinetic law 
% Predefined michaelis menten for reaction 5 ('Henri-Michaelis-Menten')
k5 = addkineticlaw(r5, 'Henri-Michaelis-Menten');
k5. Expression
% Undefined Reaction rate
r6.ReactionRate = 'VmaxGDPFuc*GDPMan/(KmGDPFuc*(1+GDPMan/KmGDPFuc+GDPMan/KmGDPFuc*GDPFuc/KiGDPFuc+GDPFuc/KiGDPFuc))';
k6.KineticLaw = 'Uknown';
% Manual eqn into simbiology *********
% Undefined Reaction rate
r7.ReactionRate = 'VmaxUDPGlcNAc *(GlcInt/(KUDPGlcNAc_glc+GlcInt))*(GlnInt/(KUDPGlcNAC_gln+GlnInt))';
k7.KineticLaw = 'Unknown';
% Predefined michaelis menten for reaction 8 ('Henri-Michaelis-Menten')
k8 = addkineticlaw(r8, 'Henri-Michaelis-Menten');
k8. Expression
% Predefined michaelis menten for reaction 9 ('Henri-Michaelis-Menten')
k9 = addkineticlaw(r9, 'Henri-Michaelis-Menten');
k9. Expression
% Predefined michaelis menten for reaction 10 ('Henri-Michaelis-Menten')
k10 = addkineticlaw(r10, 'Henri-Michaelis-Menten');
k10. Expression
% Undefined Reaction rate
r11.ReactionRate = 'VmaxCMPSA *(UDPGlcNAc)^0.45/((KCMPSA_UDPGlcNAc +UDPGlcNAc)^0.45*(1/(1+CMPSA/KiCMPSA))^4.1)';
k11.KineticLaw = 'Unknown'; 

Would you have any recommendations, the code above does not cover any extracellular activities (I am a little confused on how SimBiology can model that part as to how those reactions (r1, r2, r3, r4 - I am not sure I added in the reactions correctly either) are modeled in the toolbox)? I have attached the equations revolving the extracellular activity to be modeled. Would it be easier to write a code like the above or use building blocks in SimBiology to recreate the model for parameter estimation? Thank you for your time!

Sorry for the delayed reply. If I understand your question, the easiest way to incorporate extracellular activities is to introduce a compartment in the model that represents the extracellular space. You can create reactions that involve species in both your extracellular and cellular space. SimBiology will do the correct bookkeeping for conservation of mass, as long as you use the correct stoichiometric coefficients, volumes, and units on the model.

One way for you to check that you've set up your model correctly is to view the ODE equations for it. You can do that by calling getequations(m1).

It's difficult for me to say whether you have set up any specific reaction correctly. I would need to understand the details of how you intend to model your system. I know you've tried to provide references to address this, but I'm afraid I haven't had time to read the entire PDF you attached to an earlier comment. One way you could help is to show me specific terms from specific differential equations, and I can help you confirm that you've translated them correctly into reactions.

I also want to mention another alternative. If it's too difficult for you to reverse-engineer the reactions for your system, you can also directly enter your differential equations into SimBiology using rate rules . For example, let's consider the simplest of systems, degradation of a species with mass action kinetics. You could model it in the following two ways:

    % Create the model
    model = sbiomodel('m1');
    % I recommend turning on unit conversion
    configset = getconfigset(model);
    configset.CompileOptions.UnitConversion = true;
    % Create the compartment and species
    c = addcompartment(model, 'cell', 100, 'CapacityUnits', 'micrometer^3');
    addparameter(model, 'k', 0.1, 'ValueUnits', '1/second');
    addspecies(c, 's1', 10, 'InitialAmountUnits', 'millimole/liter');
    % Model the reaction for s1 using a reaction object
    reaction = addreaction(model, 's1 -> null');
    kl1 = addkineticlaw(reaction, 'MassAction');
    kl1.ParameterVariableNames = {'k'};
    % Model the reaction for s2 using a rate rule
    addspecies(c, 's2', 10, 'InitialAmountUnits', 'millimole/liter');
    rule = addrule(model, 's2 = -k*s2', 'rate');
    % Simulate the model and plot the results
    sd = sbiosimulate(model);

If you're still stuck, let me know. It's best if you can make your question as specific as possible.


Thank you for following up!

I will try the suggestion above and directly enter my system of ode into SimBiology using rate rules. I have attached an image of the reaction rates for extracellular concentration in an older comment called "Sample Reaction eqns." Let me know if you need that PDF from an earlier comment.

Thanks again. I will try that rate rule method.

Sign in to comment.