Stochastic Simulation in Virtual Cell

 

1. Introduction to stochastic approaches

2. A quick tour of Virtual Cell stochastic simulation

3. References

1.   Introduction to stochastic approaches

Biological systems are highly complex. Different numerical methods are required to simulate different biological systems. Deterministic simulation is valid as long as the concentrations of simulated populations are high. However, the concentrations of these factors could be very low in living cells. In these latter cases, the biology is more accurately simulated by discrete approaches. As a result, stochastic simulation has been implemented in Virtual Cell to allow users to describe the discrete nature of changes in cell systems.

 

Exact stochastic simulation was originally introduced by Gillespie [1]. However, Gillespie’s algorithms require a substantial amount of computational effort to simulate a complex system. To optimize Gillespie’s methods, Gibson and Bruck proposed the Next Reaction method [2], which uses only a single random number per simulation event and takes time proportional to the logarithm of the number of reactions. They achieved better performance also by utilizing a dependency graph and an indexed priority queue. Virtual Cell has implemented the Gibson Next Reaction method to exactly solve biological systems in a stochastic way (Version 4.4).

 

Stochastic algorithms have to deal with the same "stiffness" problem which is encountered in deterministic approaches. In order to capture the fast dynamics of the system, the entire simulation is slowed down significantly. The underlying idea is to partition the system into fast and slow sub-systems and speed up the simulation by estimating the fast system and exactly solving the slow systems.  Howard introduced a set of hybrid methods [3] [4]. It approximates the fast reactions as a continuous Markov process, using a chemical Langevin equation, and accurately describes the slow dynamics using the Gibson algorithm. Different integrators are used for approximate numerical solution of chemical Langevin equations. Three hybrid solvers from Howard will be implemented in Virtual Cell (Version 4.5). They are fixed time step Gibson-Euler-Maruyama, fixed time step Gibson-Milstein and adaptive time step Gibson-Milstein.

 

Below is an example that shows the difference between deterministic, stochastic and hybrid (using Gibson-Milstein in this example) simulations by using HMM system as an example (Figure 1, Figure 2 and Figure 3).

 

 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


               

                Figure 1. VCell Deterministic simulation for Michaelis-Menten model

 

 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


                  Figure 2. VCell Stochastic simulation for Michaelis-Menten model

 

 
 
 

 

 

 

 

 

 

 

 

 


                  Figure 3. VCell Hybrid Simulation for Michaelis-Menten model

 

2.  A quick tour of Virtual Cell stochastic simulation

In the following section, we will discuss how to create and run the stochastic simulations and view the results in Virtual Cell. We shall use Michaelis-Menten as an example.

 

·                    To create a stochastic application.

 

To stochastically simulate a system, users have to create a stochastic application based on the physiology. There are four ways to create such an application.

 

a)      Select the “Application” menu and click on the “New” item and choose “Stochastic Application” sub item (Figure 4). This will set up a new stochastic application.

 

b)      Select the model name in the tree and right click and choose “Create Stochastic Application” (Figure 5). This will set up a new stochastic application.

 

c)      Click on the existing application name and right click and choose “Copy As” item and further down choose “Stochastic Application” (Figure 6). This will copy the original application (whatever deterministic or stochastic) to a stochastic application while keeping the original parameters including structure information and initial condition etc.

 

d)      Click on the existing application name and right click and choose “New” item and further down choose “Stochastic Application” (Figure 7). This will set up a new stochastic application.

 
 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


                 Figure 4. Create a stochastic application from “Application” Menu

 

  
 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


                   Figure 5. Create a stochastic application from model tree root

 

 

 
 
 

 

 

 

 

 

 

 

 

 

 

 

 


 

 

       Figure 6. Copy a stochastic application by existing one from “Copy As” menu item

 
 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


                 Figure 7. Create a stochastic application from “New” menu item

 

 

·                    To set up structure size and initial conditions

 

Structure size and initial conditions are required by stochastic simulation. These are created in the same way that structure size and initial conditions are created in deterministic models. Users can follow instructions in the Virtual Cell user guide to set up the structure size and initial conditions. VCell user guide is located at http://www.vcell.org/login/documentation.html.

·                    To set up simulation parameters

 

After creating a simulation (see Virtual Cell user guide for instructions), users are responsible to choose preferred solver and set up the simulation parameters. There will be four stochastic solvers in Virtual Cell. One is an exact solver (available in Version 4.4) and the other three are hybrid solvers (available in Version 4.5). Each of them requires different parameters to properly run the simulation. We list the different parameter sets below according to different solvers. To set up simulation parameters, users just need to click on “Edit” button in the simulation panel.

 

a)      Gibson-Bruck Exact Solver (Figure 8)

 

STARTING TIME: the time when simulation starts.

ENDING TIME: the time when simulation ends.

RANDOM SEED: a random number generated by PC time, which is used to produce a series of uniformly distributed random numbers.

CUSTOMIZED SEED: a user specified number, which is used to produce a series of uniformly distributed random numbers.

NUMBER OF TRIALS: the number of multiple trials to be run.

OUTPUT OPTIONS: use keep every number of samples or output time interval.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


                     Figure 8. Parameters required by exact stochastic solver

b)      Hybrid (Gibson-Euler-Maruyama) and Hybrid (Gibson-Milstein) (Figure 9)

 

STARTING TIME: the time when simulation starts.

ENDING TIME: the time when simulation ends.

SDE TIME STEP: the maximum time step of the SDE numerical integrator.

RANDOM SEED: a random number generated by PC time, which is used to produce a series of uniformly distributed random numbers.

CUSTOMIZED SEED: a user specified number, which is used to produce a series of uniformly distributed random numbers.

NUMBER OF TRIALS: the number of multiple trials to be run.

EPSILON: minimum number of molecules, both reactant and product species, required for approximation as a continuous Markov process.

LAMBDA: minimum rate of reaction required for approximation to a continuous Markov process.

MSR TOLERANCE: maximum allowed effect of executing multiple slow reactions per numerical integration of the SDEs.

OUTPUT OPTIONS: use output time interval.

 

For more information, please visit http://hysss.sourceforge.net/publications.shtml.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


.

 

                  Figure 9. Parameters required by fixed time step hybrid solvers

 

c)      Hybrid (Adaptive time step Gibson-Milstein) (Figure 10)

 

STARTING TIME: the time when simulation starts.

ENDING TIME: the time when simulation ends.

SDE TIME STEP: the maximum time step of the SDE numerical integrator, it may be set for adaptive methods to decrease memory requirements

RANDOM SEED: a random number generated by PC time, which is used to produce a series of uniformly distributed random numbers.

CUSTOMIZED SEED: a user specified number, which is used to produce a series of uniformly distributed random numbers.

NUMBER OF TRIALS: the number of multiple trials to be run.

EPSILON: minimum number of molecules, both reactant and product species, required for approximation as a continuous Markov process.

LAMBDA: minimum rate of reaction required for approximation to a continuous Markov process.

MSR TOLERANCE: maximum allowed effect of executing multiple slow reactions per numerical integration of the SDEs.

SDE TOLERANCE: maximum allowed value of the drift and diffusion errors.

OUTPUT OPTIONS: use output time interval.

 

For more information, please visit http://hysss.sourceforge.net/publications.shtml.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


 

           Figure 10. Parameters required by adaptive time step hybrid solver

·                    To view simulation results

 

To intuitively display the results, Virtual Cell displays the trajectories for a single stochastic run and histograms for Monte Carlo simulation. When results status becomes “yes” in the simulation list panel, users can click the “Results’ button to view results. The results will be displayed in either plots or table by selecting the “Show plots” and “Show data” buttons respectively.

 

a)      A trajectory for a single stochastic run shows a stochastic variable or probability of a jump process changes against time (Figure 11 and Figure 12).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


                       Figure 11. Trajectories of stochastic variables

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


          Figure 12. Stochastic variables change against time listed in a table

b)      A histogram for Monte Carlo simulation shows the frequency of a stochastic variable being a certain value at a specific time point after multiple trials (Figure 13 and Figure 14).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


              Figure 13. Histogram of a stochastic variable after Monte Carlo simulation

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


              

                  Figure 14. Histograms of stochastic variables listed in a table

 

3.   References

 

[1] Gillespie, D.T., A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions, J Comp. Phys., 1976, 22:403-434.

 

[2] Gibson, M. and Bruck, J., Efficient exact stochastic simulation of chemical systems with many species and channels, J Phys. Chem. A, 2000, 104:1876-1889.

 

[3] Salis, H. and Kaznessis , Y., Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions', J. Chem. Phys., 2005, 122: 054103.

 

[4] Salis, H., Sotiropoulos, V. and Kaznessis, Y., Multiscale Hy3S: Hybrid stochastic simulation for supercomputers, BMC Bioinformatics, 2006, 7:93.