During SSA simulation of sbml model using sbiosimulate, observables are not saved.

I have a .sbml model stored in an .xml file. It contains a couple of observables. When I plto the simulation results, these are uniformly set to "1" when I use the ssa simulation, however, observables seems to work fine with ode simulations.

This is the code:

% Fetch model.

model = sbmlimport('multistate.xml');

% Run simulation.

cs = getconfigset(model);

cs.SolverType = 'sundials';

cs.StopTime = 10.0;

[t_ode, x_ode, names] = sbiosimulate(model);

% Plot observables.

plot(t_ode, x_ode(:,end-3:end))

xlabel('Time')

ylabel('States')

legend(names(end-3:end))

% Prepare model for SSA simulations (see https://uk.mathworks.com/matlabcentral/answers/1582584-simbio-running-stochastic-ssa-simulation-of-sbml-imported-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

% Run simulation.

cs = getconfigset(model);

cs.SolverType = 'ssa';

cs.StopTime = 10.0;

[t_ssa, x_ssa, names] = sbiosimulate(model);

% Plot observables.

plot(t_ssa, x_ssa(:,end-3:end))

xlabel('Time')

ylabel('States')

legend(names(end-3:end))

% Non-observables plot still work.

plot(t_ssa, x_ssa(:,1:end-4))

xlabel('Time')

ylabel('States')

legend(names(1:end-4))

% Observable plotting still work for the modified model when using ode.

cs = getconfigset(model);

cs.SolverType = 'sundials';

cs.StopTime = 10.0;

[t_ode_2, x_ode_2, names] = sbiosimulate(model);

plot(t_ode_2, x_ode_2(:,end-3:end))

xlabel('Time')

ylabel('States')

legend(names(end-3:end))

Which produces these plots:

The observables can be plotted when the mode is simulated via ode:

The observables are all "1" when simualted via ssa:

The other species can be plotted without problem though:

Finnaly, plotting the observables when the system has been simualted via ode still work, even after the model have been modified.

There is something wrong where the observables values are not set when the ssa algorithm is used.

Finally, I have attached the model file (extension changed to .txt so that mathworks website would accept it).

Answer 1

[Florian Augustin]

Hi Torkel,

The model in the SBML file defines assignment rules, which are imported into SimBiology as repeated assignment rules. During model simulation with a stochastic solver, SimBiology ignores any rate, assignment, or algebraic rules if present in the model. SimBiology issues a warning on the MATLAB command line about this behavior. So what you are seeing in the plot are the unmodified initial values.

To fix this, you could use SimBiology observables. Repeated assignment rules are continuously evaluated during model simulation and are therefore a useful tool if the states they modify (in your case A_P, A_unbound_P, A_bound_P, RLA_P) are used during simulation. This doesn't seem to be the case in your model. In this case, you can log them using SimBiology observables, which are only evaluated after the simulation is completed. The simulation result is exactly the same in your model, but you will be able to log those states this way. Here is an example how you could do this:

% Fetch model.

model = sbmlimport('multistate.xml');

% Remove observable species and rules from model

observableSpecies = sbioselect(model, "Name", ["A_P", "A_unbound_P", "A_bound_P", "RLA_P"]);

repAsgRules = sbioselect(model, "RuleType", "repeatedAssignment");

delete(observableSpecies);

delete(repAsgRules);

% Add observables (those are the same expressions as used in the assignment rules in the SBML document)

addobservable(model, "A_P", "[A(Y~P,r!1).L(r!2).R(a!1,l!2)]+[A(Y~P,r!1).R(a!1,l)]+[A(Y~P,r)]")'

addobservable(model, "A_unbound_P", "0+[A(Y~P,r)]");

addobservable(model, "A_bound_P", "[A(Y~P,r!1).L(r!2).R(a!1,l!2)]+[A(Y~P,r!1).R(a!1,l)]");

addobservable(model, "RLA_P", "0+[A(Y~P,r!1).L(r!2).R(a!1,l!2)]");

% ... continue with your script.

I hope this helps.

Best,

Florian

[Torkel Loman]

Thank you, this worked. I have however a quick and a long follow-up question:

Quick: in the line addobservable(model, "A_P", "[A(Y~P,r!1).L(r!2).R(a!1,l!2)]+[A(Y~P,r!1).R(a!1,l)]+[A(Y~P,r)]")' you end with ' and not ;. It doesn't seem to matter. Is this a typo, or is there some meaning behin it?

Long: In this case you simply manually write out the rules. However, I have a model where there are several observables, some depending on a large number of species. Is there some way of automating the observables, or would I have to write it manually somehow? (example large model attached in file).

Also, an extra question that probably shoule be an own post, but which might be related: When I simulate the model I get a warning:

Warning: Reported from Stochastic Compilation:

Currently, rules are not supported by stochastic solvers. They will be ignored.

> In sbiosimulate (line 140)

In ssa_sim_model (line 6)

Is this something I should be worried about?

[Florian Augustin]

Hi Torkel,

The ' at the end of the addobservables call is a typo, this should be a semicolon.

To your second question, you can automate the process from my previous answer by looping over all rules that are not required to be evaluated during model simulation and replace them with observables. Here is an example:

% Fetch model.
model = sbmlimport('multistate.xml');

% Remove species and rules from model and replace them with observables.
repAsgRules = sbioselect(model, "RuleType", "repeatedAssignment");
for rule = repAsgRules'
    % Get left- and right-hand sides of rules. The next line assumes that
    % the rule string only contains one equal sign. If this is 
    % not the case, you may have to resort to regexp to split the
    % strings into tokens.
    tokens = strsplit(rule.Rule, " = ");
    [ruleLhs, ruleRhs] = deal(tokens{:});
    % Delete rule and its left-hand side; 
    % this assumes the names of the component on the left-hand 
    % sides of the rules are unique within the model. If this is
    % not the case, you can add more qualifiers to sbioselect.
    component  = sbioselect(model, "Name", ruleLhs);
    ruleString = sbioselect(model, "Name", ruleRhs);
    delete(component);
    delete(rule);
    % Add the observable
    addobservable(model, ruleLhs, ruleRhs);
end

The warning that you are referencing is to alert you that SimBiology rules are ignored during stochastic simulation (see this reference page). You should not use stochastic solvers for models with active repeated assignment rules. In some cases you can replace those rules with observables as descibed above. Then you won't see these warnings when simulating the model with a stochastic solver.

I hope this helps.

-Florian