E-Cell Simulation Environment Version 3.1.100 User's Manual (Draft: Dec. 18, 2003)
PrevChapter 4. Modeling TutorialNext

4.5. A Simple Deterministic / Stochastic Composite Simulation

E-Cell 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.

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

4.5.2. 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?).

Example 4-4. composite.em


@{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;
    }
    
}
PrevHomeNext
Making the Model Switchable Between AlgorithmsUpCustom equations