Using sbiofit to estimate/fit a complex model based on multiple measurments

Hi,

I have a rather complex enzyme-kinetic model. For didactic purposes I will just give a simplified version here. The problem should be eminent anyway. I'll also drop intermediate steps in the reaction.

E, F : Enzymes. A, B : Substrates. Q, R : Products.

E + F + A + B <-> E + F + Q + B <-> E + FQB <-> E + F + Q + R
            k21,k12           k32,k23     k34,k43

I did measure the concentration of Product "R". I altered the Substrate "A" in multiple series.

My aim is to estimate the rate parameters. I created a model containing all differential equations. I used sbiofit which works, but now I got the following two problems:

1. How can I perform parameter estimation/fit on "all my measurements" at once. I do not want to fit measurement Nr. 1, then Nr. 2 etc. Reason: sbiofit will return totally different Parameter Estimates for each measurement. (Yes, this might be because my problem is "over-parametrized" if I fit each measurment at once.)
2. During the fit, sbiofit yiels negative rates. How can I avoid that?

Any hints, also beyond Matlab, are very appreciated!

Best, Vincent

Answer 1

[Arthur Goldsipe]

Hi Vincent,

You can fit multiple experiments with the same set of parameters using sbiofit, but you will need to make several small changes.

1) To fit multiple experiments with the same set of parameters, you need to call sbiofit with the name-value pair 'Pooled',true . In other words, add those to the end of the line of code that calls sbiofit.

2) You need to combine each experiment's data into one groupedData object, marking each experiment with a different ID. This ID variable needs to be set as the grouping variable, so that sbiofit can tell with time points come from the same experiment.

3) You need to set the initial value of A to 0 in the model. You will specify the initial value for each experiment separately.

4) You need to construct a vector of dose objects, one per experiment. These dose objects need to be configured to update the amount of A at time=0 to the appropriate value. The first dose object should set the initial amount of A to the value required for the first group of data, and so forth.

The example file phenobarbdemo shows how to create sample dose objects from a file containing all the data. Although that example uses sbiofitmixed, its syntax is very similar to sbiofit.

To avoid negative rate values, you probably need to constrain your kinetic parameters to be positive. When you describe the parameters you want to estimate (using estimatedInfo objects), specify a log transform. This will cause sbiofit to optimize the log of the parameters, which will force the actual SimBiology parameters to remain positive during estimation.

-Arthur

[Vincent]

Hi Arthur,

Great Answer, everything you suggested worked very fine, thank you! But when I'm starting sbiofit I get the Error

   Invalid dose target 'A' in dose 'Target1'. This object cannot be the Rule variable of any rule.

"A" is the species in my biomodel and I set the InitialValue to 0. Do you have a clue what I could have done wrongly? The Error occured already during compilation of the model (also given within the Error Message)

I gave it some tries and even when using PK-Models I do not get rid of this error…

Thank you so far! -- Vincent

[Arthur Goldsipe]

Hi Vincent,

The error message makes it sound like species A is also being set by a rule in your model. If it's being set by a rule, it's value can't be further changed by a dose.

If your model isn't using any rules, then there may be a bug, or at least a badly written error message. If that's the case, I suggest you contact Technical Support (if you have an active license) and send them your model and reproduction steps. Alternatively, I would suggest creating a new MATLAB Answers question and provide reproduction steps, including attaching any SimBiology project files or MAT files. (I can't easily see when you ask another question as a reply.)

[Vincent]

Hi Arthur,

in this case I don't think that this is a bug. The species A is set by a rule. (Just like every species in rate kinetics.) Without setting the rule, everything works fine, thank you very much!

But one question remains: isn't it some kind of "information loss", if I can't set the rule for my searched parameter?

Thanks, Vincent

[Arthur Goldsipe]

Hi Vincent,

I understand the problem now! I think the intention of this error check is to ensure that you are not using a dose and a repeated assignment rule on the same species, because then you are trying to set the value two different ways. I think it's a mistake that this limitation is also being applied to rate rules. I'm one of the SimBiology developers, so I'll see if we can remove this limitation in a future release.

The reason no one caught this earlier this that it's more common to use reactions than rate rules in SimBiology. And SimBiology will allow you to set the initial value of A using a dose, even when A participates in a reaction.

It's much more convenient to use reactions to set the rate of each species. For example, in your model you must have written 7 rate rules. Instead, you can just write 3 reactions and let SimBiology do the rest of the math. If your model is in variable m1, here's how you might create the first reaction (double-check my math on the reaction rate!):

reaction1 = m1.addreaction('E + F + A + B <-> E + F + Q + B');
reaction1.ReactionRate = 'k12*E*F*A*B - k21*E*F*Q*B';

SimBiology has two ways to set reaction rates. You can either just write the rate equation like I did above, or you can specify the kind of kinetic law (for example, mass action or Michaelis-Menten kinetics) and the parameters of the kinetic law. Here are a couple of other examples of building reactions from the documentation:

• http://www.mathworks.com/help/simbio/ug/representing-a-model-and-model-modifiers-in-simbiology.html
• http://www.mathworks.com/help/simbio/gs/simbiology-command-line-tutorial.html
• http://www.mathworks.com/help/simbio/examples/parameter-scanning-parameter-estimation-and-sensitivity-analysis-in-the-yeast-heterotrimeric-g-protein-cycle.html
• http://www.mathworks.com/help/simbio/examples/deterministic-simulation-of-a-model-containing-a-discontinuity.html

[Vincent]

Hi Arthur,

thank you very much on this answer. Now thinks start to clarify in my mind somehow too :-) Still one question (let's hope it does not go on like that forever ;-) ): I add the ReactionRate-Strings to my reaction-object similar to your description above. Only some reactants are active in each step

   k12*E*A - k21*E*Q

Now I want to set default values for k12 and k21. Therefore I have to assign a Kinetic Law Object to my Reaction-Object. I do that selecting MassAction in this case.

I see that the KineticLaw-Object now contains all 4 species in its property ParameterVariables. And my Reaction-object has forgotten now about my ReactionRate String.

How can I set default values to my model, if I do not use either MassAction, nor Michaelis-Menten-Kinetics but just my own defined ReactionRate-String?

Thank you very much for your help! I appreciate this great customer service :-)

[Arthur Goldsipe]

Hi Vincent,

Since your kinetics are not really mass action, I recommend that you set the kinetic law to Unknown. Then, you can set the reaction rate string to any expression you like.

In addition, it may be helpful to know that there are two types of parameters in SimBiology. I call them "model-scoped" and "reaction-scoped". If a parameter is model-scoped (created by calling addparameter on a model), then this parameter can be used in any reaction. If a parameter is reaction-scoped (created by calling addparameter on a kinetic law), then the parameter can only be used in the one reaction that "owns" it. This can be useful if you always want to name your reaction parameters the same way, like kf or kr, but you need different values of kf for different reactions.

I recommend using model-scoped parameters whenever possible. One reason is because they are easier to use with sbiofit. If you want to estimate a model-scoped parameter k12, then you can create an estimatedInfo object with the name 'k12'. But if you want to estimate a reaction-scoped parameter k12, then you need to specify which reaction the k12 belongs to and set the estimatedInfo's name to something like 'Reaction1.k12'.

[Vincent]

Dear Arthur,

sorry for the long delay. I thank you very much for all of your answers, I figured out a way to do handle my model by now! :-)

Best, Vincent