Modeling Tutorial¶
This chapter is a simple modeling tutorial using ECELL.
Running the model¶
All the examples in this section can be run by the following procedure.
Create and save the model file (for example,
simple-gillespie.em
) with a text editor.Convert the EM file to an EML file by
ecell3-em2eml
command.Run it in GUI mode with
gecell
command.or, in the script mode with
ecell3-session
command (see the following chapter):
Using Gillespie algorithm¶
APP comes with a set of classes for simulations using Gillespie’s stochastic algorithm.
A Trivial Reaction System¶
At the very first, let us start with the simplest possible stable system
of elementary reactions, which has two variables (in this case the
numbers of molecules of molecular species) and a couple of elementary
reaction processes. Because elementary reactions are irreversible, at
least two instances of the reactions are needed for the system to be
stable. The reaction system looks like this: – P1 –> S1 S2 <– P2 –
S1
and S2
are molecular species, and P1
and P2
are
reaction processes. Rate constants of both reactions are the same: 1.0
[1/s]. Initial numbers of molecules of S1
and S2
are 1000 and 0,
respectively. Because rate constants are the same, the system has a
steady state at S1 == S2 == 500
.
Specifying the Next Reaction method¶
NRStepper class implements Gibson’s efficient variation of the Gillespie algorithm, or the Next Reaction (NR) method.
To use the NRStepper in your simulation model, write like this in your EM file:
Stepper NRStepper( NR1 )
{
# no property
}
In this example, the NRStepper has the StepperID ‘NR1
‘. For now,
no property needs to be specified for this object.
Defining the compartment¶
Next, define a compartment, and connect it to the STEPPER NR1
.
Because this model has only a single compartment, we use the root sytem
(/
). Although this model does not depend on size of the compartment
because all reactions are first order, it is a good idea to always
define the size explicitly rather than letting it defaults to 1.0
.
Here we set it to 1e-15
[L].
System System( / )
{
StepperID NR1;
Variable Variable( SIZE ) { Value 1e-15; }
# ...
}
Defining the variables¶
Now define the main variables S1
and S2
. Use ‘Value’ properties
of the objects to set the initial values.
System System( / )
{
# ...
Variable Variable( S1 )
{
Value 1000;
}
Variable Variable( S2 )
{
Value 0;
}
# ...
}
Defining reaction processes¶
Lastly, create reaction process instances P1
and P2
.
GillespieProcess class works together with the NRStepper to simulate
elementary reactions.
Two different types of properties, k and VariableReferenceList, must be
set for each of the GillespieProcess object. k is the rate constant
parameter in [1/sec] if the reaction is first-order, or [1/sec/M] if it
is a second-order reaction. (Don’t forget to define SIZE
VARIABLE if
there is a second-order reaction.) Set VariableReferenceList properties
so that P1
consumes S1
and produce S2
, and P2
uses
S2
to create S1
.
System System( / )
{
# ...
Process GillespieProcess( P1 ) # the reaction S1 --> S2
{
VariableReferenceList [ S :.:S1 -1 ]
[ P :.:S2 1 ];
k 1.0; # the rate constant
}
Process GillespieProcess( P2 ) # the reaction S2 --> S1
{
VariableReferenceList [ S :.:S2 -1 ]
[ P :.:S1 1 ];
k 1.0;
}
}
Putting them together¶
Here is the complete EM of the model that really works. Run this model
with gecell
and open a TracerWindow to plot trajectories of S1
and S2
. You will see those two VARIABLEs immediately reaching the
steady state around 500.0. If you zoom around the trajectories, you will
be able to see stochastic fluctuations.
Stepper NRStepper( NR1 )
{
# no property
}
System System( / )
{
StepperID NR1;
Variable Variable( SIZE ) { Value 1e-15; }
Variable Variable( S1 )
{
Value 1000;
}
Variable Variable( S2 )
{
Value 0;
}
Process GillespieProcess( P1 ) # the reaction S1 --> S2
{
VariableReferenceList [ S :.:S1 -1 ]
[ P :.:S2 1 ];
k 1.0; # the rate constant
}
Process GillespieProcess( P2 ) # the reaction S2 --> S1
{
VariableReferenceList [ S :.:S2 -1 ]
[ P :.:S1 1 ];
k 1.0;
}
}
Using Deterministic Differential Equations¶
The previous section described how to create a model that runs with the stochastic Gillespie’s algorithm. ECELL is a multi-algorithm simulator, and different algorithms can be used to run the model. This section explains a way to use a deterministic differential equation solver to run the system of simple mass-action reactions.
Choosing Stepper and Process classes¶
In the current version, we recommend ODE45Stepper class as a general-purpose STEPPER for differential equation systems. This STEPPER implements an explicit fourth order numerical integration algorithm with a fifth-order error control.
MassActionFluxProcess is the continuous differential equation conterpart
of the discrete-event GillespieProcess. Unlike GillespieProcess,
MassActionFluxProcess does not have limitation in the order of the
reaction mechanism. For example, it can handle the reaction like this:
S0 + S1 + 2 S2 --> P0 + P1
.
Converting the model¶
Converting the trivial reaction system model for Gillespie to use differential equations is very easy; just replace NRStepper with ODE45Stepper, and change the classname of GillespieProcess to MassActionFluxProcess.
The following is the model of the trivial model that runs on the differential ODE45Stepper. You will get similar simulation result than the stochastic model in the previous section. However, if you zoom, you would notice that the stochastic fluctuation can no longer be observed because the model is turned from stochastic to deterministic.
Stepper ODE45Stepper( ODE1 )
{
# no property
}
System System( / )
{
StepperID ODE1;
Variable Variable( SIZE ) { Value 1e-15; }
Variable Variable( S1 )
{
Value 1000;
}
Variable Variable( S2 )
{
Value 0;
}
Process MassActionFluxProcess( P1 )
{
VariableReferenceList [ S0 :.:S1 -1 ]
[ P0 :.:S2 1 ];
k 1.0;
}
Process MassActionFluxProcess( P2 )
{
VariableReferenceList [ S0 :.:S2 -1 ]
[ P0 :.:S1 1 ];
k 1.0;
}
}
Making the Model Switchable Between Algorithms¶
Fortunately, at least as far as the model has only elementary reactions, switching between these stochastic and deterministic algorithms is just to switch between NRStepper/GillespieProcess pair and ODE45Stepper/MassActionFluxProcess pair of classnames. Both PROCESS classes takes the same property ‘k’ with the same unit.
Some use of EMPY macros makes the model generic. In the following
example, setting the PYTHON variable TYPE
to ODE
makes it run in
deterministic differential equation mode, and setting TYPE
to NR
turns it stochastic.
@{ALGORITHM='ODE'}
@{
if ALGORITHM == 'ODE':
STEPPER='ODE45Stepper'
PROCESS='MassActionFluxProcess'
elif ALGORITHM == 'NR':
STEPPER='NRStepper'
PROCESS='GillespieProcess'
else:
raise 'unknown algorithm: %s' % ALGORITHM
}
Stepper @(STEPPER)( STEPPER1 )
{
# no property
}
System System( / )
{
StepperID STEPPER1;
Variable Variable( SIZE ) { Value 1e-15; }
Variable Variable( S1 )
{
Value 1000;
}
Variable Variable( S2 )
{
Value 0;
}
Process @(PROCESS)( P1 )
{
VariableReferenceList [ S0 :.:S1 -1 ]
[ P0 :.:S2 1 ];
k 1.0;
}
Process @(PROCESS)( P2 )
{
VariableReferenceList [ S0 :.:S2 -1 ]
[ P0 :.:S1 1 ];
k 1.0;
}
}
A Simple Deterministic / Stochastic Composite Simulation¶
ECELL can drive a model with multiple STEPPER objects. Each STEPPER can implement different simulation algorithms, and have different time scales. This framework of multi-algorithm, multi-timescale simulation is quite generic, and virtually any combination of any number of different types of sub-models is possible. This section exemplifies a tiny model of coupled ODE and Gillespie reactions.
A tiny multi-timescale reactions model¶
Consider this tiny model of four VARIABLEs and six reaction PROCESSes:
– P1 –> – P4 –> S1 S2 – P3 –> S3 S4 ^ <– P2 – <– P5 – | | |
\ _______________ P6
____________________/ Although it may look
complicated at first glance, this system consists of two instances of
the ‘trivial’ system we have modeled in the previous sections coupled
together: Sub-model 1: – P1 –> S1 S2 <– P2 – and Sub-model 2: – P4
–> S3 S4 <– P5 – These two sub-models are in turn coupled by reaction
processes P3
and P6
. Because time scales of P3
and P6
are determined by S2
and S4
, respectively, P3
belongs to the
sub-model 1, and P6
is a part of the sub-model 2. Sub-model 1: S2 –
P3 –> S3 S1 <– P6 –> S4 :Sub-model 2 Rate constants of the main
reactions, P1
, P2
, P4
, and P5
are the same as the
previous model: 1.0
[1/sec]. But the ‘bridging’ reactions are slower
than the main reactions: 1e-1
for P3
and 1e-3
for P6
.
Consequently, sub-models 1 and 2 would have approximately
1e-1 / 1e-3 == 1e-2
times different steady-state levels. Because the
rate constants of the main reactions are the same, this implies time
scale of both sub-models are different.
Writing model file¶
The following code implements the multi-time scale model. The first two
lines specify algorithms to use for those two parts of the model.
ALGORITHM1
variable specifies the algorithm to use for the sub-model
1, and ALGORITHM2
is for the sub-model 2. Values of these variables
can either be 'NR'
or 'ODE'
.
For example, to try pure-stochastic simulation, set these variables like this:
@{ALGORITHM1='NR'}
@{ALGORITHM2='NR'}
Setting ALGORITHM1
to 'NR'
and ALGORITHM2
to 'ODE
would
be an ideal configuration. This runs a magnitude faster than the
pure-stochastic configuration.
@{ALGORITHM1='NR'}
@{ALGORITHM2='ODE'}
Also try pure-deterministic run.
@{ALGORITHM1='ODE'}
@{ALGORITHM2='ODE'}
In this particular model, this configuration runs very fast because the system easily reaches the steady-state and stiffness of the model is low. However, this does not necessary mean pure-ODE is always the fastest. Under some situations NR/ODE composite simulation exceeds both pure-stochastic and pure-deterministic (reference?).
@{ALGORITHM1= ['NR' or 'ODE']}
@{ALGORITHM2= ['NR' or 'ODE']}
# a function to give appropriate class names.
@{
def getClassNamesByAlgorithm( anAlgorithm ):
if anAlgorithm == 'ODE':
return 'ODE45Stepper', 'MassActionFluxProcess'
elif anAlgorithm == 'NR':
return 'NRStepper', 'GillespieProcess'
else:
raise 'unknown algorithm: %s' % ALGORITHM1
}
# get classnames
@{
STEPPER1, PROCESS1 = getClassNamesByAlgorithm( ALGORITHM1 )
STEPPER2, PROCESS2 = getClassNamesByAlgorithm( ALGORITHM2 )
}
# create appropriate steppers.
# stepper ids are the same as the ALGORITHM.
@('Stepper %s ( %s ) {}' % ( STEPPER1, ALGORITHM1 ))
# if we have two different algorithms, one more stepper is needed.
@(ALGORITHM1 != ALGORITHM2 ? 'Stepper %s( %s ) {}' % ( STEPPER2, ALGORITHM2 ))
System CompartmentSystem( / )
{
StepperID @(ALGORITHM1);
Variable Variable( SIZE ) { Value 1e-15; }
Variable Variable( S1 )
{
Value 1000;
}
Variable Variable( S2 )
{
Value 0;
}
Variable Variable( S3 )
{
Value 1000000;
}
Variable Variable( S4 )
{
Value 0;
}
Process @(PROCESS1)( P1 )
{
VariableReferenceList [ S0 :.:S1 -1 ] [ P0 :.:S2 1 ];
k 1.0;
}
Process @(PROCESS1)( P2 )
{
VariableReferenceList [ S0 :.:S2 -1 ] [ P0 :.:S1 1 ];
k 1.0;
}
Process @(PROCESS1)( P3 )
{
VariableReferenceList [ S0 :.:S2 -1 ] [ P0 :.:S3 1 ];
k 1e-1;
}
Process @(PROCESS2)( P4 )
{
StepperID @(ALGORITHM2);
VariableReferenceList [ S0 :.:S3 -1 ] [ P0 :.:S4 1 ];
k 1.0;
}
Process @(PROCESS2)( P5 )
{
StepperID @(ALGORITHM2);
VariableReferenceList [ S0 :.:S4 -1 ] [ P0 :.:S3 1 ];
k 1.0;
}
Process @(PROCESS2)( P6 )
{
StepperID @(ALGORITHM2);
VariableReferenceList [ S0 :.:S4 -1 ] [ P0 :.:S1 1 ];
k 1e-4;
}
}
Custom equations¶
Complex flux rate equations¶
The simplest way to script custom rate equations is to use
ExpressionFluxProcess. Here is an example taken from the Drosophila
sample model which you can find under
${datadir}/doc/ecell/samples/Drosophila
[1] In this expression,
Size * N_A of the supersystem of the PROCESS is used to keep the unit
of the expression [ num / second ].
Process ExpressionFluxProcess( R_toy1 )
{
vs 0.76;
KI 1;
Expression "(vs*KI) / (KI + C0.MolarConc ^ 3)
* self.getSuperSystem().SizeN_A";
VariableReferenceList [ P0 :.:M 1 ] [ C0 :.:Pn 0 ];
}
FIXME: some more examples
Algebraic equations¶
Use of ExpressionAlgebraicProcess is the easiest method to describe algebraic equations.
Be careful about the coefficients of the VARIABLEREFERENCEs. (Usually just set unities.)
FIXME: some more examples here