Adding parameters while programmatically setting up a SimBiology model

I am setting up a SimBiology model programmatically, but am running into errors adding parameters and setting up reactions. The documentation suggests that parameters can be scoped to specific kinetic laws or to the model (in which case they are available to all kinetic laws and reactions). Since some of my parameters are used in several reactions, I figured it would be good to scope them to the model. However, I get the following error among others:

Error using SBCompiler.SimulationObject/simulate
--> Error reported from KineticLaw Validation:
Cannot find parameter with name ''.

To me this reads as if there is a scoping problem - the kinetic law can't find my parameter. I made the stripped down version of my model below:

%%Model
test=sbiomodel('test');
%%Compartments
Acell=addcompartment(test,'Acell',30,'CapacityUnits','picoliter','ConstantCapacity',false);
cytoplasm=addcompartment(Acell,'cytoplasm',15,'CapacityUnits','picoliter','ConstantCapacity',false);
nucleus=addcompartment(Acell,'nucleus',15,'CapacityUnits','picoliter','ConstantCapacity',false);
%%Species
B=addspecies(cytoplasm,'B',0,'InitialAmountUnits','molecule');
A=addspecies(nucleus,'A',1943,'InitialAmountUnits','molecule');
gB=addspecies(nucleus,'gB',4,'InitialAmountUnits','molecule');
gB_A=addspecies(nucleus,'gB_A',0,'InitialAmountUnits','molecule');
%%Parameters
bB_A=addparameter(test,'bB_A',.1,'ValueUnits','1/(molecule*second)');
ubB_A=addparameter(test,'ubB_A',1.0,'ValueUnits','1/(second)');
deg_A=addparameter(test,'deg_A',0.05,'ValueUnits','1/(second)');
deg_B=addparameter(test,'deg_B',0.005,'ValueUnits','1/(second)');
synth_B=addparameter(test,'synth_B',2,'ValueUnits','1/(second)');
%%Reactions
PRODB_A=addreaction(test,'nucleus.gB_A -> cytoplasm.B + nucleus.gB_A','ReactionRate','synth_B*nucleus.gB_A');
BINDgB_A=addreaction(test,'nucleus.gB + nucleus.A <-> nucleus.gB_A','ReactionRate','bB_A*nucleus.gB*nucleus.A - ubB_A*nucleus.gB_A');
DEGA=addreaction(test,'nucleus.A -> null','ReactionRate','deg_A*nucleus.A');
DEGB=addreaction(test,'cytoplasm.B -> null','ReactionRate','deg_B*cytoplasm.B');
%%Kinetic Laws
kPRODB_A=addkineticlaw(PRODB_A,'MassAction');
kBINDgB_A=addkineticlaw(BINDgB_A,'MassAction');
kDEGA=addkineticlaw(DEGA,'MassAction');
kDEGB=addkineticlaw(DEGB,'MassAction');
%%Configuration
cs=getconfigset(test,'active');
cs.SolverType = 'ssa';
cs.StopTime = 50*60;
solver = cs.SolverOptions;
solver.LogDecimation = 10;
cs.CompileOptions.DimensionalAnalysis = false;

When I run this I get several errors:

>> [time,x,names]=sbiosimulate(test);
Error using SBCompiler.SimulationObject/simulate
--> Error reported from KineticLaw Validation:
Cannot find parameter with name ''.

--> Error reported from KineticLaw Validation:
Cannot find parameter with name ''.

--> Error reported from KineticLaw Validation:
Cannot find parameter with name ''.

--> Error reported from KineticLaw Validation:
Cannot find parameter with name ''.

--> Error reported from Expression Validation:
Invalid reaction rate '' for reaction 'nucleus.gB_A -> cytoplasm.B + nucleus.gB_A'. Reaction rates must be valid MATLAB expressions and
cannot end in semicolons, commas, comments ('%' and optional text), or line continuations ('...' and optional text).

--> Error reported from Expression Validation:
Invalid reaction rate '' for reaction 'nucleus.gB + nucleus.A <-> nucleus.geB_A'. Reaction rates must be valid MATLAB expressions and
cannot end in semicolons, commas, comments ('%' and optional text), or line continuations ('...' and optional text).

--> Error reported from Expression Validation:
Invalid reaction rate '' for reaction 'nucleus.A -> null'. Reaction rates must be valid MATLAB expressions and cannot end in semicolons,
commas, comments ('%' and optional text), or line continuations ('...' and optional text).

--> Error reported from Expression Validation:
Invalid reaction rate '' for reaction 'cytoplasm.B -> null'. Reaction rates must be valid MATLAB expressions and cannot end in semicolons,
commas, comments ('%' and optional text), or line continuations ('...' and optional text).

Error in sbiosimulate (line 143)
[t, x] = simobj.simulate(mobj, cs, variants, doses);

I'm quite puzzled as to why the model is incorrectly specified.

One solution is to reorganize everything so that I add a reaction, add a kinetic law, add parameter scoped to that kinetic law, then set the kinetic law constants. But if I do this, how do share a parameter between reactions (not in the stripped down example above, but in my larger model)?

Answer 1

[Scott]

I solved my problem and leave the question up in case it is useful to anyone else.

Each kinetic law just needs to have the parameters it uses specified:

So after the kinetic law statement block above:

%%Add kinetic law constants for each kinetic law
kPRODB_A.ParameterVariableNames={'synth_B'};
kBINDgB_A.ParameterVariableNames={'bB_A','ubB_A'};
kDEGA.ParameterVariableNames={'deg_A'};
kDEGB.ParameterVariableNames={'deg_B'};

Answer 2

[Arthur Goldsipe]

Hi Scott,

You've already found one solution: Set ParameterVariableNames. However, I wanted to provide another solution and give a little more background.

Another solution would be to remove the lines of code that call addkineticlaw.

Basically, SimBiology supports two different ways of determining reaction rates. One is to set the ReactionRate explicitly, which is what you did on the lines of code calling addreaction. The other way is to use a kinetic law. As soon as you set the kinetic law to MassAction by calling addkineticlaw, the value of ReactionRate is replaced with an automatically generated expression based on the value of ParameterVariableNames and any reactants or products.

You can use either approach. But just be aware that setting the reaction rate will change the kinetic law, and setting the kinetic law may change the reaction rate.

-Arthur

[Scott]

Arthur,

Thanks - I have a question about this approach, though. If I take the original code above and just comment out the kinetic laws part, I get the following errors:

>> [time,x,names]=sbiosimulate(test);
Error using SBCompiler.SimulationObject/simulate
--> Error reported from Stochastic Compilation:
Unnamed reaction with reaction rate 'synth_B*nucleus.gB_A' has an empty kinetic law. For stochastic simulation all kinetic laws must be
MassAction.

--> Error reported from Stochastic Compilation:
Unnamed reaction with reaction rate 'bB_A*nucleus.gB*nucleus.A - ubB_A*nucleus.gB_A' has an empty kinetic law. For stochastic simulation
all kinetic laws must be MassAction.

--> Error reported from Stochastic Compilation:
Unnamed reaction with reaction rate 'deg_A*nucleus.A' has an empty kinetic law. For stochastic simulation all kinetic laws must be
MassAction.

--> Error reported from Stochastic Compilation:
Unnamed reaction with reaction rate 'deg_B*cytoplasm.B' has an empty kinetic law. For stochastic simulation all kinetic laws must be
MassAction.

Error in sbiosimulate (line 143)
[t, x] = simobj.simulate(mobj, cs, variants, doses);

I don't get this error if I use an ODE solver (e.g. ode15s). Does SSA require an explicit declaration of kinetic law because it requires MassAction?

Thanks, Scott

[Arthur Goldsipe]

Oh, I hadn't noticed that you were using a stochastic solver. Yes, as you guessed that solver requires mass action kinetic laws. Only deterministic solvers are supported when you set the reaction rate directly.