Stochastic Solvers

When to Use Stochastic Solvers

The stochastic simulation algorithms provide a practical method for simulating reactions that are stochastic in nature. Models with a small number of molecules can realistically be simulated stochastically, that is, allowing the results to contain an element of probability, unlike a deterministic solution.

Model Prerequisites for Simulating with a Stochastic Solver

Model prerequisites include:

• All reactions in the model must have their KineticLaw property set to MassAction.

• If your model contains events, you can simulate using the stochastic ssa solver. Other stochastic solvers do not support events.

• Your model must not contain doses. No stochastic solvers support doses.

Additionally, if you perform an individual or population fitting on a model whose Configset object specifies a stochastic solver and options, be aware that during the fitting SimBiology® temporarily changes:

• SolverType property to the default solver of ode15s

• SolverOptions property to the options last configured for a deterministic solver

What Happens During a Stochastic Simulation?

During a stochastic simulation of a model, the software ignores any rate, assignment, or algebraic rules if present in the model. Depending on the model, stochastic simulations can require more computation time than deterministic simulations.

Tip

When simulating a model using a stochastic solver, you can increase the LogDecimation property of the configset object to record fewer data points and decrease run time.

Stochastic Simulation Algorithm (SSA)

The Chemical Master Equation (CME) describes the dynamics of a chemical system in terms of the time evolution of probability distributions. Directly solving for this distribution is impractical for most realistic problems. The stochastic simulation algorithm (SSA) instead efficiently generates individual simulations that are consistent with the CME, by simulating each reaction using its propensity function. Thus, analyzing multiple stochastic simulations to determine the probability distribution is more efficient than directly solving the CME.

• This algorithm is exact.

• Because this algorithm evaluates one reaction at a time, it might be too slow for models with a large number of reactions.

• If the number of molecules of any reactants is huge, it might take a long time to complete the simulation.

Explicit Tau-Leaping Algorithm

Because the stochastic simulation algorithm might be too slow for many practical problems, this algorithm was designed to speed up the simulation at the cost of some accuracy. The algorithm treats each reaction as being independent of the others. It automatically chooses a time interval such that the relative change in the propensity function for each reaction is less than your error tolerance. After selecting the time interval, the algorithm computes the number of times each reaction occurs during the time interval and makes the appropriate changes to the concentration of various chemical species involved.

• This algorithm can be orders of magnitude faster than the SSA.

• You can use this algorithm for large problems (if the problem is not numerically stiff).

• This algorithm sacrifices some accuracy for speed.

• This algorithm is not good for stiff models.

• You need to specify the error tolerance so that the resulting time steps are of the order of the fastest time scale.

Implicit Tau-Leaping Algorithm

Like the explicit tau-leaping algorithm, the implicit tau-leaping algorithm is also an approximate method of simulation designed to speed up the simulation at the cost of some accuracy. It can handle numerically stiff problems better than the explicit tau-leaping algorithm. For deterministic systems, a problem is said to be numerically stiff if there are “fast” and “slow” time scales present in the system. For such problems, the explicit tau-leaping method performs well only if it continues to take small time steps that are of the order of the fastest time scale. The implicit tau-leaping method can potentially take much larger steps and still be stable. The algorithm treats each reaction as being independent of others. It automatically selects a time interval such that the relative change in the propensity function for each reaction is less than the user-specified error tolerance. After selecting a time interval, the algorithm computes the number of times each reaction occurs during the time interval and makes the appropriate changes to the concentration of various chemical species involved.

• This algorithm can be much faster than the SSA. It is also usually faster than the explicit tau-leaping algorithm.

• You can use this algorithm for large problems and also for numerically stiff problems.

• The total number of steps taken is usually less than the explicit-tau-leaping algorithm.

• This algorithm sacrifices some accuracy for speed.

• There is a higher computational burden for each step as compared to the explicit tau-leaping algorithm. This leads to a larger CPU time per step.

• This method often dampens perturbations of the slow manifold leading to a reduced state variance about the mean.

References

 Gibson M.A., Bruck J. (2000), “Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels,” Journal of Physical Chemistry, 105:1876–1899.

 Gillespie D. (1977), “Exact Stochastic Simulation of Coupled Chemical Reactions,” The Journal of Physical Chemistry, 81(25): 2340–2361.

 Gillespie D. (2000), “The Chemical Langevin Equation,” Journal of Chemical Physics, 113(1): 297–306.

 Gillespie D. (2001), “Approximate Accelerated Stochastic Simulation of Chemically Reacting Systems,” Journal of Chemical Physics,115(4):1716–1733.

 Gillespie D., Petzold L. (2004), “Improved Leap-Size Selection for Accelerated Stochastic Simulation,” Journal of Chemical Physics, 119:8229–8234

 Rathinam M., Petzold L., Cao Y., Gillespie D. (2003), “Stiffness in Stochastic Chemically Reacting Systems: The Implicit Tau-Leaping Method,” Journal of Chemical Physics, 119(24):12784–12794.

 Moler, C. (2003), “Stiff Differential Equations Stiffness is a subtle, difficult, and important concept in the numerical solution of ordinary differential equations,” MATLAB News & Notes.

Related Examples 