SimBio: Running Stochastic SSA Simulation of SBML imported model

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?

Answer 1

[Arthur Goldsipe]

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.

[Torkel Loman]

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