SBML Observables set to 0 when non-default solver is used

[TorkelE]

I have an SBML model which I want to simulate with the SSA using Gillespy2 (in attached .txt file, was .xml but changed it so that GitHub would accept). However, I noticed that the output plots were weird. It turns out that the observable fields of the result were all set to 0 at all time points. However, this only happened when a non-default solver was used:

import sys, os, numpy, libsbml, gillespy2
result = gillespy2.core.import_SBML('multistate.txt')[0].run();
result['A_P']

gives:

array([ 0.21219319,  0.21219319,  0.65106383,  1.14426936,  1.65573076,
        2.17529173,  3.22380921,  3.74917415,  4.27395738,  4.79757012,
        5.31959544,  5.83973699,  6.35778071,  6.35778071,  6.87357039,
        7.38699153,  7.89796027,  8.40641555,  8.91231352,  9.41562344,
        9.91632461, 10.41440406, 10.9098548 , 11.40267448, 11.89286439,
...

while

import sys, os, numpy, libsbml, gillespy2
from gillespy2 import NumPySSASolver
result = gillespy2.core.import_SBML('multistate.txt')[0].run(solver=NumPySSASolver);
result['A_P']

gives:

array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
...

similar things happen when I try solver=ODESolver, solver=TauLeapingSolver. However, specifying solver=TauHybridSolver (the default, right?) does give non-zero outputs.
multistate.txt

[jtcooper10]

Hello there!

The default solver is selected based on what elements are present in the model. In the case of this particular model, the TauHybridSolver is selected by default because it includes assignment rules.

Assignment rules are, unfortunately, only supported by the TauHybridSolver at the moment. The reason it looks like your Observables are "zeroed out" when using other solvers is because the assignment rules that modify them are not applied to the simulation.

(I will admit, though, that it's not particularly clear from the documentation that this is the case, which is worth making a note of and addressing)

Hope this helps!

[TorkelE]

Just double checking here, you mean that

import sys, os, numpy, libsbml, gillespy2
from gillespy2 import NumPySSASolver
result = gillespy2.core.import_SBML('multistate.txt')[0].run(solver=NumPySSASolver);

is bad code, because the model uses "assignment rules", which is unsupported by the solver. When I run it the results look fine, it is just the observables that are weird.

Do you mean that the combination of assignment rules and solver:

• Messes up observables, but everything is otherwise fine.
• Even though it looks that non-observables stuff is fine, the entire solver is messed up because of this incompatibility.
• ?

[jtcooper10]

Yes, that's correct. I opened an issue, #721, to make this behavior clearer, but essentially the behavior of the Observables is being modeled using an SBML feature that only works when using the TauHybridSolver. Using any other solver (including the NumPySSASolver) for this will have unexpected results (in this case, everything works except the Observables).