Tellurium Methods

Installing Packages

Tellurium provides utility methods for installing Python packages from PyPI. These methods simply delegate to pip, and are usually more reliable than running !pip install xyz.

How to install additional packages

If you are using Tellurium notebook or Tellurium Spyder, you can install additional package using installPackage function. In Tellurium Spyder, you can also install packages using included command Prompt. For more information, see Running Command Prompt for Tellurium Spyder.

import tellurium as te
# install cobra (https://github.com/opencobra/cobrapy)
te.installPackage('cobra')
# update cobra to latest version
te.upgradePackage('cobra')
# remove cobra
# te.removePackage('cobra')

Utility Methods

The most useful methods here are the notices routines. Roadrunner will offen issue warning or informational messages. For repeated simulation such messages will clutter up the outputs. noticesOff and noticesOn can be used to turn on an off the messages.

Tellurium’s version can be obtained via te.__version__. .printVersionInfo() also returns information from certain constituent packages.

import tellurium as te

# to get the tellurium version use
print('te.__version__')
print(te.__version__)
# or
print('te.getTelluriumVersion()')
print(te.getTelluriumVersion())

# to print the full version info use
print('-' * 80)
te.printVersionInfo()
print('-' * 80)
te.__version__
2.1.0
te.getTelluriumVersion()
2.1.0
--------------------------------------------------------------------------------
tellurium : 2.1.0
roadrunner : 1.5.1
antimony : 2.9.4
libsbml : 5.15.0
libsedml : 0.4.3
phrasedml : 1.0.9
--------------------------------------------------------------------------------
from builtins import range
# Load SBML file
r = te.loada("""
model test
    J0: X0 -> X1; k1*X0;
    X0 = 10; X1=0;
    k1 = 0.2
end
""")

import matplotlib.pyplot as plt

# Turn off notices so they don't clutter the output
te.noticesOff()
for i in range(0, 20):
    result = r.simulate (0, 10)
    r.reset()
    r.plot(result, loc=None, show=False,
           linewidth=2.0, linestyle='-', color='black', alpha=0.8)
    r.k1 = r.k1 + 0.2
# Turn the notices back on
te.noticesOn()
# create tmp file
import tempfile
ftmp = tempfile.NamedTemporaryFile(suffix=".xml")
# load model
r = te.loada('S1 -> S2; k1*S1; k1 = 0.1; S1 = 10')
# save to file
te.saveToFile(ftmp.name, r.getMatlab())

# or easier via
r.exportToMatlab(ftmp.name)

# load file
sbmlstr = te.readFromFile(ftmp.name)
print('%' + '*'*80)
print('Converted MATLAB code')
print('%' + '*'*80)
print(sbmlstr)
%********************************************************************************
Converted MATLAB code
%********************************************************************************
%  How to use:
%
%  __main takes 3 inputs and returns 3 outputs.
%
%  [t x rInfo] = __main(tspan,solver,options)
%  INPUTS:
%  tspan - the time vector for the simulation. It can contain every time point,
%  or just the start and end (e.g. [0 1 2 3] or [0 100]).
%  solver - the function handle for the odeN solver you wish to use (e.g. @ode23s).
%  options - this is the options structure returned from the MATLAB odeset
%  function used for setting tolerances and other parameters for the solver.
%
%  OUTPUTS:
%  t - the time vector that corresponds with the solution. If tspan only contains
%  the start and end times, t will contain points spaced out by the solver.
%  x - the simulation results.
%  rInfo - a structure containing information about the model. The fields
%  within rInfo are:
%     stoich - the stoichiometry matrix of the model
%     floatingSpecies - a cell array containing floating species name, initial
%     value, and indicator of the units being inconcentration or amount
%     compartments - a cell array containing compartment names and volumes
%     params - a cell array containing parameter names and values
%     boundarySpecies - a cell array containing boundary species name, initial
%     value, and indicator of the units being inconcentration or amount
%     rateRules - a cell array containing the names of variables used in a rate rule
%
%  Sample function call:
%     options = odeset('RelTol',1e-12,'AbsTol',1e-9);
%     [t x rInfo] = __main(linspace(0,100,100),@ode23s,options);
%
function [t x rInfo] = __main(tspan,solver,options)
    % initial conditions
    [x rInfo] = model();

    % initial assignments

    % assignment rules

    % run simulation
    [t x] = feval(solver,@model,tspan,x,options);

    % assignment rules

function [xdot rInfo] = model(time,x)
%  x(1)        S1
%  x(2)        S2

% List of Compartments
vol__default_compartment = 1;               %default_compartment

% Global Parameters
rInfo.g_p1 = 0.1;           % k1

if (nargin == 0)

    % set initial conditions
   xdot(1) = 10*vol__default_compartment;           % S1 = S1 [Concentration]
   xdot(2) = 0*vol__default_compartment;            % S2 = S2 [Concentration]

   % reaction info structure
   rInfo.stoich = [
      -1
      1
   ];

   rInfo.floatingSpecies = {                % Each row: [Species Name, Initial Value, isAmount (1 for amount, 0 for concentration)]
      'S1' , 10, 0
      'S2' , 0, 0
   };

   rInfo.compartments = {           % Each row: [Compartment Name, Value]
      'default_compartment' , 1
   };

   rInfo.params = {         % Each row: [Parameter Name, Value]
      'k1' , 0.1
   };

   rInfo.boundarySpecies = {                % Each row: [Species Name, Initial Value, isAmount (1 for amount, 0 for concentration)]
   };

   rInfo.rateRules = {               % List of variables involved in a rate rule
   };

else

    % calculate rates of change
   R0 = rInfo.g_p1*(x(1));

   xdot = [
      - R0
      + R0
   ];
end;


%listOfSupportedFunctions
function z = pow (x,y)
    z = x^y;


function z = sqr (x)
    z = x*x;


function z = piecewise(varargin)
            numArgs = nargin;
            result = 0;
            foundResult = 0;
            for k=1:2: numArgs-1
                    if varargin{k+1} == 1
                            result = varargin{k};
                            foundResult = 1;
                            break;
                    end
            end
            if foundResult == 0
                    result = varargin{numArgs};
            end
            z = result;


function z = gt(a,b)
   if a > b
      z = 1;
   else
      z = 0;
   end


function z = lt(a,b)
   if a < b
      z = 1;
   else
      z = 0;
   end


function z = geq(a,b)
   if a >= b
      z = 1;
   else
      z = 0;
   end


function z = leq(a,b)
   if a <= b
      z = 1;
   else
      z = 0;
   end


function z = neq(a,b)
   if a ~= b
      z = 1;
   else
      z = 0;
   end


function z = and(varargin)
            result = 1;
            for k=1:nargin
               if varargin{k} ~= 1
                  result = 0;
                  break;
               end
            end
            z = result;


function z = or(varargin)
            result = 0;
            for k=1:nargin
               if varargin{k} ~= 0
                  result = 1;
                  break;
               end
            end
            z = result;


function z = xor(varargin)
            foundZero = 0;
            foundOne = 0;
            for k = 1:nargin
                    if varargin{k} == 0
                       foundZero = 1;
                    else
                       foundOne = 1;
                    end
            end
            if foundZero && foundOne
                    z = 1;
            else
              z = 0;
            end



function z = not(a)
   if a == 1
      z = 0;
   else
      z = 1;
   end


function z = root(a,b)
    z = a^(1/b);

Model Loading

There are a variety of methods to load models into libRoadrunner.

Antimony files can be read with te.loada or te.loadAntimonyModel. For SBML te.loadSBMLModel, for CellML te.loadCellMLModel is used. All the functions accept either model strings or respective model files.

import tellurium as te

# Load an antimony model
ant_model = '''
    S1 -> S2; k1*S1;
    S2 -> S3; k2*S2;

    k1= 0.1; k2 = 0.2;
    S1 = 10; S2 = 0; S3 = 0;
'''
# At the most basic level one can load the SBML model directly using libRoadRunner
print('---  load using roadrunner ---')
import roadrunner
# convert to SBML model
sbml_model = te.antimonyToSBML(ant_model)
r = roadrunner.RoadRunner(sbml_model)
result = r.simulate(0, 10, 100)
r.plot(result)

# The method loada is simply a shortcut to loadAntimonyModel
print('---  load using tellurium ---')
r = te.loada(ant_model)
result = r.simulate (0, 10, 100)
r.plot(result)

# same like
r = te.loadAntimonyModel(ant_model)
---  load using roadrunner ---
_images/tellurium_model_loading_2_11.png
---  load using tellurium ---
_images/tellurium_model_loading_2_31.png

Interconversion Utilities

Use these routines interconvert verious standard formats

Tellurium can convert between Antimony, SBML, and CellML.

import tellurium as te

# antimony model
ant_model = """
    S1 -> S2; k1*S1;
    S2 -> S3; k2*S2;

    k1= 0.1; k2 = 0.2;
    S1 = 10; S2 = 0; S3 = 0;
"""

# convert to SBML
sbml_model = te.antimonyToSBML(ant_model)
print('sbml_model')
print('*'*80)
# print first 10 lines
for line in list(sbml_model.splitlines())[:10]:
    print(line)
print('...')

# convert to CellML
cellml_model = te.antimonyToCellML(ant_model)
print('cellml_model (from Antimony)')
print('*'*80)
# print first 10 lines
for line in list(cellml_model.splitlines())[:10]:
    print(line)
print('...')

# or from the sbml
cellml_model = te.sbmlToCellML(sbml_model)
print('cellml_model (from SBML)')
print('*'*80)
# print first 10 lines
for line in list(cellml_model.splitlines())[:10]:
    print(line)
print('...')
sbml_model
****************************************************************************
<?xml version="1.0" encoding="UTF-8"?>
<!-- Created by libAntimony version v2.9.4 with libSBML version 5.15.0. -->
<sbml xmlns="http://www.sbml.org/sbml/level3/version1/core" level="3" version="1">
  <model id="__main" name="__main">
    <listOfCompartments>
      <compartment sboTerm="SBO:0000410" id="default_compartment" spatialDimensions="3" size="1" constant="true"/>
    </listOfCompartments>
    <listOfSpecies>
      <species id="S1" compartment="default_compartment" initialConcentration="10" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
      <species id="S2" compartment="default_compartment" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
...
cellml_model (from Antimony)
****************************************************************************
<?xml version="1.0"?>
<model xmlns:cellml="http://www.cellml.org/cellml/1.1#" xmlns="http://www.cellml.org/cellml/1.1#" name="__main">
<component name="__main">
<variable initial_value="10" name="S1" units="dimensionless"/>
<variable initial_value="0" name="S2" units="dimensionless"/>
<variable initial_value="0.1" name="k1" units="dimensionless"/>
<variable name="_J0" units="dimensionless"/>
<variable initial_value="0" name="S3" units="dimensionless"/>
<variable initial_value="0.2" name="k2" units="dimensionless"/>
<variable name="_J1" units="dimensionless"/>
...
cellml_model (from SBML)
****************************************************************************
<?xml version="1.0"?>
<model xmlns:cellml="http://www.cellml.org/cellml/1.1#" xmlns="http://www.cellml.org/cellml/1.1#" name="__main">
<component name="__main">
<variable initial_value="10" name="S1" units="dimensionless"/>
<variable initial_value="0" name="S2" units="dimensionless"/>
<variable initial_value="0" name="S3" units="dimensionless"/>
<variable initial_value="0.1" name="k1" units="dimensionless"/>
<variable initial_value="0.2" name="k2" units="dimensionless"/>
<variable name="_J0" units="dimensionless"/>
<variable name="_J1" units="dimensionless"/>
...

Export Utilities

Use these routines to convert the current model state into other formats, like Matlab, CellML, Antimony and SBML.

Given a RoadRunner instance, you can get an SBML representation of the current state of the model using getCurrentSBML. You can also get the initial SBML from when the model was loaded using getSBML. Finally, exportToSBML can be used to export the current model state to a file.

import tellurium as te
import tempfile

# load model
r = te.loada('S1 -> S2; k1*S1; k1 = 0.1; S1 = 10')
# file for export
f_sbml = tempfile.NamedTemporaryFile(suffix=".xml")

# export current model state
r.exportToSBML(f_sbml.name)

# to export the initial state when the model was loaded
# set the current argument to False
r.exportToSBML(f_sbml.name, current=False)

# The string representations of the current model are available via
str_sbml = r.getCurrentSBML()

# and of the initial state when the model was loaded via
str_sbml = r.getSBML()
print(str_sbml)
<?xml version="1.0" encoding="UTF-8"?>
<!-- Created by libAntimony version v2.9.4 with libSBML version 5.15.0. -->
<sbml xmlns="http://www.sbml.org/sbml/level3/version1/core" level="3" version="1">
  <model id="__main" name="__main">
    <listOfCompartments>
      <compartment sboTerm="SBO:0000410" id="default_compartment" spatialDimensions="3" size="1" constant="true"/>
    </listOfCompartments>
    <listOfSpecies>
      <species id="S1" compartment="default_compartment" initialConcentration="10" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
      <species id="S2" compartment="default_compartment" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
    </listOfSpecies>
    <listOfParameters>
      <parameter id="k1" value="0.1" constant="true"/>
    </listOfParameters>
    <listOfReactions>
      <reaction id="_J0" reversible="true" 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> k1 </ci>
              <ci> S1 </ci>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
    </listOfReactions>
  </model>
</sbml>

Similar to the SBML functions above, you can also use the functions getCurrentAntimony and exportToAntimony to get or export the current Antimony representation.

import tellurium as te
import tempfile

# load model
r = te.loada('S1 -> S2; k1*S1; k1 = 0.1; S1 = 10')
# file for export
f_antimony = tempfile.NamedTemporaryFile(suffix=".txt")

# export current model state
r.exportToAntimony(f_antimony.name)

# to export the initial state when the model was loaded
# set the current argument to False
r.exportToAntimony(f_antimony.name, current=False)

# The string representations of the current model are available via
str_antimony = r.getCurrentAntimony()

# and of the initial state when the model was loaded via
str_antimony = r.getAntimony()
print(str_antimony)
// Created by libAntimony v2.9.4
// Compartments and Species:
species S1, S2;

// Reactions:
_J0: S1 -> S2; k1*S1;

// Species initializations:
S1 = 10;
S2 = ;

// Variable initializations:
k1 = 0.1;

// Other declarations:
const k1;

Tellurium also has functions for exporting the current model state to CellML. These functionalities rely on using Antimony to perform the conversion.

import tellurium as te
import tempfile

# load model
r = te.loada('S1 -> S2; k1*S1; k1 = 0.1; S1 = 10')
# file for export
f_cellml = tempfile.NamedTemporaryFile(suffix=".cellml")

# export current model state
r.exportToCellML(f_cellml.name)

# to export the initial state when the model was loaded
# set the current argument to False
r.exportToCellML(f_cellml.name, current=False)

# The string representations of the current model are available via
str_cellml = r.getCurrentCellML()

# and of the initial state when the model was loaded via
str_cellml = r.getCellML()
print(str_cellml)
<?xml version="1.0"?>
<model xmlns:cellml="http://www.cellml.org/cellml/1.1#" xmlns="http://www.cellml.org/cellml/1.1#" name="__main">
<component name="__main">
<variable initial_value="10" name="S1" units="dimensionless"/>
<variable name="S2" units="dimensionless"/>
<variable initial_value="0.1" name="k1" units="dimensionless"/>
<variable name="_J0" units="dimensionless"/>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<eq/>
<ci>_J0</ci>
<apply>
<times/>
<ci>k1</ci>
<ci>S1</ci>
</apply>
</apply>
</math>
<variable name="time" units="dimensionless"/>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<eq/>
<apply>
<diff/>
<bvar>
<ci>time</ci>
</bvar>
<ci>S1</ci>
</apply>
<apply>
<minus/>
<ci>_J0</ci>
</apply>
</apply>
</math>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<eq/>
<apply>
<diff/>
<bvar>
<ci>time</ci>
</bvar>
<ci>S2</ci>
</apply>
<ci>_J0</ci>
</apply>
</math>
</component>
<group>
<relationship_ref relationship="encapsulation"/>
<component_ref component="__main"/>
</group>
</model>

To export the current model state to MATLAB, use getCurrentMatlab.

import tellurium as te
import tempfile

# load model
r = te.loada('S1 -> S2; k1*S1; k1 = 0.1; S1 = 10')
# file for export
f_matlab = tempfile.NamedTemporaryFile(suffix=".m")

# export current model state
r.exportToMatlab(f_matlab.name)

# to export the initial state when the model was loaded
# set the current argument to False
r.exportToMatlab(f_matlab.name, current=False)

# The string representations of the current model are available via
str_matlab = r.getCurrentMatlab()

# and of the initial state when the model was loaded via
str_matlab = r.getMatlab()
print(str_matlab)
%  How to use:
%
%  __main takes 3 inputs and returns 3 outputs.
%
%  [t x rInfo] = __main(tspan,solver,options)
%  INPUTS:
%  tspan - the time vector for the simulation. It can contain every time point,
%  or just the start and end (e.g. [0 1 2 3] or [0 100]).
%  solver - the function handle for the odeN solver you wish to use (e.g. @ode23s).
%  options - this is the options structure returned from the MATLAB odeset
%  function used for setting tolerances and other parameters for the solver.
%
%  OUTPUTS:
%  t - the time vector that corresponds with the solution. If tspan only contains
%  the start and end times, t will contain points spaced out by the solver.
%  x - the simulation results.
%  rInfo - a structure containing information about the model. The fields
%  within rInfo are:
%     stoich - the stoichiometry matrix of the model
%     floatingSpecies - a cell array containing floating species name, initial
%     value, and indicator of the units being inconcentration or amount
%     compartments - a cell array containing compartment names and volumes
%     params - a cell array containing parameter names and values
%     boundarySpecies - a cell array containing boundary species name, initial
%     value, and indicator of the units being inconcentration or amount
%     rateRules - a cell array containing the names of variables used in a rate rule
%
%  Sample function call:
%     options = odeset('RelTol',1e-12,'AbsTol',1e-9);
%     [t x rInfo] = __main(linspace(0,100,100),@ode23s,options);
%
function [t x rInfo] = __main(tspan,solver,options)
    % initial conditions
    [x rInfo] = model();

    % initial assignments

    % assignment rules

    % run simulation
    [t x] = feval(solver,@model,tspan,x,options);

    % assignment rules

function [xdot rInfo] = model(time,x)
%  x(1)        S1
%  x(2)        S2

% List of Compartments
vol__default_compartment = 1;               %default_compartment

% Global Parameters
rInfo.g_p1 = 0.1;           % k1

if (nargin == 0)

    % set initial conditions
   xdot(1) = 10*vol__default_compartment;           % S1 = S1 [Concentration]
   xdot(2) = 0*vol__default_compartment;            % S2 = S2 [Concentration]

   % reaction info structure
   rInfo.stoich = [
      -1
      1
   ];

   rInfo.floatingSpecies = {                % Each row: [Species Name, Initial Value, isAmount (1 for amount, 0 for concentration)]
      'S1' , 10, 0
      'S2' , 0, 0
   };

   rInfo.compartments = {           % Each row: [Compartment Name, Value]
      'default_compartment' , 1
   };

   rInfo.params = {         % Each row: [Parameter Name, Value]
      'k1' , 0.1
   };

   rInfo.boundarySpecies = {                % Each row: [Species Name, Initial Value, isAmount (1 for amount, 0 for concentration)]
   };

   rInfo.rateRules = {               % List of variables involved in a rate rule
   };

else

    % calculate rates of change
   R0 = rInfo.g_p1*(x(1));

   xdot = [
      - R0
      + R0
   ];
end;


%listOfSupportedFunctions
function z = pow (x,y)
    z = x^y;


function z = sqr (x)
    z = x*x;


function z = piecewise(varargin)
            numArgs = nargin;
            result = 0;
            foundResult = 0;
            for k=1:2: numArgs-1
                    if varargin{k+1} == 1
                            result = varargin{k};
                            foundResult = 1;
                            break;
                    end
            end
            if foundResult == 0
                    result = varargin{numArgs};
            end
            z = result;


function z = gt(a,b)
   if a > b
      z = 1;
   else
      z = 0;
   end


function z = lt(a,b)
   if a < b
      z = 1;
   else
      z = 0;
   end


function z = geq(a,b)
   if a >= b
      z = 1;
   else
      z = 0;
   end


function z = leq(a,b)
   if a <= b
      z = 1;
   else
      z = 0;
   end


function z = neq(a,b)
   if a ~= b
      z = 1;
   else
      z = 0;
   end


function z = and(varargin)
            result = 1;
            for k=1:nargin
               if varargin{k} ~= 1
                  result = 0;
                  break;
               end
            end
            z = result;


function z = or(varargin)
            result = 0;
            for k=1:nargin
               if varargin{k} ~= 0
                  result = 1;
                  break;
               end
            end
            z = result;


function z = xor(varargin)
            foundZero = 0;
            foundOne = 0;
            for k = 1:nargin
                    if varargin{k} == 0
                       foundZero = 1;
                    else
                       foundOne = 1;
                    end
            end
            if foundZero && foundOne
                    z = 1;
            else
              z = 0;
            end



function z = not(a)
   if a == 1
      z = 0;
   else
      z = 1;
   end


function z = root(a,b)
    z = a^(1/b);

The above examples rely on Antimony as in intermediary between formats. You can use this functionality directly using e.g. antimony.getCellMLString. A comprehensive set of functions can be found in the Antimony API documentation.

import antimony
antimony.loadAntimonyString('''S1 -> S2; k1*S1; k1 = 0.1; S1 = 10''')
ant_str = antimony.getCellMLString(antimony.getMainModuleName())
print(ant_str)
<?xml version="1.0"?>
<model xmlns:cellml="http://www.cellml.org/cellml/1.1#" xmlns="http://www.cellml.org/cellml/1.1#" name="__main">
<component name="__main">
<variable initial_value="10" name="S1" units="dimensionless"/>
<variable name="S2" units="dimensionless"/>
<variable initial_value="0.1" name="k1" units="dimensionless"/>
<variable name="_J0" units="dimensionless"/>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<eq/>
<ci>_J0</ci>
<apply>
<times/>
<ci>k1</ci>
<ci>S1</ci>
</apply>
</apply>
</math>
<variable name="time" units="dimensionless"/>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<eq/>
<apply>
<diff/>
<bvar>
<ci>time</ci>
</bvar>
<ci>S1</ci>
</apply>
<apply>
<minus/>
<ci>_J0</ci>
</apply>
</apply>
</math>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<eq/>
<apply>
<diff/>
<bvar>
<ci>time</ci>
</bvar>
<ci>S2</ci>
</apply>
<ci>_J0</ci>
</apply>
</math>
</component>
<group>
<relationship_ref relationship="encapsulation"/>
<component_ref component="__main"/>
</group>
</model>

Stochastic Simulation

Use these routines to carry out Gillespie style stochastic simulations.

Stochastic simulations can be run by changing the current integrator type to ‘gillespie’ or by using the r.gillespie function.

import tellurium as te
import numpy as np

r = te.loada('S1 -> S2; k1*S1; k1 = 0.1; S1 = 40')
r.integrator = 'gillespie'
r.integrator.seed = 1234

results = []
for k in range(1, 50):
    r.reset()
    s = r.simulate(0, 40)
    results.append(s)
    r.plot(s, show=False, alpha=0.7)
te.show()
_images/tellurium_stochastic_2_01.png

Setting the identical seed for all repeats results in identical traces in each simulation.

results = []
for k in range(1, 20):
    r.reset()
    r.setSeed(123456)
    s = r.simulate(0, 40)
    results.append(s)
    r.plot(s, show=False, loc=None, color='black', alpha=0.7)
te.show()
_images/tellurium_stochastic_4_01.png

You can combine two timecourse simulations and change e.g. parameter values in between each simulation. The gillespie method simulates up to the given end time 10, after which you can make arbitrary changes to the model, then simulate again.

When using the r.plot function, you can pass the parameter labels, which controls the names that will be used in the figure legend, and tag, which ensures that traces with the same tag will be drawn with the same color (each species within each trace will be plotted in its own color, but these colors will match trace to trace).

import tellurium as te

r = te.loada('S1 -> S2; k1*S1; k1 = 0.02; S1 = 100')
r.setSeed(1234)
for k in range(1, 20):
    r.resetToOrigin()
    res1 = r.gillespie(0, 10)
    r.plot(res1, show=False) # plot first half of data

    # change in parameter after the first half of the simulation
    # We could have also used an Event in the antimony model,
    # which are described further in the Antimony Reference section
    r.k1 = r.k1*20
    res2 = r.gillespie (10, 20)

    r.plot(res2, show=False) # plot second half of data

te.show()
_images/tellurium_stochastic_6_01.png

Math

Use these routines to perform various calculations.

ODE Extraction Methods

Routines to extract ODEs.

Plotting

Tellurium has a plotting engine which can target either Plotly (when used in a notebook environment) or Matplotlib. To specify which engine to use, use te.setDefaultPlottingEngine().

NOTE: When loading a model with r = te.loada(‘antimony_string’) and calling r.plot(), it is the below tellerium.ExtendedRoadRunner.plot() method below that is called not te.plot().

The function tellurium.plotArray assumes that the first column in the array is the x-axis and the second and subsequent columns represent curves on the y-axis.

NOTE: When loading a model with r = te.loada(‘antimony_string’) and calling r.plot(), it is the above tellerium.ExtendedRoadRunner.plot() method below that is called not te.plot().

Add plot elements

Example showing how to embelish a graph - change title, axes labels, set axis limit. Example also uses an event to pulse S1.

import tellurium as te, roadrunner

r = te.loada ('''
   $Xo -> S1; k1*Xo;
   S1 -> $X1; k2*S1;

   k1 = 0.2; k2 = 0.4; Xo = 1; S1 = 0.5;
   at (time > 20): S1 = S1 + 0.35
''')

# Simulate the first part up to 20 time units
m = r.simulate (0, 50, 100, ["time", "S1"])
                                                            # using latex syntax to render math
r.plot(m, ylim=(0.,1.), xtitle='Time', ytitle='Concentration', title='My First Plot ($y = x^2$)')
_images/tellurium_examples_9_01.png

Saving plots

To save a plot, use r.plot and the savefig parameter. Use dpi to specify image quality. Pass in the save location along with the image name.

import tellurium as te, os
r = te.loada('S1 -> S2; k1*S1; k1 = 0.1; S1 = 10')
result = r.simulate(0, 50, 100)
currentDir = os.getcwd() # gets the current directory
r.plot(title='My plot', xtitle='Time', ytitle='Concentration', dpi=150,
       savefig=currentDir + '\\test.png')  # save image to current directory as "test.png"

The path can be specified as a written out string. The plot can also be saved as a pdf instead of png.

savefig='C:\\Tellurium-Winpython-3.6\\settings\\.spyder-py3\\test.pdf'

Logarithmic axis

The axis scale can be adapted with the xscale and yscale settings.

import tellurium as te

r = te.loadTestModel('feedback.xml')
r.integrator.variable_step_size = True
s = r.simulate(0, 50)
r.plot(s, logx=True, xlim=[10E-4, 10E2],
      title="Logarithmic x-Axis with grid", ylabel="concentration");
_images/tellurium_plotting_4_01.png

Plotting multiple simulations

All plotting is done via the r.plot or te.plotArray functions. To plot multiple curves in one figure use the show=False setting.

import tellurium as te

import numpy as np
import matplotlib.pylab as plt

# Load a model and carry out a simulation generating 100 points
r = te.loada ('S1 -> S2; k1*S1; k1 = 0.1; S1 = 10')
r.draw(width=100)

# get colormap
# Colormap instances are used to convert data values (floats) from the interval [0, 1]
cmap = plt.get_cmap('Blues')

k1_values = np.linspace(start=0.1, stop=1.5, num=15)
max_k1 = max(k1_values)
for k, value in enumerate(k1_values):
    r.reset()
    r.k1 = value
    s = r.simulate(0, 30, 100)

    color = cmap((value+max_k1)/(2*max_k1))
    # use show=False to plot multiple curves in the same figure
    r.plot(s, show=False, title="Parameter variation k1", xtitle="time", ytitle="concentration",
          xlim=[-1, 31], ylim=[-0.1, 11])

te.show()

print('Reference Simulation: k1 = {}'.format(r.k1))
print('Parameter variation: k1 = {}'.format(k1_values))
_images/tellurium_plotting_2_01.png _images/tellurium_plotting_2_11.png
Reference Simulation: k1 = 1.5
Parameter variation: k1 = [0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.  1.1 1.2 1.3 1.4 1.5]

Using Tags and Names

Tags can be used to coordinate the color, opacity, and legend names between several sets of data. This can be used to highlight certain features that these datasets have in common. Names allow you to give a more meaningful description of the data in the legend.

import tellurium as te
import numpy as np

for i in range(1, 10):
    x = np.linspace(0, 10, num = 10)
    y = i*x**2 + 10*i

    if i % 2 == 0:
        next_tag = "positive slope"
    else:
        next_tag = "negative slope"
        y = -1*y

    next_name = next_tag + " (i = " + str(i) + ")"
    te.plot(x, y, show = False, tag = next_tag, name = next_name)

te.show()
_images/tellurium_plotting_3_01.png

Note that only two items show up in the legend, one for each tag used. In this case, the name found in the legend will match the name of the last set of data plotted using that specific tag. The color and opacity for each tagged groups will also be chosen from the last dataset inputted with that given tag.

Subplots

te.plotArray can be used in conjunction with matplotlib functions to create subplots.

import tellurium as te
import numpy as np
import matplotlib.pylab as plt

r = te.loada ('S1 -> S2; k1*S1; k1 = 0.1; S1 = 20')
r.setIntegrator('gillespie')
r.integrator.seed = '1234'
kValues = np.linspace(0.1, 0.9, num=9) # generate k1 values

plt.gcf().set_size_inches(10, 10) # size of figure
plt.subplots_adjust(wspace=0.4, hspace=0.4) # adjust the space between subplots
plt.suptitle('Variation in k1 value', fontsize=16) # main title

for i in range(1, len(kValues) + 1):
    r.k1 = kValues[i - 1]
    # designates number of subplots (row, col) and spot to plot next
    plt.subplot(3, 3, i)
    for j in range(1, 30):
        r.reset()
        s = r.simulate(0, 10)
        t = "k1 = " + '{:.1f}'.format(kValues[i - 1])
        # plot each subplot, use show=False to save multiple traces
        te.plotArray(s, show=False, title=t, xlabel='Time',
                     ylabel='Concentration', alpha=0.7)
_images/tellurium_plotting_1_01.png

External Plotting

For those more familiar with plotting in Python, other libraries such as matplotlib.pylab offer a wider range of plotting options. To use these external libraries, extract the simulation timecourse data returned from r.simulate. Data is returned in the form of a dictionary/NamedArray, so specific elements can easily be extracted using the species name as the key.

import tellurium as te
import matplotlib.pylab as plt

antimonyString = ('''
model feedback()
// Reactions:
J0: Nan1 + Mol -> Nan1Mol; (K1*Nan1*Mol);
J1: Nan1Mol -> Nan1 + Mol; (K_1*Nan1Mol);
J2: Nan1Mol + Nan2 -> Nan1MolNan2; (K2*Nan1Mol*Nan2)
J3: Nan1MolNan2 + GeneOff -> GeneOn; (K3*Nan1MolNan2*GeneOff);
J4: GeneOn -> Nan1MolNan2 + GeneOff; (K_3*GeneOn);

// Species initializations:
Nan1 = 0.0001692; Mol = 0.0001692/2; Nan2 = 0.0001692; Nan1Mol = 0;
Nan1MolNan2 = 0; GeneOff = 5*10^-5; GeneOn = 0;

// Variable initialization:
K1 = 6.1*10^5; K_1 = 8*10^-5; K2 = 3.3*10^5; K_2 = 5.7*10^-8;  K3 = 1*10^5; K_3 = 0;
end''')

r = te.loada(antimonyString)
results = r.simulate(0,0.5,1000)
r.plot()

plt.figure(figsize=(30,10));
plt.rc('font', size=30);

plt.subplot(1,2,1);
plt.plot(results['time'], results['[Nan2]'], 'r', results['time'], results['[Nan1MolNan2]'], 'b');
plt.legend({'Nan2', 'Nan1MolNan2'});

plt.subplot(1,2,2);
plt.plot(results['time'], results['[GeneOff]'], 'r', results['time'], results['[GeneOn]'], 'b');
plt.legend({'GeneOff', 'GeneOn'});
_images/tellurium_plotting_extendedplotting1.png

Note that we can extract all the time course data for a specific species such as Nan2 by calling results['[Nan2]']. The extract brackets [ ] around Nan2 may or may not be required depending on if the units are in terms of concentration or just a count. To check, simply print out results and you can see the names of each species.

Draw diagram

This example shows how to draw a network diagram, requires graphviz.

import tellurium as te


r = te.loada('''
model feedback()
   // Reactions:http://localhost:8888/notebooks/core/tellurium_export.ipynb#
   J0: $X0 -> S1; (VM1 * (X0 - S1/Keq1))/(1 + X0 + S1 +   S4^h);
   J1: S1 -> S2; (10 * S1 - 2 * S2) / (1 + S1 + S2);
   J2: S2 -> S3; (10 * S2 - 2 * S3) / (1 + S2 + S3);
   J3: S3 -> S4; (10 * S3 - 2 * S4) / (1 + S3 + S4);
   J4: S4 -> $X1; (V4 * S4) / (KS4 + S4);

  // Species initializations:
  S1 = 0; S2 = 0; S3 = 0;
  S4 = 0; X0 = 10; X1 = 0;

  // Variable initialization:
  VM1 = 10; Keq1 = 10; h = 10; V4 = 2.5; KS4 = 0.5;
end''')

# simulate using variable step size
r.integrator.setValue('variable_step_size', True)
s = r.simulate(0, 50)
# draw the diagram
r.draw(width=200)
# and the plot
r.plot(s, title="Feedback Oscillations", ylabel="concentration", alpha=0.9);
_images/tellurium_plotting_6_01.png _images/tellurium_plotting_6_11.png

Parameter Scans

To study the consequences of varying a specific parameter value or initial concentration on a simulation, iteratively adjust the given parameter over a range of values of interest and re-run the simulation. Using the show parameter and te.show we can plot all these simulations on a single figure.

import tellurium as te
import roadrunner
import numpy as np

r = te.loada("""
     $Xo -> A; k1*Xo;
      A -> B; kf*A - kr*B;
      B -> ; k2*B;

      Xo = 5
      k1 = 0.1; k2 = 0.5;
      kf = 0.3; kr = 0.4
""")

for Xo in np.arange(1.0, 10, 1):
    r.reset()
    r.Xo = Xo
    m = r.simulate (0, 50, 100, ['time', 'A'])
    te.plotArray (m, show=False, labels=['Xo='+str(Xo)], resetColorCycle=False)
te.show()
_images/tellurium_plotting_parameter_scans1.png

Parameter Uncertainty Modeling

In most systems, some parameters are more sensitve to perturbations than others. When studying these systems, it is important to understand which parameters are highly sensitive, as errors (i.e. measurement error) introduced to these variables can create drastic differences between experimental and simulated results. To study the sensitivity of these parameters, we can sweep over a range of values as we did in the parameter scan example above. These ranges represent our uncertainty in the value of the parameter, and those parameters that create highly variable results in some measure of an output variable are deemed to be sensitive.

import numpy as np
import tellurium as te
import roadrunner
import antimony
import matplotlib.pyplot as plt
import math

antimonyString = ('''
model feedback()
  // Reactions:
  J0: Nan1 + Mol -> Nan1Mol; (K1*Nan1*Mol);
  J1: Nan1Mol -> Nan1 + Mol; (K_1*Nan1Mol);
  J2: Nan1Mol + Nan2 -> Nan1MolNan2; (K2*Nan1Mol*Nan2)
  J3: Nan1MolNan2 + GeneOff -> GeneOn; (K3*Nan1MolNan2*GeneOff);
  J4: GeneOn -> Nan1MolNan2 + GeneOff; (K_3*GeneOn);

  // Species initializations:
  Nan1 = 0.0001692; Mol = 0.0001692/2; Nan2 = 0.0001692; Nan1Mol = 0;
  Nan1MolNan2 = 0; GeneOff = 5*10^-5; GeneOn = 0;

  // Variable initialization:
  K1 = 6.1*10^5; K_1 = 8*10^-5; K2 = 3.3*10^5; K_2 = 5.7*10^-8;  K3 = 1*10^5; K_3 = 0;
end''')

r = te.loada (model.antimonyString)

def plot_param_uncertainty(model, startVal, name, num_sims):
    stdDev = 0.6

    # assumes initial parameter estimate as mean and iterates 60% above and below.
    vals = np.linspace((1-stdDev)*startVal, (1+stdDev)*startVal, 100)
    for val in vals:
        r.resetToOrigin()
        exec("r.%s = %f" % (name, val))
        result = r.simulate(0,0.5,1000, selections = ['time', 'GeneOn'])
        plt.plot(result[:,0],result[:,1])
        plt.title("uncertainty in " + name)
    plt.legend(["GeneOn"])
    plt.xlabel("Time (hours)")
    plt.ylabel("Concentration")

startVals = r.getGlobalParameterValues();
names = list(enumerate([x for x in r.getGlobalParameterIds() if ("K" in x or "k" in x)]));

n = len(names) + 1;
dim = math.ceil(math.sqrt(n))
for i,next_param in enumerate(names):
    plt.subplot(dim,dim,i+1)
    plot_param_uncertainty(r, startVals[next_param[0]], next_param[1], 100)

plt.tight_layout()
plt.show()
_images/tellurium_plotting_parameter_uncertainty1.png

In the above code, the exec command is used to set the model parameters to their given value (i.e. r.K1 = 1.5) and the code sweeps through all the given parameters of interests (names). Above, we see that the K3 parameter produces the widest distribution of outcomes, and is thus the most sensitive under the given model, taking into account its assumptions and approximate parameter values. On the other hand, variations in K_1, K1, and K_3 seem to have very little effect on the outcome of the system.

Model Reset

The reset function of a RoadRunner instance reset the system variables (usually species concentrations) to their respective initial values. resetAll resets variables to their CURRENT initial as well as resets parameters. resetToOrigin completely resets the model.

import tellurium as te

r = te.loada ('S1 -> S2; k1*S1; k1 = 0.1; S1 = 10')
r.integrator.setValue('variable_step_size', True)

# simulate model
sim1 = r.simulate(0, 5)
print('*** sim1 ***')
r.plot(sim1)

# continue from end concentration of sim1
r.k1 = 2.0
sim2 = r.simulate(0, 5)
print('-- sim2 --')
print('continue simulation from final concentrations with changed parameter')
r.plot(sim2)

# Reset initial concentrations, does not affect the changed parameter
r.reset()
sim3 = r.simulate(0, 5)
print('-- sim3 --')
print('reset initial concentrations but keep changed parameter')
r.plot(sim3)

# Change CURRENT initial of k1, resetAll clears parameter but
# resets to CURRENT initial
r.setValue('init(k1)', 0.3)
r.resetAll()
sim4 = r.simulate(0, 5)
print('-- sim4 --')
print('reset to CURRENT initial of k1, reset to initial parameters')
print('k1 = ' + str(r.k1))
r.plot(sim4)

# Reset model to the state it was loaded
r.resetToOrigin()
sim5 = r.simulate(0, 5)
print('-- sim5 --')
print('reset all to origin')
r.plot(sim5);
* sim1 *
_images/tellurium_reset_2_11.png
-- sim2 --
continue simulation from final concentrations with changed parameter
_images/tellurium_reset_2_31.png
-- sim3 --
reset initial concentrations but keep changed parameter
_images/tellurium_reset_2_51.png
-- sim4 --
reset to CURRENT initial of k1, reset to initial parameters
k1 = 0.3
_images/tellurium_reset_2_61.png
-- sim5 --
reset all to origin
_images/tellurium_reset_2_71.png

jarnac Short-cuts

Routines to support the Jarnac compatibility layer

Test Models

RoadRunner has built into it a number of predefined models that can be use to easily try and test tellurium.

Test models

import tellurium as te

# To get the builtin models use listTestModels
print(te.listTestModels())
['EcoliCore.xml', 'linearPathwayClosed.xml', 'test_1.xml', 'linearPathwayOpen.xml', 'feedback.xml']

Load test model

# To load one of the test models use loadTestModel:
# r = te.loadTestModel('feedback.xml')
# result = r.simulate (0, 10, 100)
# r.plot (result)

# If you need to obtain the SBML for the test model, use getTestModel
sbml = te.getTestModel('feedback.xml')

# To look at one of the test model in Antimony form:
ant = te.sbmlToAntimony(te.getTestModel('feedback.xml'))
print(ant)
// Created by libAntimony v2.9.4
model *feedback()

  // Compartments and Species:
  compartment compartment_;
  species S1 in compartment_, S2 in compartment_, S3 in compartment_, S4 in compartment_;
  species $X0 in compartment_, $X1 in compartment_;

  // Reactions:
  J0: $X0 => S1; J0_VM1*(X0 - S1/J0_Keq1)/(1 + X0 + S1 + S4^J0_h);
  J1: S1 => S2; (10*S1 - 2*S2)/(1 + S1 + S2);
  J2: S2 => S3; (10*S2 - 2*S3)/(1 + S2 + S3);
  J3: S3 => S4; (10*S3 - 2*S4)/(1 + S3 + S4);
  J4: S4 => $X1; J4_V4*S4/(J4_KS4 + S4);

  // Species initializations:
  S1 = 0;
  S2 = 0;
  S3 = 0;
  S4 = 0;
  X0 = 10;
  X1 = 0;

  // Compartment initializations:
  compartment_ = 1;

  // Variable initializations:
  J0_VM1 = 10;
  J0_Keq1 = 10;
  J0_h = 10;
  J4_V4 = 2.5;
  J4_KS4 = 0.5;

  // Other declarations:
  const compartment_, J0_VM1, J0_Keq1, J0_h, J4_V4, J4_KS4;
end

Running external tools

Routines to run external tools.

Model Methods

Routines flattened from model, saves typing and easier for finding the methods