4.3. Using Deterministic Differential Equations

The previous section described how to create a model that runs with the stochastic Gillespie's algorithm. E-Cell 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.

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

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

Example 4-2. A simple deterministic mass-action model


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