Model Simulation

Simulating a Model

Simulate a model by providing the model object as an input argument to the sbiosimulate function.

When you simulate a model, you can return results (time points, state data, and state names) in two forms:

A SimData object also includes metadata such as the types and names for the logged states, the configuration set used during simulation, and the date of the simulation. It is a convenient way of keeping time data, state data, and associated metadata together. A SimData object has associated properties and methods, which you can use to access and manipulate the data.

For more information on simulating a model, see Simulate a Model and View Results.

Plotting Simulation Results

If you return time and state data from a simulation in three output arguments, you can use these arguments as inputs to the plot function to view your results. For more information, see sbiosimulate.

If you return time and state data from a simulation in a SimData object, you can use the SimData object as an input to the sbioplot function to view your results.

For more information on plotting simulation results, see Simulate a Model and View Results.

Interpreting Simulation Results

After running a simulation, you may see negative amounts or concentrations for species in the results plot or data array. These negative values can be either:

  • Slightly negative due to numerical noise introduced by the simulation process. In this case, you can interpret these values as 0.

  • Significantly negative due to the dynamics in your model not being physical, that is, the dynamics in the system are driving a particular species to be negative. In this case, examine your reaction rate expressions to ensure they implement correct dynamics.

Configuring Stop Time and Other Simulation Settings

A model has a configuration set (Configset object) associated with it to control the simulation. You can edit the properties of a Configset object to control all aspects of the simulation, including:

To view the Configset object, provide the model object as an input argument to the getconfigset method.

To edit the properties of a Configset object, use the set method.

For more information on viewing and editing the stop time and other simulation settings, see Simulate a Model and View Results.

Choosing a Simulation Solver

To simulate a model, the SimBiology® software converts a model to a system of differential equations. It then uses a solver function to compute solutions for these equations at different time intervals, giving the model's states and outputs over a span of time.

Available solvers are:

  • ODE Solvers — These include Nonstiff Deterministic Solvers and Stiff Deterministic Solvers. The solver functions implement numerical integration methods for solving initial value problems for ordinary differential equations (ODEs). Beginning at the initial time with initial conditions, they step through the time interval, computing a solution at each time step. If the solution for a time step satisfies the solver's error tolerance criteria, it is a successful step. Otherwise, it is a failed attempt; the solver shrinks the step size and tries again. For more information, see ODE Solvers in the MATLAB® Mathematics documentation.

  • SUNDIALS Solvers — At a fundamental level the core algorithms for the SUNDIALS solvers are similar to those for some of the solvers in the MATLAB ODE suite and work as described above in ODE Solvers. For more information, see SUNDIALS Solvers.

  • Stochastic Solvers — Use with models containing a small number of molecules. Stochastic solvers include stochastic simulation algorithm, explicit tau-leaping algorithm, and implicit tau-leaping algorithm. For more information, see Stochastic Solvers.

SUNDIALS Solvers

SUNDIALS (Suite of Nonlinear and Differential/Algebraic Equation Solvers) are part of a freely available third-party package developed at Lawrence Livermore National Laboratory. All other ODE solvers used for simulation of SimBiology models, such as ode45 and ode15s, are part of the MATLAB ODE suite.

When you specify sundials for the solver, the software chooses one of two SUNDIALS solvers, CVODE or IDA, as appropriate for your model:

  • CVODE is a solver for systems of ODEs, both nonstiff and stiff. This is used when a model has no algebraic rules.

  • IDA is a differential-algebraic equation (DAE) solver, used when one or more algebraic rules are present.

For more information on the SUNDIALS solvers, see http://www.llnl.gov/casc/sundials/description/description.html.

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.

Advantage

  • This algorithm is exact.

Disadvantages

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

Advantages

  • 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).

Disadvantages

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

Advantages

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

Disadvantages

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

More About

References

[1] 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.

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

[3] Gillespie D. (2000), "The Chemical Langevin Equation," Journal of Chemical Physics, 113(1): 297–306.

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

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

[6] 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.

[7] 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.

Ensemble Runs of Stochastic Simulations

Because stochastic simulations rely on an element of probability, sequential runs produce different results. Therefore, multiple stochastic runs are needed to determine the probability distribution of the simulation results.

Ensemble runs perform multiple simulations of a model using a stochastic solver. They let you gather data from multiple stochastic runs of the model so you can compare and analyze fluctuations in the behavior of a model over repeated stochastic simulations.

Running Ensemble Simulations

The following functions let you perform and analyze ensemble runs at the command line:

  • sbioensemblerun — Perform a stochastic ensemble run of the MATLAB model object.

  • sbioensembleplot — Show a 2-D distribution plot or a 3-D shaded plot of the time varying distribution of one or more specified species.

  • sbioensemblestats — Get mean and variance as a function of time for all the species in the model used to generate ensemble data by running sbioensemblerun.

Was this topic helpful?