1 / 26

Stochastic Simulation of Biological Systems

Stochastic Simulation of Biological Systems. Chemical Reactions. Reactants  Products m 1 R 1 + m 2 R 2 + ··· + m r R r – ! n 1 P 1 + n 2 P 2 + ··· + n p P p. Example. 3 reactions, 2 species: Y 1 ! 2 Y 1 (1) Y 1 + Y 2 ! 2 Y 2 (2)

astra-santos
Download Presentation

Stochastic Simulation of Biological Systems

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Stochastic Simulation of Biological Systems

  2. Chemical Reactions • Reactants  Products • m1R1 + m2R2 + ···+ mrRr–! n1P1 + n2P2 + ···+ npPp

  3. Example • 3 reactions, 2 species: Y1 ! 2Y1 (1) Y1+Y2! 2Y2 (2) Y2 ! 0 (3) • Lotka-Volterra

  4. Representations (I) • Y1 Y2 1 1 0 = A 2 -1 1 3 0 -1 • A = Net Effect Matrix

  5. Example • x(t0) = (10 5) • Reaction (2) fires at time t1 > t0 • Then x(t1) = x(t0) + r.A = x(t0) + (0 1 0).A = (10 5) + (-1 1) = (9 6)

  6. Representations (II) • SBML • Electronic XML-based system • Represent models in a standard format • Sharing models • Transferring models between different software systems

  7. <reaction name="Protein dimerisation" reversible="false"> <listOfReactants> <specieReference specie="P" stoichiometry="2"/> </listOfReactants> <listOfProducts> <specieReference specie="P2" stoichiometry="1"/> </listOfProducts> <kineticLaw formula="k4*0.5*P*(P-1)"/> <listOfParameters> <parameter name="k4" value="1"/> </listOfParameters> </reaction>

  8. Simulating the LV model • ODEs • Use 2 assumptions (i) Continuous quantities (not integers) (ii) System behaves deterministically • Assumptions are often appropriate….. • ….but not always.

  9. Y1 ! 2Y1 (1) Y1+Y2! 2Y2 (2) Y2 ! 0 (3) • c1, c2, and c3 are the rate constants for reactions 1, 2, and 3 respectively.

  10. Start at time t0 and state x0 • Simulate time of first reaction event • ‘Hazards’ of the three reactions are: h1=c1y1 h2=c2y1y2 h3=c3y2 • Total hazard function is h0=h1+h2+h3 • h0 is constant during inter-event time periods • T ~ Exponential(h0)

  11. Now select which of the 3 reactions has occurred. • Probability of choosing reaction i is given by: hi / h0

  12. Update system time: t := t0 + T • Update state vector: x := x0 + ri.A

  13. Continue moving forwards in time • Gillespie’s Stochastic Simulation Algorithm (the SSA) • A single run generates a time vector (tvec) and a state matrix (xmat) • Algorithm must be run many times • Calculate statistical quantities • E.g. mean and standard deviation of Y1 species over the first 50 seconds

  14. SSA is an exact algorithm • Each reaction event simulated individually • Algorithm must run thousands of times • Very time consuming for large, complex systems • How to speed things up?

  15. Time Discretization • Divide time interval into discrete chunks of width Δt. • Go from t to t+Δt in one “leap” • More than one reaction may occur in the interval [t, t+Δt) • Simulate which reactions have occurred in this interval.

  16. Time interval [0, t] • Let X be the number of times reaction 1 fires in [0, t] • Divide interval into N discrete chunks of width δt • N large, δt small • Prob(Reaction 1 fires once in time δt)  h1.δt • Prob(Reaction 1 fires more than once in time δt) is negligible (relative to δt) • Then X ~ Bin(N, h1.δt) = Bin(N, h1.(t/N)) • E(X) = h1.t • As N  , δt  0 , X~Poisson(h1.t)

  17. “Tau-Leap” Algorithm • Start at time t0, state x0 • Select discrete step size,  • Move forwards to time t0+  • Simulate the number of firings of reaction i during the interval as a Poisson(hi.) (i=1,…n) • Update the system state according to the reactions which have occurred. • Now advance to t0+ 2, t0+ 3, ……

  18. Tau-Leap is faster than the SSA • But some accuracy is sacrificed • Relies on a sensible choice of 

  19. Hybrid Algorithms • Divide the system species and reactions into discrete and continuous regimes • The discrete regime contains species with low concentrations and reactions which fire infrequently • The continuous regime contains species with high concentrations and reactions which fire frequently

  20. Use different methods for simulating the two regimes. • Simulate the continuous regime using deterministic or stochastic differential equations • Simulate the discrete regime using the SSA • The two regimes must be syncronized • Some species will be involved in both regimes

  21. Measuring Accuracy • The accuracy of an approximate algorithm (leaping or hybrid) can be assessed by comparing results with SSA • Run SSA and approximate algorithm 1,000 times for the same system, using the same rate constants and initial concentrations • Do the results from the approximate algorithm closely resemble the results from the SSA?

  22. The Test Suite for Stochastic Simulators • Available at www.calibayes.ncl.ac.uk • Contains stochastic simulation results for some simple biological systems • Results generated using analytic expressions for species means and standard deviations. • Can be used to measure the correctness of a stochastic simulator

  23. Testing stochastic simulators is not straightforward. • Even if a simulator is working perfectly, it will not produce exactly the results contained in the test suite. • Need statistical method to compare the simulator results with the test suite results. • E.g. if the concentration of species X after 10 seconds varies by 1% from the test suite value, is this evidence of a problem with the simulator?

More Related