|On this page…|
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:
Three separate arrays
One SimData object
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.
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.
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.
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:
Time units (TimeUnits property)
Maximum time step size (MaxStep property)
Data to record (RuntimeOptions property)
Dimensional analysis and unit conversion (CompileOptions property)
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.
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 (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.
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 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
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.
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.
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.
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.
 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.
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.
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.
For examples of simulating models, see: