search

sys-bio / roadrunner

libRoadRunner: A high-performance SBML simulator
http://libroadrunner.org/
Other
39 stars 24 forks source link

Negative values running Gillespie stochastic simulations. #342

Closed J-Revell closed 8 years ago

J-Revell commented 8 years ago

In my models, running Gillespie stochastic simulations has been producing negative values (!) for species.

I have traced the issue back to reactions of the form X + X -> Y and produced a simple SBML for this one reaction, 2 species case. For example, starting with a species value X(t=0) = 9, will produce a trajectory of states of X = 9, 7, 5, 3, 1, -1, -3... and so forth, while starting with X(t=0) = 8 reaches a minimum at 0 as expected.

Is there a problem with the way RoadRunner is calculating the stochastic propensities for the reactants (that should be of binomial coefficient form --- [X]*([X]-1)/2 in this case)?

Models were produced in COPASI. Could it also be an issue with how COPASI exports model information to SBML? I note that the produced SBML gives the rate law [X]^2. But even so, the stoichiometry of the reactants should be sufficient information to be able to produce correct stochastic simulations.

Many thanks,

<?xml version="1.0" encoding="UTF-8"?>
<!-- Created by COPASI version 4.16 (Build 104) on 2016-08-03 14:02 with libSBML version 5.11.1. -->
<sbml xmlns="http://www.sbml.org/sbml/level2/version4" level="2" version="4">
  <model metaid="COPASI0" id="ExampleCase" name="ExampleCase">
    <annotation>
      <COPASI xmlns="http://www.copasi.org/static/sbml">
        <rdf:RDF xmlns:dcterms="http://purl.org/dc/terms/" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
          <rdf:Description rdf:about="#COPASI0">
            <dcterms:created>
              <rdf:Description>
                <dcterms:W3CDTF>2016-08-03T13:57:27Z</dcterms:W3CDTF>
              </rdf:Description>
            </dcterms:created>
          </rdf:Description>
        </rdf:RDF>
      </COPASI>
    </annotation>
    <listOfUnitDefinitions>
      <unitDefinition id="volume" name="volume">
        <listOfUnits>
          <unit kind="dimensionless" exponent="1" scale="0" multiplier="1"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="substance" name="substance">
        <listOfUnits>
          <unit kind="dimensionless" exponent="1" scale="0" multiplier="1"/>
        </listOfUnits>
      </unitDefinition>
    </listOfUnitDefinitions>
    <listOfCompartments>
      <compartment id="compartment" name="compartment" spatialDimensions="3" size="1" constant="true"/>
    </listOfCompartments>
    <listOfSpecies>
      <species metaid="COPASI1" id="Y" name="Y" compartment="compartment" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
        <annotation>
          <COPASI xmlns="http://www.copasi.org/static/sbml">
            <rdf:RDF xmlns:dcterms="http://purl.org/dc/terms/" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
              <rdf:Description rdf:about="#COPASI1">
                <dcterms:created>
                  <rdf:Description>
                    <dcterms:W3CDTF>2016-08-03T13:58:41Z</dcterms:W3CDTF>
                  </rdf:Description>
                </dcterms:created>
              </rdf:Description>
            </rdf:RDF>
          </COPASI>
        </annotation>
      </species>
      <species metaid="COPASI2" id="X" name="X" compartment="compartment" initialConcentration="9" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false">
        <annotation>
          <COPASI xmlns="http://www.copasi.org/static/sbml">
            <rdf:RDF xmlns:dcterms="http://purl.org/dc/terms/" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
              <rdf:Description rdf:about="#COPASI2">
                <dcterms:created>
                  <rdf:Description>
                    <dcterms:W3CDTF>2016-08-03T13:58:22Z</dcterms:W3CDTF>
                  </rdf:Description>
                </dcterms:created>
              </rdf:Description>
            </rdf:RDF>
          </COPASI>
        </annotation>
      </species>
    </listOfSpecies>
    <listOfReactions>
      <reaction metaid="COPASI3" id="Reaction_1" name="Reaction_1" reversible="false">
        <annotation>
          <COPASI xmlns="http://www.copasi.org/static/sbml">
            <rdf:RDF xmlns:dcterms="http://purl.org/dc/terms/" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
              <rdf:Description rdf:about="#COPASI3">
                <dcterms:created>
                  <rdf:Description>
                    <dcterms:W3CDTF>2016-08-03T13:57:48Z</dcterms:W3CDTF>
                  </rdf:Description>
                </dcterms:created>
              </rdf:Description>
            </rdf:RDF>
          </COPASI>
        </annotation>
        <listOfReactants>
          <speciesReference species="X" stoichiometry="2"/>
        </listOfReactants>
        <listOfProducts>
          <speciesReference species="Y" stoichiometry="1"/>
        </listOfProducts>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML">
            <apply>
              <times/>
              <ci> compartment </ci>
              <ci> k1 </ci>
              <apply>
                <power/>
                <ci> X </ci>
                <cn> 2 </cn>
              </apply>
            </apply>
          </math>
          <listOfParameters>
            <parameter id="k1" name="k1" value="0.1"/>
          </listOfParameters>
        </kineticLaw>
      </reaction>
    </listOfReactions>
  </model>
</sbml>
hsauro commented 8 years ago

Will take a look. Gillespie should never produce negative numbers assuming the model is ok.

Herbert Sauro

Herbert

On Wednesday, August 3, 2016, Jeremy Revell notifications@github.com wrote:

In my models, running Gillespie stochastic simulations has been producing negative values (!) for species.

I have traced the issue back to reactions of the form X + X -> Y and produced a simple SBML for this one reaction, 2 species case. For example, starting with a species value X(t=0) = 9, will produce a trajectory of states of X = 9, 7, 5, 3, 1, -1, -3... and so forth, while starting with X(t=0) = 8 reaches a minimum at 0 as expected.

Is there a problem with the way RoadRunner is calculating the stochastic propensities for the reactants (that should be of binomial coefficient form --- [X]*([X]-1)/2 in this case)?

Models were produced in COPASI. Could it also be an issue with how COPASI exports model information to SBML? I note that the produced SBML gives the rate law [X]^2. But even so, the stoichiometry of the reactants should be sufficient information to be able to produce correct stochastic simulations.

Many thanks,

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

dcterms:created rdf:Description dcterms:W3CDTF2016-08-03T13:57:27Z/dcterms:W3CDTF /rdf:Description /dcterms:created /rdf:Description /rdf:RDF dcterms:created rdf:Description dcterms:W3CDTF2016-08-03T13:58:41Z/dcterms:W3CDTF /rdf:Description /dcterms:created /rdf:Description /rdf:RDF dcterms:created rdf:Description dcterms:W3CDTF2016-08-03T13:58:22Z/dcterms:W3CDTF /rdf:Description /dcterms:created /rdf:Description /rdf:RDF dcterms:created rdf:Description dcterms:W3CDTF2016-08-03T13:57:48Z/dcterms:W3CDTF /rdf:Description /dcterms:created /rdf:Description /rdf:RDF compartment k1 X 2

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/342, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDhy8PwEbvqxA1Rlh2UCfqAfR7F0Hks5qcJ9vgaJpZM4Jbq5L .

0u812 commented 8 years ago

Linked #97

hsauro commented 8 years ago

I had a look and can confirm the negative concentration. The reason has to do with the rate law which is currently k_X_X. and the stoichoimetry which is 2X ->

Assume for the moment that X = 1, this means it is possible, according to the rate law (X*X), for the reaction to fire. If it does, the stoichiometry tell us to take 2 molecules away resulting in -1 which is what you're seeing.

The solution to the problem is to use the more proper bimoleular rate law kx(x-1)/2. This is because the probability that this reaction will take place is proportional to the number of ways two molecules of type X can be grouped as an unordered pair, that is n!/[(n-k)k!] where k is 2 and n the number of molecules of x. When X = 1 the popensity function is equal to zero, since it is impossible for one species X to react on its own. You have to specifiy the appropriate rate law, roadrunner won't convert say x^2 to x*(x-1)/2.

J-Revell commented 8 years ago

Many thanks for the swift response. I thought that this might be the case.

Jeremy