4.2. Using Gillespie algorithm

E-Cell Simulation Environment comes with a set of classes for simulations using Gillespie's stochastic algorithm.

4.2.1. 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.

4.2.2. 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.

4.2.3. 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; }

    # ...
}

4.2.4. 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;
    }
        
    # ...
}

4.2.5. 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;
    }
}

4.2.6. 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.

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;
    }
}