3.4. Modeling Schemes

E-Cell is a multi-algorithm simulator. It can run any kind of simulation algorithms, both discrete and continuous, and these simulation algorithms can be used in any combinations. This section exlains how you can find appropriate set of object classes for your modeling and simulation projects. This section does not give a complete list of available object classes nor detailed usage of those classes. Read the chapter "Standard Dynamic Module Library" for more info.

3.4.1. Discrete Or Continuous ?

E-Cell can model both discrete and continuous processes, and these can be mixed in simulation. The system models discrete and continuous systems by discriminating two different types of Process and Stepper objects: discrete Process / Stepper and continuous Process / Stepper.

NoteVariable is discrete and continous
 

Variable and System do not have special discrete and continuous classes. The base Variable class supports both discrete and continous operations, because it can be connected to any types of Process and Stepper objects. System objects do not do any computation that needs to discriminate discrete and continuos.

3.4.1.1. Discrete classes

A Process object that models discrete changes of one or more Variable objects is called a discrete Process, and it must be used in conjunction with a discrete Stepper. A discrete Process directly changes the values of related Variable objects when its Stepper requests to do so.

There are two types of discrete Process / Stepper classes: discrete and discrete event.

  • Discrete

    A discrete Process changes values of connected Variable objects (i.e. appear in its VariableReferenceList property) discretely. In the current version, there is no special class named DiscreteProcess, because the base Process class is already a discrete Process by default. The manner of the change of Variable values is determined from values of its accessor VariableReferences, its property values, and sometimes the current time of the Stepper. Unlike discrete event Process, which is explained in the next item, it does not necessary specify when the discrete changes of Variable values occur. Instead, it is unilaterally determined and fired by a discrete Stepper.

    A Stepper that requires all Process objects connected is discrete Process objects is call a discrete Stepper. The current version has no special class DiscreteStepper, because the base Stepper class is already discrete.

  • Discrete event

    Discrete event is a special case of discreteness. The system provides DiscreteEventStepper and DiscreteEventProcess classes for discrete-event modeling. In addition to the ordinary firing method (fire() method) of the base Process class, the DiscreteEventProcess defines a method to calculate when is the next occurrence of the event (the discrete change of Variable values that this discrete event Process models) from values of its accessor VariableReferences, its property values, and the current time of the Stepper. DiscreteEventStepper uses information given by this method to determine when each of discrete event Process should be fired. DiscreteEventStepper is instantiatable. See the chapter Standard Dynamic Module Library for more detailed description of how DiscreteEventStepper works.

3.4.1.2. Continuous classes

On the other hand, a Process that calculates continuous changes of Variable objects is called a continuous Process, and is used in combination with a continuous Stepper. Continuous Process objects simulate the phenomena that represents by setting velocities of connected Variable objects, rather than directly changing their values in the case of discrete Process objects. A continuous Stepper integrates the values of Variable objects from the velocities given by the continuous Process objects, and determines when the velocities should be recalculated by the Process objects. A typical application of continuous Process and Stepper objects is to implement differential equations and differential equation solvers, respectively, to form a simulation system of the system of differential equations.

3.4.2. Some Available Discrete Classes

Followings are some available discrete classes.

3.4.2.1. NRStepper and GillespieProcess (Gillespie-Gibson pair)

An example of discrete-event simulation method provided by E-Cell is a variant of Gillespie's stochastic algorithm, the Next Reaction Method, or Gillespie-Gibson algorithm. NRStepper class implements this algorithm. When this Stepper is used in conjunction with GillespieProcess objects, which is a subclass of DiscreteEventProcess and calculates a time of the next occurence of the reaction using Gillespie's reaction probability equation and a random number, E-Cell conducts a Gillespie-Gibson stochastic simulation of elementary chemical reactions. In fact, the Next Reaction Method is nothing but a standard discrete event simulation algorithm, and NRStepper is just an alias of the DiscreteEventStepper class.

Usage of this pair of classes of objects is simple: just set the StepperID, VariableReferenceList and the rate constant property k of those GillespieProcess objects.

3.4.2.2. DiscreteTimeStepper

A type of discrete Stepper that is provided by the system is DiscreteTimeStepper. This class of Stepper, when instantiated, calls all discrete Process objects with a fixed user-specified time-interval. For example, if the model has a DiscreteTimeStepper with 0.001 (second) of StepInterval property, it fires all of its Process objects every milli-second. DiscreteTimeStepper is discrete time because it does not have time between steps; it ignores a signal from other Stepper objects (Stepper interruption) that notifies a change of system state (values of Variable objects) that may affect its Process objects. Such a change is reflected in the next step.

3.4.2.3. PassiveStepper

Another class of discrete Stepper is PassiveStepper. This can partially be seen as a DiscreteTimeStepper with an infinite StepInterval, but there is a difference. Unlike DiscreteTimeStepper, this does not ignore Stepper interruptions, which notify change in the system state that may affect this Stepper's Process objects.

This Stepper is used when some special procedures (coded in discrete Process objects) must be invoked when other Stepper object may have changed a value or a velocity of at least one Variable that this Stepper's Process objects accesses.

3.4.2.4. PythonProcess

PythonProcess allows users to script a Process object in full Python syntax.

initialize() and fire() methods can be scripted with InitializeMethod and FireMethod properties, respectively.

PythonProcess can be either discrete or continuous. This 'operation mode' can be specified by setting IsContinuous property. The default is false (0), or discrete. To switch to the continuous mode, set 1 to the property:

Process PythonProcess( PY1 )
{
    IsContinuous 1;
}

In addition to regular Python constructs, the following objects, methods, and attributes are available in both of the method properties (InitializeMethod and FireMethod):

  • Properties

    PythonProcess accepts arbitrary names of properties. For example, the following code creates two new properties.

    Process PythonProcess( PY1 )
    {
        NewProperty "new property";
        KK          3.0;
    }
    These properties can be use in Python methods:
    Process PythonProcess( PY1 )
    {
        # ... NewProperty and KK are set
    
        InitializeMethod "print NewProperty";
    
        FireMethod '''
    KK += 1.0
    print KK 
    ''';
    }
    A new property can also be created within Python methods.
    
    InitializeMethod "A = 3.0"; # A is created
        FireMethod "print A * 2";   # A can be used here
    
    These properties are treated as a global variable.

  • Objects

    • self

      This is the Process object itself. This has the following attributes:

      • Activity

        The current value of Activity property of this Process.

      • addValue( value )

        Add each VariableReference the value multiplied by the coefficient of the VariableReference.

        Using this method implies that this Process is discrete. Check that IsContinuous property is false.

      • getSuperSystem()

        This method gets the super system of this Process. See below for the attributes of System objects.

      • Priority

        The Priority property of this Process.

      • setFlux( value )

        Add each VariableReference's velocity the value multiplied by the coefficient of the VariableReference.

        Using this method implies that this Process is continuous. Check that IsContinuous property is true.

      • StepperID

        StepperID of this Process.

    • VariableReference

      VariableReference instances given in the VariableReferenceList property of this Process can be used in the Python methods. Each instance has the following attributes:

      • addFlux( value )

        Multiply the value by the Coefficient of this VariableReference, and add that to the Variable's velocity.

      • addValue( value )

        Add the value to the Value property of the Variable.

      • addVelocity( value )

        Add the value to the Velocity property of the Variable.

      • Coefficient

        The coefficient of the VariableReference

      • getSuperSystem()

        Get the super system of the Variable. A System object is returned.

      • MolarConc

        The concentration of the Variable in Molar [M].

      • Name

        The name of the VariableReference.

      • NumberConc

        The concentration in number [ num / size of the Variable's super system. ]

      • IsFixed

        Zero if the Fixed property of the Variable is false. Otherwise a non-zero integer.

      • IsAccessor

        Zero if the IsAccessor flag of the VariableReference is false. Otherwise a non-zero integer.

      • TotalVelocity

        The total current velocity. Usually of no use.

      • Value

        The value of the Variable

      • Velocity

        The provisional velocity given by the currently stepping Stepper. Usually of no use.

    • System

      A System object has the following attributes.

      • getSuperSystem()

        Get the super system of the System. A System object is returned.

      • Size

        The size of the System.

      • SizeN_A

        Equivalent to Size * N_A, where N_A is a Avogadro's number.

      • StepperID

        The StepperID of the System.

Here is an example uses of PythonProcess.

Example 3-2. An assignment rule by using PythonProcess


Process PythonProcess( PY1 )
{
    # IsContinuous 0; -- default
    FireMethod "S1.Value = S2.Value + S3.Value";
    VariableReferenceList [(S1)] [(S2)] [(S3)];
}

3.4.2.5. PythonEventProcess

This class enables users Python scripting of time-events. In addition to initialize() and fire(), updateStepInterval() method can be scripted with this class. Use UpdateStepIntervalMethod property to set this.

In addition to those of PythonProcess, the self object of PythonEventProcess has some more attributes:

  • StepInterval

    The most recent StepInterval calculated by the updateStepInterval() method.

  • DependentProcessList

    This attribute holds a tuple of IDs of dependent Processes of this Process.

This class of objects must be used with a DiscreteEventStepper.

This class is under development.

3.4.2.6. Other discrete classes

Stepper classes for explicit and implicit tau leaping algorithms are under development.

A flux-distribution method for hybrid dynamic/static simulation of biochemical pathways is available with the following classes: FluxDistributionStepper, FluxDistributionProcess, QuasiDynamicFluxProcess. Usage of this scheme is to be described.

3.4.3. Some Available Continuous Classes

E-Cell supports both Ordinary Differential Equation (ODE) and Differential-Algebraic Equation (DAE) models, and has Stepper classes for each type of formalisms.

Also, the system is shipped with some continuous Process classes. For example, MassActionFluxProcess calculates a reaction rate according to the law of mass action. ExpressionFluxProcess allows users to describe arbitrary rate equations in model files. PythonProcess and PythonFluxProcess are used to script Process objects in Python. Some enzyme kinetics rate laws are also available.

3.4.3.1. Generic ordinary differential Steppers

If your model is a system of ODEs, then in this version of the software (version 3.1.100) the recommended choice is ODE45Stepper. A current drawback of this Stepper is that it is not good at solving stiff systems. In near future, a high-performance replacement of this, ODEStepper, will be included.

Some other available ODE Stepper classes are ODE23Stepper, which employes a lower (the second) order integration algorithm, and FixedODE1Stepper that implements the simplest Euler algorithm without an adaptive step sizing mechanism.

These ODE Stepper classes except for the FixedODE1Stepper have some common property slots for user-specifiable parameters. Here is a partial list:

  • Tolerance

    An error tolerance in local truncation error. Giving this smaller numbers forces the Stepper to take smaller step sizes, and slows down the simulation. Greater numbers results in faster run with sacrifice of accuracy. A typical number is 1e-6.

  • MinStepInterval / MaxStepInterval

    Minimum and maximum limits in step size. These limits preceeds the Tolerance property above.

    These properties can also be useful to completely disable the adaptive step size control mechanism: set the same number to both of the property slots.

3.4.3.2. MassActionFluxProcess

MassActionFluxProcess is a class of Process for simple mass-actions. This class calculates a flux rate according to the irreversible mass-action. Use a property k to specify a rate constant.

3.4.3.3. ExpressionFluxProcess

ExpressionFluxProcess is designed for easy and efficient representations of continuous flux rate equations.

Expression property of this class accepts a plain text rate expression. The expression must be evaluated to give a flux rate in [ number / second ]. (Note that this is a number per second, not concentration per second.) Here is an example use of ExpressionFluxProcess:


Process ExpressionFluxProcess( P1 )
{
    k 0.1;
    Expression "k * S.Value";

    VariableReferenceList [ S :.:S -1 ] [ P :.:P 1 ];
}

Compared to PythonProcess or PythonFluxProcess below, it runs significantly faster with sacrifice of some flexibility in scripting.

The following shows elements those can be used in the Expression property. The set of available arithmetic operators and mathematical functions are meant to be equivalent to SBML level 2, except control structures.

  • Constants

    Numbers (e.g. 10, 10.33, 1.33e-5), true, false (equivalent to zero), pi (Pi), NaN (Not-a-Number), INF (Infinity), N_A (Avogadro's number), exp (the base of natural logarithms).

  • Arithmetic operators

    +, -, *, /, ^ (power; this can equivalently be written as pow( x, y )).

  • Built-in functions

    abs, ceil, exp, *fact, floor, log, log10, pow sqrt, *sec, sin, cos, tan, sinh, cosh, tanh, coth, *csch, *sech, *asin, *acos, *atan, *asec, *acsc, *acot, *asinh, *acosh, *atanh, *asech, *acsch, *acoth. (Functions with astarisk '*' are currently not available on the Windows version.)

    All functions but pow are unary functions. pow is a binary function.

  • Properties

    Similar to PythonProcess, ExpressionFluxProcess accepts arbitrary name properties in the model. Unlike PythonProcess, however, these properties of this class can hold only Real values.

  • Objects

    • self

      This Process object itself. This has the following attribute which is a sub set of that of PythonProcess:

      • getSuperSystem()

    • VariableReference

      VariableReference instances given in the VariableReferenceList property of this Process can be used in the expression. Each instance has the following set of attributes, which is a sub set of that of PythonProcess:

      • Value

      • MolarConc

      • NumberConc

      • TotalVelocity

      • Velocity

    • System

      A System object has the following two attributes.

      • Size

      • SizeN_A

Below is an example of the basic Michaelis-Menten reaction programmed with the ExpressionFluxProcess.

Example 3-3. Michaelis-Menten reaction with ExpressionFluxProcess


Process ExpressionFluxProcess( P )
{
    Km    1.0;
    Kcat  10;

    Expression "E.Value * Kcat * S.MolarConc / ( S.MolarConc + Km )";

    VariableReferenceList [ S :.:S -1 ] [ P :.:P 1 ] [ E :.:E 0 ];
}

3.4.3.4. Some pre-defined reaction rate classes

See the standard dynamic module library reference for availability of some enzyme kinetics Process classes.

3.4.3.5. PythonFluxProcess

PythonFluxProcess is almost the same as PythonProcess, except that (1) it takes just a Python expression (instead of statements) to its Expression property, and (2) similar to ExpressionFluxProcess, the evaluated value of the expression is implicitly passed to the setFlux() method.

3.4.3.6. Generic differential-algebraic Steppers

For DAE models, use DAEStepper. The model must form a valid index-1 DAE system. When a DAE Stepper detects one or more discrete Process objects, it assumes that these are algebraic Process objects. Thus, all discrete Process objects in a DAE Stepper must be algebraic. See below for what is algebraic Process.

NoteDAE Stepper objects can be used for ODE systems
 

Because it can be viewed that ODE is a special case of DAE problems which does not have a algebraic equations, but only differential equations, a DAE Stepper can be used to run an ODE model. However, ODE Steppers are specialized for ODE problems, in terms of both the selection of integration algorithms and implementation issues, and generally use of an ODE Stepper benefits in performance when the model is a system of ODEs.

Those properties of ODE Stepper classes described above (such as the Tolerance property) are also available for DAE Stepper classes.

3.4.3.7. Algebraic Processes

This is a type of discrete Process, but placed here because it is used with a DAE Stepper, which is continuous.

In principle, continuous Process objects must be connected with continuous Stepper instances, and a discrete Stepper is assumed to take only discrete Process objects. However, there are some exceptions. One of such is the algebraic processes. Strangely enough, in DAE simulations, seemingly discrete algebraic equations are solved continuously in conjunction with other differential equations.

Algebraic equations in E-Cell has the following form:

0 = g( t, x )
where t is the time and x is a vector of variable references.

The DAE solver system of E-Cell uses Activity property of Process objects to represent the value of the algebraic function g( t, x ). An algebraic Process must not change values of Variable objects explicitly. The DAE Stepper does this job of finding a point where the equation g() holds.

When modeling, be careful about coefficients of VariableReferences of an algebraic Process. In most cases, simply set unities. The solver respects these numbers when solving the system. For example, if the coefficient of A is zero, it does not change the variable when trying to find the solution, while it is used in the calculation of the equation.

As a means of describing algebraic equations, ExpressionAlgebraicProcess is available. The usage is the same as ExpressionFluxProcess, except that the evaluation of its expression is interpreted as the value of the algebraic function g().

The following examble describes an equation


a * A + B = 10,  a = 1.5

Example 3-4. A simple algebraic equation using ExpressionAlgebraicProcess


Stepper DAEStepper( DAE1 ) {}

Process ExpressionFluxProcess( P )
{
    StepperID DAE1;

    a    1.5;

    Expression "( a * A + B ) - 10";

    VariableReferenceList [ A :.:A 1 ] [ B :.:B 1 ];
}

To use C++ or PythonProcess for algebraic equations, call setActivity() method to set the value of the equation. The following is an example with a PythonProcess:

Example 3-5. A simple algebraic equation using PythonProcess


Process PythonProcess( PY )
{
    a    1.5;

    FireMethod "self.setActivity( ( a * A + B ) - 10 )";

    VariableReferenceList [ A :.:A 1 ] [ B :.:B 1 ];
}

3.4.3.8. Power-law canonical DEs (S-System and GMA)

ESSYNSStepper supports S-System and GMA simulations by using the ESSYNS algorithm. A ESSYNSStepper must be connected with either a SSystemProcess or a GMAProcess. Use SSystemMatrix or GMAMatrix property to set the system parameters.

A sample model under the directory doc/sample/ssystem/ gives an example usage.

These modules are still under development. More descriptions to come....