Gillespie integrator ignores substance unit of model

[matthiaskoenig]

Hi all, I just tried to create some simple stochastic model and reproduce the results in COPASI (see below).

It seems like the roadrunner gillespie integrator is not calculating the particle number if the substance unit of the model is not item, but just uses the substance as particle number! The particle number has to be calculated with the avogadro constant if the substance unit is not item, for instance if it is mole (which it is in many models).

Both simulations below (one with item as substance unit, the other as mole) give exactly the same output, but should be completely different. COPASI complaints correctly about the too large particle number in the system in case of mole.

%matplotlib inline
from __future__ import print_function
import tellurium as te
import numpy as np

r = te.loada("""
model stochastic()
    unit substance = 1 item;

    J1_f : S1 => S2; kf*S1; 
    J1_b : S2 => S1; kb*S2; 
    kf = 0.2; kb = 0.1;
    S1 = 100; S2 = 0;
end
""")

# run some stochastic simulations with the model
r.setIntegrator('gillespie')
r.setSeed(1234)
results = []
for k in range(1, 50):
    r.reset()
    s = r.simulate(0, 40)
    results.append(s)
    r.plot(s, show=False, loc=None, color='black', alpha=0.7)

# run the determenistic version
r.setIntegrator('cvode')
r.reset()
s = r.simulate(0, 40)
r.plot(s, show=True, alpha=0.7)

# save SBML
r.exportToSBML("StochasticExample.xml", current=False)

# If the substance unit is not item, than the particle number is false
r = te.loada("""
model stochastic()
    unit substance = 1 mole;

    J1_f : S1 => S2; kf*S1; 
    J1_b : S2 => S1; kb*S2; 
    kf = 0.2; kb = 0.1;
    S1 = 100; S2 = 0;
end
""")

# run some stochastic simulations with the model
r.setIntegrator('gillespie')
r.setSeed(1234)
results = []
for k in range(1, 50):
    r.reset()
    s = r.simulate(0, 40)
    results.append(s)
    r.plot(s, show=False, loc=None, color='black', alpha=0.7)

# run the determenistic version
r.setIntegrator('cvode')
r.reset()
s = r.simulate(0, 40)
r.plot(s, show=True, alpha=0.7)

<?xml version="1.0" encoding="UTF-8"?>
<!-- Created by libAntimony version v2.9.0 with libSBML version 5.12.1. -->
<sbml xmlns="http://www.sbml.org/sbml/level3/version1/core" level="3" version="1">
  <model id="stochastic" name="stochastic">
    <listOfCompartments>
      <compartment sboTerm="SBO:0000410" id="default_compartment" spatialDimensions="3" size="1" constant="true"/>
    </listOfCompartments>
    <listOfSpecies>
      <species id="S1" compartment="default_compartment" initialConcentration="100" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
      <species id="S2" compartment="default_compartment" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
    </listOfSpecies>
    <listOfParameters>
      <parameter id="kf" value="0.2" constant="true"/>
      <parameter id="kb" value="0.1" constant="true"/>
    </listOfParameters>
    <listOfReactions>
      <reaction id="J1_f" reversible="false" fast="false">
        <listOfReactants>
          <speciesReference species="S1" stoichiometry="1" constant="true"/>
        </listOfReactants>
        <listOfProducts>
          <speciesReference species="S2" stoichiometry="1" constant="true"/>
        </listOfProducts>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML">
            <apply>
              <times/>
              <ci> kf </ci>
              <ci> S1 </ci>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
      <reaction id="J1_b" reversible="false" fast="false">
        <listOfReactants>
          <speciesReference species="S2" stoichiometry="1" constant="true"/>
        </listOfReactants>
        <listOfProducts>
          <speciesReference species="S1" stoichiometry="1" constant="true"/>
        </listOfProducts>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML">
            <apply>
              <times/>
              <ci> kb </ci>
              <ci> S2 </ci>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
    </listOfReactions>
  </model>
</sbml>

[hsauro]

This s very likely true. Should be a simple fix though.

[0u812]

I think part of the reason no one has implemented this yet is that the SBML documentation is clear as mud when it comes to species units. In the XML posted above, the word "mole" doesn't even show up, so I presume it is the default unit for amount but even after reading this documentation I still don't see where that default choice is explained.

[luciansmith]

What exactly is unclear about species units in the documentation?

It was decided very early on in SBML's development that units were always arbitrary, and that no conversion would ever be necessary in any model for any reason. If people are converting moles to items using avogadro, and avogadro's number doesn't show up in the model itself, that's a use of SBML that was not intended. If it's a common convention, it's something that libroadrunner may or may not want to follow, but no simulator should ever have to look at the units of a model for any reason. Units are there for cross-checking and validation only, and should not change any math.

[0u812]

Like I said in my original message, the unit 'moles' does not appear in the XML Matthias posted, but it does appear in the Antimony. So where is roadrunner supposed to get the information that describes the species units?

[luciansmith]

Hmm, I'm not sure where the XML that Matthias posted was obtained from. When I translate

model stochastic()
    unit substance = 1 mole;

   J1_f : S1 => S2; kf*S1;
   J1_b : S2 => S1; kb*S2;
   kf = 0.2; kb = 0.1;
   S1 = 100; S2 = 0;
end

to SBML, I get:

<?xml version="1.0" encoding="UTF-8"?>

<!-- Created by libAntimony version v2.9.0 with libSBML version 5.12.1. -->

<sbml xmlns="http://www.sbml.org/sbml/level3/version1/core" xmlns:comp="
http://www.sbml.org/sbml/level3/version1/comp/version1" level="3"
version="1" comp:required="true">
  <model id="stochastic" name="stochastic" substanceUnits="substance">
    <listOfUnitDefinitions>
      <unitDefinition id="substance">
        <listOfUnits>
          <unit kind="mole" exponent="1" scale="0" multiplier="1"/>
        </listOfUnits>
      </unitDefinition>
    </listOfUnitDefinitions>
    <listOfCompartments>
      <compartment sboTerm="SBO:0000410" id="default_compartment"
spatialDimensions="3" size="1" constant="true"/>
    </listOfCompartments>
    <listOfSpecies>
      <species id="S1" compartment="default_compartment"
initialConcentration="100" hasOnlySubstanceUnits="false"
boundaryCondition="false" constant="false"/>
      <species id="S2" compartment="default_compartment"
initialConcentration="0" hasOnlySubstanceUnits="false"
boundaryCondition="false" constant="false"/>
    </listOfSpecies>
    <listOfParameters>
      <parameter id="kf" value="0.2" constant="true"/>
      <parameter id="kb" value="0.1" constant="true"/>
    </listOfParameters>
    <listOfReactions>
      <reaction id="J1_f" reversible="false" fast="false">
        <listOfReactants>
          <speciesReference species="S1" stoichiometry="1" constant="true"/>
        </listOfReactants>
        <listOfProducts>
          <speciesReference species="S2" stoichiometry="1" constant="true"/>
        </listOfProducts>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML">
            <apply>
              <times/>
              <ci> kf </ci>
              <ci> S1 </ci>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
      <reaction id="J1_b" reversible="false" fast="false">
        <listOfReactants>
          <speciesReference species="S2" stoichiometry="1" constant="true"/>
        </listOfReactants>
        <listOfProducts>
          <speciesReference species="S1" stoichiometry="1" constant="true"/>
        </listOfProducts>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML">
            <apply>
              <times/>
              <ci> kb </ci>
              <ci> S2 </ci>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
    </listOfReactions>
  </model>
</sbml>

which is awkward, but at least does mention moles.

[0u812]

Okay, that makes a lot more sense. If that's the XML I should be able to fix this pretty easily. Thanks!

[0u812]

Btw everyone here should be aware that this issue has nothing to do with stochastic models. If I make this change it will affect EVERY model that uses moles as a substance unit, regardless of the timecourse solver (so all amounts will get multiplied by Avogadro's constant). I suspect roadrunner also ignores the hasOnlySubstanceUnits attribute, which is a separate bug and would be much harder to fix.

[hsauro]

Can we just consider the units for stochastic simulations because this is where the difference was found between copasi and roadrunner?

[0u812]

That's not possible because you choose the solver separately from the model. The model gets compiled once and used across solvers.

[luciansmith]

'hasOnlySubstanceUnits' doesn't have anything to do with moles or avogadro or stochastics or anything. It means 'when the symbol for this species appears in a formula, use its amount, not its concentration'.

Again, units are designed to only be annotations, and should never change the behavior of a simulator. If Copasi is multiplying things by avogadro or not based on the units, that is a bug in Copasi, as far as I can tell.

[0u812]

Yes but roadrunner will internally multiply the initial amount by Avogadro's constant when substance units is set to moles if I implement this. Sorry for the confusion re hasOnlySubstanceUnits.

[hsauro]

Maybe Matthias can give more info on what's going on. He'll be up in 8 hours or so. I've not looked at the model he sent as I've been doing other work.

[matthiaskoenig]

Sorry for the confusion, I did not post the latest SBML. Of course it contains the substance=mole or item respectively. Here updated code and models (SBML and copasi files and zip)

%matplotlib inline
from __future__ import print_function
import tellurium as te
import numpy as np

r_item = te.loada("""
model stochastic_item()
    unit substance = 1 item;

    J1_f : S1 => S2; kf*S1; 
    J1_b : S2 => S1; kb*S2; 
    kf = 0.2; kb = 0.1;
    S1 = 100; S2 = 0;
end
""")

r_mole = te.loada("""
model stochastic_mole()
    unit substance = 1 mole;

    J1_f : S1 => S2; kf*S1; 
    J1_b : S2 => S1; kb*S2; 
    kf = 0.2; kb = 0.1;
    S1 = 100; S2 = 0;
end
""")

for name in ['r_item', 'r_mole']:
    r = globals()[name]
    # run some stochastic simulations with the model
    r.setIntegrator('gillespie')
    r.setSeed(1234)
    results = []
    for k in range(1, 50):
        r.reset()
        s = r.simulate(0, 40)
        results.append(s)
        r.plot(s, show=False, loc=None, color='black', alpha=0.7)

    # run the determenistic version
    r.setIntegrator('cvode')
    r.reset()
    s = r.simulate(0, 40)
    r.plot(s, show=True, alpha=0.7, title=name)

    # save SBML
    r.exportToSBML("{}.xml".format(name), current=False)

item_mole_sbml.zip

Yes, Copasi gives different results depending on the substance units on the model element. They calculate the particle number depending on the the units. Not sure if bug or feature in Copasi. As a consequence all simulations with stochastic models which do not use item as substance unit result in different simulation results between roadrunner and Copasi.

@luciansmith It would be very nice to have the units as simple as possible on the model element, i.e. directly have item or mole instead of going via an additional unitDefinition (see https://sourceforge.net/p/antimony/bug-reports/37/). This would make it also much easier to implement the conversion.

[matthiaskoenig]

It was decided very early on in SBML's development that units were always arbitrary, and that no conversion would ever be necessary in any model for any reason. If people are converting moles to items using avogadro, and avogadro's number doesn't show up in the model itself, that's a use of SBML that was not intended. If it's a common convention, it's something that libroadrunner may or may not want to follow, but no simulator should ever have to look at the units of a model for any reason. Units are there for cross-checking and validation only, and should not change any math.

This is my opinion also. There have to be tests for this in the the SBML stochastic test suite! It has to be clear what is the right simulation result for a given model and in my opinion the current implementation of COPASI is wrong. It is convenient & probably what the user wanted when writing mole as substance unit, but from the perspective of "units are only annotations" it is wrong. It must be possible to remove all unit information from the model and still get the same simulation results.

There should be tests for that in the stochastic SBML test suite & COPASI should fail this tests. It must be clear and agreed on in the community what the way is to simulate models, so that different tools can generate the identical numerical results (in the stochastic simulation the same mean trajectories with same variance on all curves). There should be a test which tests that the two models I attached give the same results.

For roadrunner

• roadrunner should give a warning that one tries to do a stochastic simulation with substance_units != item, because it is probably not what the user intended
• there must be a way to get the conversion to items in roadrunner (or tell roadrunner to do this conversion somehow). Almost all my models are in mole, but I want to be able to run stochastic simulations with them without the requirement to create a second model with all units/numerical values converted to items.

[0u812]

I'm not clear on what fix should be. Based on the above posts by Lucian and Matthias, it sounds like I shouldn't multiply the initial amount by Avogadro's number. Yet, when you run a stochastic simulation the state vector is usually expressed in number of molecules (i.e. items if you prefer). Can someone clarify exactly what the fix should be?

[matthiaskoenig]

Here is what I would find a good solution:

• default behavior is like now (this is in agreement with SBML, i.e. units should not matter)
• there should be a warning if the substance/extent unit is not item but mole ("Warning: The model substance/extent units are not items in stochastic simulation, to convert the model to items use the setting ...)
• There should be an option to convert a model on loading from mole to item, which allows to do what COPASI is doing right now (I am not sure what has to be adapted, and if this is easy to do. Probably everything has to be scaled in the model, i.e. species, but also reactions) . This could be one of the SBML loading options.

[luciansmith]

I think Matthias's proposal could work, too. It's probably worth a a discussion on sbml-discuss that points out the discrepancy about using units to affect the mathematics in this way--maybe everyone who does stochastic modeling just accepts this, or something. If there is consensus, adding tests to the stochastic test suite is definitely a good idea.

As far as specific implementation, I think I would implement a special simulator interface and/or passed-in option to the stochastic simulator itself that indicated that the model was to be converted from moles to items before simulation. If that's possible.