E-Cell Simulation Environment Version 3.1.100 User's Manual (Draft: Dec. 18, 2003) | ||
---|---|---|
Prev | Chapter 4. Modeling Tutorial | Next |
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.
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 2Rate 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.
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; } }