E-Cell Simulation Environment Version 3.1.100 User's Manual (Draft: Dec. 18, 2003) | ||
---|---|---|
Prev | Chapter 4. Modeling Tutorial | Next |
E-Cell Simulation Environment comes with a set of classes for simulations using Gillespie's stochastic algorithm.
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.
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.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; } # ... }
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; } # ... }
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; } }
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.
Example 4-1. The simplest Gillespie-Gibson model
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; } }