SimBio: Running Stochastic SSA Simulation of SBML imported model

1 view (last 30 days)
Hello,
I have a biochemical reaction network model defined in the SBML format (attached). I can load the model using the sbmlimport:
sbml_mod = sbmlimport('multistate.xml');
I can perform ODE simulations of this model using sbiosimulate:
[time,x,names] = sbiosimulate(sbml_mod);
However, if I want to perform stochastic SSA (GIllespie) simulations (https://uk.mathworks.com/help/simbio/ug/stochastic-solvers.html), I get problems.
Running
cs = getconfigset(sbml_mod,'active');
cs.SolverType = 'ssa';
cs.StopTime = 30;
[t_ssa, x_ssa] = sbiosimulate(sbml_mod);
I get errors:
--> Error reported from Stochastic Compilation:
Reaction named 'Reaction_2' has an empty kinetic law. For stochastic simulation all kinetic laws must be MassAction.
(one for each of my 18 reactions)
It turns out that when loading SBML models each reactions get no KineticLaw, which is needed for SSA simulations. E.g. sbml_mod.Reactions(1).KineticLaw simply returns []. However, all reactions in the SBML file are simple MassAction reactions, e.g. sbml_mod.Reactions(1).ReactionRate is 'kon*[R(a,l)]*[L(r)]'. I tried simply looping through all reactions and setting their KineticLaw field to 'MassAction'. But this still yields errors like
--> Error reported from KineticLaw Validation:
Parameter variable name on kinetic law '' is empty. The number of parameters on the kinetic law must match the number in its definition, and all
parameter names must be set.
So I somehow need to ftech the right parameters and att that to the MassAction kinetic law. It seems like this should be fairly straightforward (e.g. many other packages, like Copasi can do SSA simulations directly from SBML files). Is there some automatic way of setting all reactions to MassAction kinetic law, or some other way which would enable me to do an SSA simulation of the model?

Accepted Answer

Arthur Goldsipe
Arthur Goldsipe on 31 Jan 2022
First, a little background on why you're seeing this behavior. SBML does not indicate whether a particular reaction follows mass action kinetics. And SimBiology needs to support more general reaction rates and kinetic laws. So SimBiology simply imports SBML reactions by recreating the specifed math as written. And this is the first time I personally have heard anyone request the product to handle this. I don't think we'd want to do the analysis required for this by default, since it could slow down the import of large models. But maybe we could make it an option for import. I'll log this in our database of enhancement requests.
In the mean time, I would probably update the reactions programmatically. To figure out which parameters are used with which reaction, I would take advantage of findUsages. Here's some prototype code to do that, in case it's helpful. (But please note that I have not thoroughly tested it, so use it at your own risk.)
model = sbmlimport('lotka.xml');
% Find all parameters in the model
parameters = sbioselect(model, 'Type', 'parameter');
% See how each parameter is used.
for paramIndex = 1:numel(parameters)
parameter = parameters(paramIndex);
[~, usageTable] = findUsages(parameter);
% Update any reaction that uses this parameter in its rate.
reactions = usageTable.Component(usageTable.Property == "ReactionRate");
for reactionIndex = 1:numel(reactions)
reaction = reactions(reactionIndex);
oldRate = reaction.ReactionRate;
kineticLaw = reaction.KineticLaw;
if isempty(kineticLaw)
% Add a kinetic law
kineticLaw = addkineticlaw(reaction, 'MassAction');
else
% Update the existing kinetic law to mass action
kineticLaw.KineticLawName = 'MassAction';
end
kineticLaw.ParameterVariableNames = parameter.Name;
newRate = reaction.ReactionRate;
if ~strcmp(oldRate, newRate)
warning("Reaction rate for reaction %s changed from %s to %s.", ...
reaction.Name, oldRate, newRate);
end
end
end
Warning: Reaction rate for reaction Reaction1 changed from c1*y1*x to c1*x*y1.
  1 Comment
Torkel Loman
Torkel Loman on 26 Feb 2022
This worked and I can now simualte the model. Loads of thanks, for the help, it was realyl useful.
For future reference, I made a small modification since some paramters gave empty usageTable's, which caused and error:
function [] = clean_ssa_sbml_model(model)
parameters = sbioselect(model, 'Type', 'parameter');
for paramIndex = 1:numel(parameters)
parameter = parameters(paramIndex);
[~, usageTable] = findUsages(parameter);
% Update any reaction that uses this parameter in its rate.
if isempty(usageTable) % if check added by Torkel.
continue
end
reactions = usageTable.Component(usageTable.Property == "ReactionRate");
for reactionIndex = 1:numel(reactions)
reaction = reactions(reactionIndex);
oldRate = reaction.ReactionRate;
kineticLaw = reaction.KineticLaw;
if isempty(kineticLaw)
% Add a kinetic law
kineticLaw = addkineticlaw(reaction, 'MassAction');
else
% Update the existing kinetic law to mass action
kineticLaw.KineticLawName = 'MassAction';
end
kineticLaw.ParameterVariableNames = parameter.Name;
newRate = reaction.ReactionRate;
if ~strcmp(oldRate, newRate)
warning("Reaction rate for reaction %s changed from %s to %s.", ...
reaction.Name, oldRate, newRate);
end
end
end
end

Sign in to comment.

More Answers (0)

Communities

More Answers in the  SimBiology Community

Categories

Find more on Scan Parameter Ranges in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!