SBML Solver¶
When you study biology, sooner or later, you encounter pathway diagrams,
gene expression networks, Physiologicaly Based Pharmacokinetics (PBPK)
whole body diagrams, etc… Often, these can be mathematically represented
in the form of Ordinary Differential Equations (ODEs). There are many
ODE solvers available and you can write your own. However the solution
we like most is called SBML Solver
. Before going any further let us
explain briefly what SBML itself is. SBML stands for Systems Biology Markup Language.
It was proposed around year 2000 by few scientists from
Caltech (Mike Hucka, Herbert Sauro, Andrew Finney). According to
Wikipedia, SBML is a representation format, based on XML, for
communicating and storing computational models of biological processes.
In practice SBML focuses on reaction kinetics models but can also be
used to code these models that can be described in the form of ODEs such
as e.g. PBPK, population models etc…
Being a multi-cell modeling platform CC3D allows users to associate
multiple SBML model solvers with a single cell or create “free floating”
SBML model solvers. The CC3D Python syntax that deals with the SBML
models is referred to as SBML Solver. Internally SBML Solver relies on a
C++ RoadRunnerLib developed by Herbert Sauro team. RoadRunnerLib in turn
is based on the C# code written by Frank Bergmann. CC3D uses
RoadRunnerLib as the engine to solve systems of ODEs. All SBML Solver
functionality is available via SteppableBasePy member functions.
Twedit++ provides nice shortcuts that help users write valid code
tapping into SBML Solver functionality. See CC3DPython->SBML Solver
menu
for options available.
Let us look at the example steppable that uses SBML Solver:
from cc3d.core.PySteppables import *
class SBMLSolverSteppable(SteppableBasePy):
def __init__(self, frequency=1):
SteppableBasePy.__init__(self, frequency)
def start(self):
# adding options that setup SBML solver integrator
# these are optional but useful when encountering integration instabilities
options = {'relative': 1e-10, 'absolute': 1e-12}
self.set_sbml_global_options(options)
model_file = 'Simulation/test_1.xml'
initial_conditions = {}
initial_conditions['S1'] = 0.00020
initial_conditions['S2'] = 0.000002
self.add_sbml_to_cell_ids(model_file=model_file, model_name='dp',
cell_ids=list(range(1, 11)), step_size=0.5, initial_conditions=initial_conditions)
self.add_free_floating_sbml(model_file=model_file, model_name='Medium_dp', step_size=0.5,
initial_conditions=initial_conditions)
self.add_free_floating_sbml(model_file=model_file, model_name='Medium_dp1', step_size=0.5,
initial_conditions=initial_conditions)
self.add_free_floating_sbml(model_file=model_file, model_name='Medium_dp2')
self.add_free_floating_sbml(model_file=model_file, model_name='Medium_dp3')
self.add_free_floating_sbml(model_file=model_file, model_name='Medium_dp4')
cell_20 = self.fetch_cell_by_id(20)
self.add_sbml_to_cell(model_file=model_file, model_name='dp', cell=cell_20)
def step(self, mcs):
self.timestep_sbml()
cell_20 = self.fetch_cell_by_id(20)
print('cell_20, dp=', cell_20.sbml.dp.values())
print('Free Floating Medium_dp2', self.sbml.Medium_dp2.values())
if mcs == 3:
Medium_dp2 = self.sbml.Medium_dp2
Medium_dp2['S1'] = 10
Medium_dp2['S2'] = 0.5
if mcs == 5:
self.delete_sbml_from_cell_ids(model_name='dp', cell_ids=list(range(1, 11)))
if mcs == 7:
cell_25 = self.fetch_cell_by_id(25)
self.copy_sbml_simulators(from_cell=cell_20, to_cell=cell_25)
In the start function we specify path to the SBML model (here we use
partial path Simulation/test_1.xml
because test_1.xml
is in our CC3D
Simulation project directory) and also create python dictionary that has
initial conditions for the SBML model. This particular model has two
floating species : S1
and S2
and our dictionary – initialConditions
stores the initial concentration of these species to 0.0002 and 0.000002
respectively:
model_file = 'Simulation/test_1.xml'
initial_conditions = {}
initial_conditions['S1'] = 0.00020
initial_conditions['S2'] = 0.000002
Note
We can initialize each SBML Solver using different initial conditions. When we forget to specify initial conditions the SBML code usually has initial conditions defined and they will be used as starting values.
Before we discuss add_sbml_to_cell_ids
function let us focus on statements
that open the start function:
options = {'relative': 1e-10, 'absolute': 1e-12}
self.set_sbml_global_options(options)
We set here SBML integrator options. These statements are optional,
however when your SBML model crashes with e.g. CVODE error, it often
means that your numerical tolerances (relative and absolute) or number
of integration steps in each integration interval (steps) should be
changed. Additionally you may want to enable stiff ODE solver by setting
stiff to True
:
options = {'relative': 1e-10, 'absolute': 1e-12, 'stiff': False}
self.set_sbml_global_options(options)
After defining options dictionary we inform CC3D to use these settings
. We do it by using as shown above. A thing to remember that new options
will apply to all SBML model that were added after calling
set_sbml_global_options
. This means that usually you want to ensure that
SBML integration optin setting should be first thing you do in your
Python steppable file. If you want ot retrieve options simply type:
options = self.get_sbml_global_options()
notice that options can be None indicating that options have not been set (this is fine) and the default SBML integrator options will be applied.
Let us see how we associate SBML model with several cells using add_sbml_to_cell_ids
:
self.add_sbml_to_cell_ids(model_file=model_file, model_name='dp',
cell_ids=list(range(1, 11)), step_size=0.5, initial_conditions=initial_conditions)
This function looks relatively simple but it does quite a lot if you look under the hood. The first argument is path to SBML models file. The second one is model alias - it is a name you choose for model. It is an arbitrary model identifier that you use to retrieve model values. The third argument is a Python list that contains cell ids to which CC3D wil attach an instance of the SBML Solver.
Note
Each cell will get separate SBML solver object. SBML Solver objects associated with cells or free floating SBML Solvers are independent.
The fourth argument specifies the size of the integration step – here we
use value of 0.5
time unit. The fifth argument passes initial conditions
dictionary. Integration step size and initial conditions arguments are
optional.
Each SBMLSolver function that associates models with a cell or adds
free floating model calls libroadrunner
functions that parse SBML,
translate it to to very fast LLVM code. Everything happens
automatically and produces optimized solvers which are much faster than
solvers that rely on some kind of interpreters.
Next five function calls to self.add_free_floating_sbml
create instances of
SBML solvers which are not associated with cells but, as you can see,
have distinct names. This is required because when we want to refer to
such solver to extract model values we will do it using model name. The
reason all models attached to cells have same name was that when we
refer to such model we pass cell object and a name and this uniquely
identifies the model. Free floating models need to have distinct names
to be uniquely identified. Notice that last 3 calls to
self.add_free_floating_sbml
do not specify step size (we use default step
size 1.0 time unit) nor initial conditions (we use whatever defaults are
in the SBML code).
Finally, last two lines of start functions demonstrate how to add SBML Solver object to a single cell:
cell_20 = self.fetch_cell_by_id(20)
self.add_sbml_to_cell(model_file=model_file, model_name='dp', cell=cell_20)
Instead of passing list of cell ids we pass cell object (cell_20
).
We can also associate SBML model with certain cell types using the following syntax:
self.add_sbml_to_cell_types(model_file=model_file, model_name='dp', cell_types=[self.NONCONDENSING],
step_size=step_size, initial_conditions=initial_conditions)
This time instead of passing list of cell ids we pass list of cell types.
Let us move on to step function. First call we see there, is
self.timestep_sbml
. This function carries out integration of all SBML
Solver instances defined in the simulation. The integration step can be
different for different SBML Solver instances (as shown in our example).
To check the values of model species after integration step we can call e.g.
print('Free Floating Medium_dp2', self.sbml.Medium_dp2.values())
These functions check and print model variables for free floating model
called Medium_dp2
.
The next set of function calls:
if mcs == 3:
Medium_dp2 = self.sbml.Medium_dp2
Medium_dp2['S1'] = 10
Medium_dp2['S2'] = 0.5
set new state for for free floating model called Medium_dp2
. If we
wanted to print state of the model dp belonging to cell object called
cell20
we would use the following syntax:
print('cell_20, dp=', cell_20.sbml.dp.values())
To assign new values to dp model variables for cell20 we use the following syntax:
cell_20.sbml.dp['S1'] = 10
cell_20.sbml.dp['S2'] = 0.5
Note
We access free-floating SBML solved via self.sbml.MODEL_ALIAS
syntax whereas SBML solvers associated with a particular cell are accessed using reference to cell objects e.g. cell_20.sbml.MODEL_ALIAS
Another useful operation within SBML Solver capabilities is deletion of models. This comes handy when at certain point in your simulation you no longer need to solve ODE’s described in the SBML model. This is the syntax that deletes SBML from cell ids:
self.delete_sbml_from_cell_ids(model_name='dp', cell_ids=list(range(1, 11)))
As you probably suspect, we can delete SBML Solver instance from cell types:
self.delete_sbml_from_cell_types(model_name='dp' ,cell_types=range[self.A,self.B])
from single cell:
self.delete_sbml_from_cell(model_name='dp',cell=cell20)
or delete free floating SBML Solver object:
self.delete_free_floating_sbml(model_name='Medium_dp2'))
Note
When cells get deleted all SBML Solver models are deleted automatically. You do not need to call delete_sbml
functions in such a case.
Sometimes you may encounter a need to clone all SBML models from one cell to another (e.g. in the mitosis updateAttributes function where you clone SBML Solver objects from parent cell to a child cell). SBML Solver lets you do that very easily:
cell_10 = self.fetch_cell_by_id(10)
cell_25 = self.fetch_cell_by_id(25)
self.copy_sbml_simulators(from_cell=cell_10, to_cell=cell_25)
What happens here is that source cell (from_cell
) provides SBML Solver
object templates and based on these templates new SBML Solver objects
are gets created and CC3D assigns them to target cell (to_cell
). All
the state variables in the target SBML Solver objects are the same as
values in the source objects.
If you want to copy only select models you would use the following syntax:
cell_10 = self.fetch_cell_by_id(10)
cell_25 = self.fetch_cell_by_id(25)
self.copy_sbml_simulators(from_cell=cell_10, to_cell=cell_25, sbml_names=['dp'])
As you can see there is third argument - a Python list that specifies which models to copy. Here we are copying only dp models. All other models associated with parent cells will not be copied.
This example demonstrates most important capabilities of SBML Solver.
The next example shows slightly more complex simulation where we reset
initial condition of the SBML model before each integration step
(Demos/SBMLSolverExamples/DeltaNotch
).
Full description of the Delta-Notch simulation is in the introduction to CompuCell3D Manual. The Delta-Notch example demonstrates multi-cellular implementation of Delta-Notch mutual inhibitory coupling. In this juxtacrine signaling process, a cell’s level of membrane-bound Delta depends on its intracellular level of activated Notch, which in turn depends on the average level of membrane-bound Delta of its neighbors. In such a situation, the Delta-Notch dynamics of the cells in a tissue sheet will depend on the rate of cell rearrangement and the fluctuations it induces. While the example does not explore the richness due to the coupling of sub-cellular networks with inter-cellular networks and cell behaviors, it already shows how different such behaviors can be from those of their non-spatial simplifications. We begin with the Ordinary Differential Equation (ODE) Delta-Notch patterning model of Collier in which juxtacrine signaling controls the internal levels of the cells’ Delta and Notch proteins. The base model neglects the complexity of the interaction due to changing spatial relationships in a real tissue:
where and are the concentrations of activated Delta and Notch proteins inside a cell, is the average concentration of activated Delta protein at the surface of the cell’s neighbors, and are saturation constants, and are Hill coefficients, and is a constant that gives the relative lifetimes of Delta and Notch proteins.
Figure 18 Diagram of Delta-Notch feedback regulation between and within cells.
For the sake of simplicity let us assume that we downloaded SBML model implementing Delta-Notch ODE’s. How do we use such SBML model in CC3D? Here is the code:
from random import uniform
from cc3d.core.PySteppables import *
class DeltaNotchClass(SteppableBasePy):
def __init__(self, frequency=1):
SteppableBasePy.__init__(self, frequency)
def start(self):
# adding options that setup SBML solver integrator
# these are optional but useful when encounteting integration instabilities
options = {'relative': 1e-10, 'absolute': 1e-12}
self.set_sbml_global_options(options)
model_file = 'Simulation/DN_Collier.sbml'
self.add_sbml_to_cell_types(model_file=model_file, model_name='DN', cell_types=[self.TYPEA], step_size=0.2)
for cell in self.cell_list:
dn_model = cell.sbml.DN
dn_model['D'] = uniform(0.9, 1.0)
dn_model['N'] = uniform(0.9, 1.0)
cell.dict['D'] = dn_model['D']
cell.dict['N'] = dn_model['N']
def step(self, mcs):
for cell in self.cell_list:
delta_tot = 0.0
nn = 0
for neighbor, commonSurfaceArea in self.get_cell_neighbor_data_list(cell):
if neighbor:
nn += 1
delta_tot += neighbor.sbml.DN['D']
if nn > 0:
D_avg = delta_tot / nn
cell.sbml.DN['Davg'] = D_avg
cell.dict['D'] = D_avg
cell.dict['N'] = cell.sbml.DN['N']
self.timestep_sbml()
In the start function we add SBML model (Simulation/DN_Collier.sbml
) to
all cells of type A
(it is the only cell type in this simulation besides
Medium
). Later in the for loop we initialize D
and N
species from the
SBML using random values so that each cell has different SBML starting
state. We also store the initial SBML in cell dictionary for
visualization purposes – see full code in the
Demos/SBMLSolverExamples/DeltaNotch
. In the step function for each
cell we visit its neighbors and sum value of Delta in the neighboring
cells. We divide this value by the number of neighbors (this gives
average Delta concentration in the neighboring cells - D_avg
). We pass
Davg to the SBML Solver for each cell and then carry out integration for
the new time step. Before calling self.timestep_sbm
function we store
values of Delta and Notch concentration in the cell dictionary, but we
do it for the visualization purposes only. As you can see from this
example SBML Solver programing interface is convenient to use, not to
mention SBML Solver itself which is very powerful tool which allows
coupling cell-level and sub-cellular scales.