search

sys-bio / roadrunner

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

Summation theorem MCA #127

Closed 0u812 closed 9 years ago

0u812 commented 9 years ago

Via Matthias: A] I found some problems in the Metabolic Control Analysis of RoadRunner. Somehow the scaled flux control coefficients to not add up to 1 for transporters. I have an example attached. bA, bB and bC are transporters. For the internal reactions v1, v2, v3, v4 summation works (mca.py)

** Summation theorem flux control (=1) **
bA    0.000003
bB    0.000024
bC   -0.000033
v1    1.000000
v2    1.000000
v3    1.000000
v4    1.000000

bA    False
bB    False
bC    False
v1     True
v2     True
v3     True
v4     True

B] Unable to calculate CC I could not calculate the following CC in the example mca.py:

# Scaled control coefficient with respect to a global parameter.
# ! This should give a numerical value, but is nan ???
print('CC:', rr.getCC('bA', 'Vmax_bA'))

In my opinion this should give a numerical value, but I am not sure if I misunderstood something.

0u812 commented 9 years ago

Duplicate of #94?

0u812 commented 9 years ago

mca.py:

"""
Performing metabolic control analysis with a simple model.
"""
from __future__ import print_function
import pandas as pd
import numpy as np
from pandas import DataFrame

import roadrunner
roadrunner.Config.setValue(roadrunner.Config.LOADSBMLOPTIONS_CONSERVED_MOIETIES, True)

# equidistant time course for 20 [s]
rr = roadrunner.RoadRunner("Koenig_demo_v02.xml")
s = rr.simulate(0, 40, 100, plot=True)

rr.getSteadyStateValues()

# Scaled control coefficient with respect to a global parameter.
# ! This should give a numerical value, but is nan ???
print('CC:', rr.getCC('bA', 'Vmax_bA'))

# The n by n matrix of scaled flux control coefficients 
# where n is the number of reactions.
print('** Scaled Flux Control Coefficients **')
C_J = DataFrame(rr.getScaledFluxControlCoefficientMatrix(), \
                index=rr.model.getReactionIds(), \
                columns=rr.model.getReactionIds())
print(C_J)

# The m by n matrix of scaled concentration control coefficients where m is the number
# of floating species and n the number of reactions.
print('** Scaled Concentration Control Coefficients **')
C_S = DataFrame(rr.getScaledConcentrationControlCoefficientMatrix(), \
                index=rr.model.getFloatingSpeciesIds(), \
                columns=rr.model.getReactionIds())
print(C_S)

# Flux control coefficients have to sum to 1
print('** Summation theorem flux control (=1) **')
print(np.sum(C_J, axis=1))
print(abs(np.sum(C_J, axis=1)-1) < 1E-6)
# ! The three transporters do not sum up to 1 ???

# Concentration control coefficients have to sum to 0
print('** Summation theorem concentration control (=0) **')
print(np.sum(C_S, axis=1))
print(np.sum(C_S, axis=1) < 1E-6)
0u812 commented 9 years ago

Koenig_demo_v02.xml:

<?xml version="1.0" encoding="UTF-8"?>
<!-- Created by libAntimony version v2.7.0 on 2015-05-18 10:03 with libSBML version 5.11.4. -->
<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="Koenig_demo_v02" name="Koenig_demo_v02" substanceUnits="substance" timeUnits="time_unit" volumeUnits="volume" areaUnits="area" lengthUnits="length" extentUnits="extent">
    <listOfUnitDefinitions>
      <unitDefinition id="mM">
        <listOfUnits>
          <unit kind="mole" exponent="1" scale="0" multiplier="1"/>
          <unit kind="metre" exponent="-3" scale="0" multiplier="1"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="m3">
        <listOfUnits>
          <unit kind="metre" exponent="3" scale="0" multiplier="1"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="m2">
        <listOfUnits>
          <unit kind="metre" exponent="2" scale="0" multiplier="1"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="mole_per_s">
        <listOfUnits>
          <unit kind="mole" exponent="1" scale="0" multiplier="1"/>
          <unit kind="second" exponent="-1" scale="0" multiplier="1"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="substance">
        <listOfUnits>
          <unit kind="mole" exponent="1" scale="0" multiplier="1"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="volume">
        <listOfUnits>
          <unit kind="metre" exponent="3" scale="0" multiplier="1"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="area">
        <listOfUnits>
          <unit kind="metre" exponent="2" scale="0" multiplier="1"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="length">
        <listOfUnits>
          <unit kind="metre" exponent="1" scale="0" multiplier="1"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="extent">
        <listOfUnits>
          <unit kind="mole" exponent="1" scale="0" multiplier="1"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="time_unit">
        <listOfUnits>
          <unit kind="second" exponent="1" scale="0" multiplier="1"/>
        </listOfUnits>
      </unitDefinition>
    </listOfUnitDefinitions>
    <listOfCompartments>
      <compartment id="extern" name="External Compartment" spatialDimensions="3" size="1e-06" units="m3" constant="true"/>
      <compartment id="cell" name="Cell Compartment" spatialDimensions="3" size="1e-06" units="m3" constant="true"/>
      <compartment id="membrane" name="Cell Membrane" spatialDimensions="3" size="1" units="m2" constant="true"/>
    </listOfCompartments>
    <listOfSpecies>
      <species id="A" name="A cell" compartment="cell" initialConcentration="0" substanceUnits="mole" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
      <species id="B" name="B cell" compartment="cell" initialConcentration="0" substanceUnits="mole" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
      <species id="C" name="C cell" compartment="cell" initialConcentration="0" substanceUnits="mole" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
      <species id="A_ext" name="A extern" compartment="extern" initialConcentration="10" substanceUnits="mole" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
      <species id="B_ext" name="B extern" compartment="extern" initialConcentration="0" substanceUnits="mole" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
      <species id="C_ext" name="C extern" compartment="extern" initialConcentration="0" substanceUnits="mole" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
    </listOfSpecies>
    <listOfParameters>
      <parameter id="scale_f" value="1e-06" units="dimensionless" constant="true"/>
      <parameter id="Vmax_bA" value="5" units="mole_per_s" constant="true"/>
      <parameter id="Km_A" value="1" units="mM" constant="true"/>
      <parameter id="Vmax_bB" value="2" units="mole_per_s" constant="true"/>
      <parameter id="Km_B" value="0.5" units="mM" constant="true"/>
      <parameter id="Vmax_bC" value="2" units="mole_per_s" constant="true"/>
      <parameter id="Km_C" value="3" units="mM" constant="true"/>
      <parameter id="Vmax_v1" value="1" units="mole_per_s" constant="true"/>
      <parameter id="Keq_v1" value="10" units="dimensionless" constant="true"/>
      <parameter id="Vmax_v2" value="0.5" units="mole_per_s" constant="true"/>
      <parameter id="Vmax_v3" value="0.5" units="mole_per_s" constant="true"/>
      <parameter id="Vmax_v4" value="0.5" units="mole_per_s" constant="true"/>
      <parameter id="Keq_v4" value="2" units="dimensionless" constant="true"/>
    </listOfParameters>
    <listOfReactions>
      <reaction id="bA" name="A import" reversible="false" fast="false" compartment="membrane">
        <listOfReactants>
          <speciesReference species="A_ext" stoichiometry="1" constant="true"/>
        </listOfReactants>
        <listOfProducts>
          <speciesReference species="A" stoichiometry="1" constant="true"/>
        </listOfProducts>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML" xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
            <apply>
              <divide/>
              <apply>
                <times/>
                <ci> scale_f </ci>
                <apply>
                  <divide/>
                  <ci> Vmax_bA </ci>
                  <ci> Km_A </ci>
                </apply>
                <apply>
                  <minus/>
                  <ci> A_ext </ci>
                  <ci> A </ci>
                </apply>
              </apply>
              <apply>
                <plus/>
                <cn sbml:units="dimensionless" type="integer"> 1 </cn>
                <apply>
                  <divide/>
                  <ci> A_ext </ci>
                  <ci> Km_A </ci>
                </apply>
                <apply>
                  <divide/>
                  <ci> A </ci>
                  <ci> Km_A </ci>
                </apply>
              </apply>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
      <reaction id="bB" name="B export" reversible="false" fast="false" compartment="membrane">
        <listOfReactants>
          <speciesReference species="B" stoichiometry="1" constant="true"/>
        </listOfReactants>
        <listOfProducts>
          <speciesReference species="B_ext" stoichiometry="1" constant="true"/>
        </listOfProducts>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML" xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
            <apply>
              <divide/>
              <apply>
                <times/>
                <ci> scale_f </ci>
                <apply>
                  <divide/>
                  <ci> Vmax_bB </ci>
                  <ci> Km_B </ci>
                </apply>
                <apply>
                  <minus/>
                  <ci> B </ci>
                  <ci> B_ext </ci>
                </apply>
              </apply>
              <apply>
                <plus/>
                <cn sbml:units="dimensionless" type="integer"> 1 </cn>
                <apply>
                  <divide/>
                  <ci> B_ext </ci>
                  <ci> Km_B </ci>
                </apply>
                <apply>
                  <divide/>
                  <ci> B </ci>
                  <ci> Km_B </ci>
                </apply>
              </apply>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
      <reaction id="bC" name="C export" reversible="false" fast="false" compartment="membrane">
        <listOfReactants>
          <speciesReference species="C" stoichiometry="1" constant="true"/>
        </listOfReactants>
        <listOfProducts>
          <speciesReference species="C_ext" stoichiometry="1" constant="true"/>
        </listOfProducts>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML" xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
            <apply>
              <divide/>
              <apply>
                <times/>
                <ci> scale_f </ci>
                <apply>
                  <divide/>
                  <ci> Vmax_bC </ci>
                  <ci> Km_C </ci>
                </apply>
                <apply>
                  <minus/>
                  <ci> C </ci>
                  <ci> C_ext </ci>
                </apply>
              </apply>
              <apply>
                <plus/>
                <cn sbml:units="dimensionless" type="integer"> 1 </cn>
                <apply>
                  <divide/>
                  <ci> C_ext </ci>
                  <ci> Km_C </ci>
                </apply>
                <apply>
                  <divide/>
                  <ci> C </ci>
                  <ci> Km_C </ci>
                </apply>
              </apply>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
      <reaction id="v1" name="A -&gt; B" reversible="true" fast="false" compartment="cell">
        <listOfReactants>
          <speciesReference species="A" stoichiometry="1" constant="true"/>
        </listOfReactants>
        <listOfProducts>
          <speciesReference species="B" stoichiometry="1" constant="true"/>
        </listOfProducts>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML" xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
            <apply>
              <times/>
              <apply>
                <divide/>
                <apply>
                  <times/>
                  <ci> scale_f </ci>
                  <ci> Vmax_v1 </ci>
                </apply>
                <ci> Km_A </ci>
              </apply>
              <apply>
                <minus/>
                <ci> A </ci>
                <apply>
                  <times/>
                  <apply>
                    <divide/>
                    <cn sbml:units="dimensionless" type="integer"> 1 </cn>
                    <ci> Keq_v1 </ci>
                  </apply>
                  <ci> B </ci>
                </apply>
              </apply>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
      <reaction id="v2" name="A -&gt; C" reversible="true" fast="false" compartment="cell">
        <listOfReactants>
          <speciesReference species="A" stoichiometry="1" constant="true"/>
        </listOfReactants>
        <listOfProducts>
          <speciesReference species="C" stoichiometry="1" constant="true"/>
        </listOfProducts>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML">
            <apply>
              <times/>
              <apply>
                <divide/>
                <apply>
                  <times/>
                  <ci> scale_f </ci>
                  <ci> Vmax_v2 </ci>
                </apply>
                <ci> Km_A </ci>
              </apply>
              <ci> A </ci>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
      <reaction id="v3" name="C -&gt; A" reversible="true" fast="false" compartment="cell">
        <listOfReactants>
          <speciesReference species="C" stoichiometry="1" constant="true"/>
        </listOfReactants>
        <listOfProducts>
          <speciesReference species="A" stoichiometry="1" constant="true"/>
        </listOfProducts>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML">
            <apply>
              <times/>
              <apply>
                <divide/>
                <apply>
                  <times/>
                  <ci> scale_f </ci>
                  <ci> Vmax_v3 </ci>
                </apply>
                <ci> Km_A </ci>
              </apply>
              <ci> C </ci>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
      <reaction id="v4" name="C -&gt; B" reversible="true" fast="false" compartment="cell">
        <listOfReactants>
          <speciesReference species="C" stoichiometry="1" constant="true"/>
        </listOfReactants>
        <listOfProducts>
          <speciesReference species="B" stoichiometry="1" constant="true"/>
        </listOfProducts>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML" xmlns:sbml="http://www.sbml.org/sbml/level3/version1/core">
            <apply>
              <times/>
              <apply>
                <divide/>
                <apply>
                  <times/>
                  <ci> scale_f </ci>
                  <ci> Vmax_v4 </ci>
                </apply>
                <ci> Km_A </ci>
              </apply>
              <apply>
                <minus/>
                <ci> C </ci>
                <apply>
                  <times/>
                  <apply>
                    <divide/>
                    <cn sbml:units="dimensionless" type="integer"> 1 </cn>
                    <ci> Keq_v4 </ci>
                  </apply>
                  <ci> B </ci>
                </apply>
              </apply>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
    </listOfReactions>
  </model>
</sbml>
0u812 commented 9 years ago

Can run these tests via

import roadrunner; import roadrunner.testing; roadrunner.testing.runTester()
kirichoi commented 9 years ago

COPASI reports following flux CC matrix:

image

0u812 commented 9 years ago

Notes: need to bring steady state tolerance and modulation factor out and allow to be configured by config file or Python

0u812 commented 9 years ago

Per @kirichoi we can only reproduce the error with this specific model. He tried others in the Python test suite and the bug did not show up. May need to examine the model in detail.

0u812 commented 9 years ago

COPASI also reports zero for bA, bB, and bC from the first post.

0u812 commented 9 years ago

Here is the summary of this issue:

In summary, I think the calculations are correct after all, although we may want to add some tests for how we handle units at some point. For this issue, I am simply going to bring out the modulation factor so that it can be adjusted from Python.

hsauro commented 9 years ago

Note that the coefficeints are unitless so I'm not sure how units can have an effect. Are you saying however that our coefficeints are out by a factor of ten?

Herbert

On Wed, Sep 16, 2015 at 10:15 AM, Kyle Medley notifications@github.com wrote:

Here is the summary of this issue:

  • Unscaled elasticities and unscaled concentration control coefficients are identical to COPASI down to reordering and scaling by a factor of ten to a power. The digits are exactly the same, just off by orders of magnitude. I suspect this may be due to the use of units in the model.
  • Unscaled flux control coefficients, scaled conc. CC's, and scaled elasticities are identical to COPASI down to reordering.

In summary, I think the calculations are correct after all, although we may want to add some tests for how we handle units at some point. For this issue, I am simply going to bring out the modulation factor so that it can be adjusted from Python.

— Reply to this email directly or view it on GitHub https://github.com/sys-bio/roadrunner/issues/127#issuecomment-140809072.

0u812 commented 9 years ago

The unscaled easlticities and conc. CC's seem to be off by several powers of ten.

On 09/16/2015 10:35 AM, Herbert Sauro wrote:

Note that the coefficeints are unitless so I'm not sure how units can have an effect. Are you saying however that our coefficeints are out by a factor of ten?

Herbert

On Wed, Sep 16, 2015 at 10:15 AM, Kyle Medley notifications@github.com wrote:

Here is the summary of this issue:

  • Unscaled elasticities and unscaled concentration control coefficients are identical to COPASI down to reordering and scaling by a factor of ten to a power. The digits are exactly the same, just off by orders of magnitude. I suspect this may be due to the use of units in the model.
  • Unscaled flux control coefficients, scaled conc. CC's, and scaled elasticities are identical to COPASI down to reordering.

In summary, I think the calculations are correct after all, although we may want to add some tests for how we handle units at some point. For this issue, I am simply going to bring out the modulation factor so that it can be adjusted from Python.

— Reply to this email directly or view it on GitHub

https://github.com/sys-bio/roadrunner/issues/127#issuecomment-140809072.

— Reply to this email directly or view it on GitHub https://github.com/sys-bio/roadrunner/issues/127#issuecomment-140814124.

0u812 commented 9 years ago

RoadRunner:

** Unscaled Elasticities **
                 A,          C,            B,        C_ext,        B_ext,       A_ext
bA [[ -2.28916e-06,          0, -1.85426e-20, -2.12089e-20, -5.02195e-21, 2.28916e-06],
bB  [            0,          0,  2.82528e-07,            0, -2.82528e-07,           0],
bC  [            0, 3.8191e-07,            0,  -3.8191e-07,            0,           0],
v1  [        1e-06,          0,       -1e-07,            0,            0,           0],
v2  [        5e-07,          0,            0,            0,            0,           0],
v3  [            0,      5e-07,            0,            0,            0,           0],
v4  [            0,      5e-07,     -2.5e-07,            0,            0,           0]]

COPASI:

Unscaled elasticities
Unscaled elasticity matrix
Rows:    Reactions (reduced system)
Columns: Species (reduced system)
    A   C   B   A_ext   C_ext   B_ext
(bA)    -2.28916    0   0   2.28916 0   0
(bB)    0   0   0.282528    0   0   -0.282528
(bC)    0   0.38191 0   0   -0.38191    0
(v1)    1   0   -0.1    0   0   0
(v2)    0.5 0   0   0   0   0
(v3)    0   0.5 0   0   0   0
(v4)    0   0.5 -0.25   0   0   0
hsauro commented 9 years ago

If the model is the same (ie same numbers) then there shouldn't be a scaling difference. Mayby copasi is rescaling the concentrations and fluxes so some standard units before doing the computation.

H

On Wed, Sep 16, 2015 at 10:43 AM, Kyle Medley notifications@github.com wrote:

The unscaled easlticities and conc. CC's seem to be off by several powers of ten.

On 09/16/2015 10:35 AM, Herbert Sauro wrote:

Note that the coefficeints are unitless so I'm not sure how units can have an effect. Are you saying however that our coefficeints are out by a factor of ten?

Herbert

On Wed, Sep 16, 2015 at 10:15 AM, Kyle Medley notifications@github.com wrote:

Here is the summary of this issue:

  • Unscaled elasticities and unscaled concentration control coefficients are identical to COPASI down to reordering and scaling by a factor of ten to a power. The digits are exactly the same, just off by orders of magnitude. I suspect this may be due to the use of units in the model.
  • Unscaled flux control coefficients, scaled conc. CC's, and scaled elasticities are identical to COPASI down to reordering.

In summary, I think the calculations are correct after all, although we may want to add some tests for how we handle units at some point. For this issue, I am simply going to bring out the modulation factor so that it can be adjusted from Python.

— Reply to this email directly or view it on GitHub

<https://github.com/sys-bio/roadrunner/issues/127#issuecomment-140809072 .

— Reply to this email directly or view it on GitHub <https://github.com/sys-bio/roadrunner/issues/127#issuecomment-140814124 .

— Reply to this email directly or view it on GitHub https://github.com/sys-bio/roadrunner/issues/127#issuecomment-140816073.

0u812 commented 9 years ago

I made the step size and steady state threshold exposed via the RoadRunner.diffstep and RoadRunner.steadyStateThresh properties in Python. Can we close this and make another issue for the scaling difference that isn't due by COMBINE? I have my Simons app due by the end of the month as well...

hsauro commented 9 years ago

Good you brought out the diffsep and threshold. I won't add another issue, I'll run some test first but what ever happens we won't deal with it now until after COMBINE. I'll close the issue.

On Wed, Sep 16, 2015 at 11:39 AM, Kyle Medley notifications@github.com wrote:

I made the step size and steady state threshold exposed via the RoadRunner.diffstep and RoadRunner.steadyStateThresh properties in Python. Can we close this and make another issue for the scaling difference that isn't due by COMBINE? I have my Simons app due by the end of the month as well...

— Reply to this email directly or view it on GitHub https://github.com/sys-bio/roadrunner/issues/127#issuecomment-140837255.