1 / 50

Simulating the Dynamics of Complex Biological Systems

Simulating the Dynamics of Complex Biological Systems. Maria J Schilstra Biocomputation Research Group University of Hertfordshire, Hatfield, UK Stephen R. Martin Physical Biochemistry MRC National Instritute for Medical Resaerch, London, UK. Computer simulation.

chanel
Download Presentation

Simulating the Dynamics of Complex 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. Simulating the Dynamics of Complex Biological Systems Maria J Schilstra Biocomputation Research Group University of Hertfordshire, Hatfield, UK Stephen R. Martin Physical Biochemistry MRC National Instritute for Medical Resaerch, London, UK

  2. Computer simulation • “An attempt to model a real-life or hypothetical situation on a computer, so that it can be studied to see how the system works. By changing variables, predictions may be made about the behaviour of the system.” (Wikipedia)

  3. Running example

  4. Microtubule Dynamic Instability • Michison & Kirschner (1984): “Microtubules in vitro coexist in growing and shrinking populations which interconvert rather infrequently”

  5. G-end (net growth) S-end (rapid shrinkage) “Rescue” “Catastrophe” Cartoon

  6. Microtubule Dynamic Instability • Michison & Kirschner (1984): “microtubules in vitro coexist in growing and shrinking populations which interconvert rather infrequently” • Horio & Hotani (1986): • Average lifetime of growing microtubules: 3 min • Growth rate: 0.6 m/min (16 subunits/s) • Average lifetime of shrinking microtubules: 18 s • Shrinkage rate: 0.13 m/s (220 subunits/s)

  7. Model Notation

  8. G Growing end S Shrinking end TuT Tubulin-GTP TuD Tubulin-GDP n Number of Tu in MT (TuM) MT Microtubule TuD TuT 3 G S Model notation Chemical Reaction Notation “Petri-net” notation Gn + TuT  Gn+1 Gn Sn Sn Sn-1 + TuD Sn + 3 TuT  Gn+3

  9. Petri-net? A  B “Place” (pool, store, state) Contains “items” (tokens, molecules, particles) A B “Transition” (reaction, process) Can “fire”; has an average firing rate

  10. Transition firing? A  B A B X

  11. Transition firing? 2A  B A  B2 2 2 A B A B A  B+X A+X  B B A A B A A

  12. k (1/s) A  B k (L/(mole.s)) A B A+X  B A B A Transition firing rate? J = rate at which B appears J = k.[A] (mol/(L.s)) J = k.nA(items/s) nA = [A].NA.Vol J = k.[A][X] mol/(L.s)) J = k.nA.[X] (items/s) nA = [A].NA.Vol

  13. TuD TuT 3 G S Microtubule dynamic instability model

  14. Transition firing TuD TuT 3 G S

  15. Transition firing TuD TuT 3 G S

  16. Counting TuM (bound Tu) TuD TuT 3 G S TuM

  17. Counting TuM (bound Tu) TuD TuT 3 G S TuM

  18. Firing rates TuD TuT J = kSG nS[TuT]3 3 kSG G S kGG kSS J = kGG nG[TuT] J = kSS nS kGS J = kGS nG TuM

  19. Parameter values Initial conditions TuD TuT J = kSG [S][TuT]3 3 kSG G S kGG kSS J = kGG [G][TuT] J = kSS [S] kGS J = kGS [G] TuM

  20. Stochastic Simulation

  21. Stochastic? • “Random or probabilistic but with some direction. For example the arrival of people at a post office might be random but average properties (such as the queue length) can be predicted.” (http://www.cs.ucl.ac.uk/staff/W.Langdon/gpdata/glossary.html) • Stochastic simulation: uses a random number generator to produce one or more possible time courses.

  22. Stochastic simulation • Assess which transitions are “enabled” (can fire) • Use a “weighted lottery” to determine which enabled transition will fire first (Tfirst), and when (tnext) • Let transition Tfirst fire at time tnext • Repeat from 1 as long as there are enabled transitions, and tnext < tstop

  23. Firing rates JGG = 1.6x106 x 50x10-9 x 1 = 0.08 s-1 (i.e. lifetime Gn = 12.5 s) JGS = 0.0056 x 1 = 0.0056 s-1 (i.e. lifetime G = 180 s) TuD TuT 3 kSG G S kGG kSS J = kGG [G][TuT] kGS J = kGS [G] TuM

  24. Probability that reaction has occurred(Cumulative distribution functions) Gn Gn+1 (lifetime Gn: 12.5 s) G  S (lifetime G: 180 s)

  25. The weighted lottery rGG = 0.76 rGS = 0.40 tGG = 18 s tGS = 92 s

  26. G  G transition fires 0 20 40 60 18 s time (s) TuD TuT 3 kSG G S kGG kSS kGS TuM

  27. The weighted lottery (continued) rGG = 0.98 rGS = 0.15 tGG = 49 s tGS = 30 s

  28. G  S transition fires 0 20 40 60 18 s 30 s time (s) TuD TuT 3 kSG G S kGG kSS kGS TuM

  29. Other transitions enabled 0 20 40 60 18 s 30 s time (s) TuD TuT 3 kSG G S kGG kSS kGS TuM

  30. Trajectory of TuM[TuT]constant (10 M) - 150,000 events GG Delta TuM SS

  31. [TuT] variable2 million events TuD TuT MT length MT length (m) [TuT], [TuD] (M)

  32. G S G S G S n n n More than one MT end TuT TuD Method 1: modelling all items individually

  33. Trajectories(6 million events) MT length [TuT] (M) MT length (m) [TuT]

  34. Improving efficiency

  35. After most recent transition firing at time t: For each transition Ti: Compute current firing rate Ji for Ti Draw two numbers, r1 and r2 Compute J = sum of all J Compute firing time tfire from J,r1, and t Construct roulette wheel (pie chart) for all J Spin wheel over r2x 360º, find associated T Set t to tnext and fire T “First Reaction Algorithm” (Gillespie,1976)

  36. Algorithms • Accurate (no further assumptions) • First Reaction (Gillespie, 1976) • Next Reaction (Gibson & Bruck, 2000) • Priority queue; exploits independence in system • Approximate • Tau-leaping (Gillespie, 2001) • Assumes that “leap condition” is satisfied (leap condition: negligible changes in firing probability over interval Tau) • Multiple firings per step • Chemical Langevin • Assumes that 1) leap condition is satisfied; 2) All current firing rates are much larger than 1/Tau (“Large yet small”) • Like deterministic trajectory with superimposed noise

  37. Deterministic Simulation

  38. Deterministic? Average of many trajectories converge to a smooth path, that is entirely determined by the initial conditions

  39. Assessing the average trajectory • Stochastic • Many trajectories • Many particles (larger volume) • Deterministic • Numerical integration of coupled ODEs (ordinary differential equations)

  40. ODEs • ODE: defines how system variables change with another variable • General form: • d[X]/dt = v+(X) – v-(X) • v+(X) = J+(X) and v-(X) = J-(X) Rate at which X is consumed Change in [X] over very small time interval dt Rate at which X is produced

  41. Composing ODEs from firing rates TuD TuT JSG = kSG [S][TuT]3 v+(S) = JSS + JGS 3 kSG G S kGG kSS v(S) = JSG + JSS kGS JGG = kGG [G][TuT] JSS= kSS [S] JGS = kGS [G] TuM

  42. The full system Growing ends Shrinking ends Free Tu-GTP Free Tu-GDP Bound Tu

  43. Predicting average trajectories:Numerical integration • Set component concentrations to their initial values e.g. [TuT] = 50.0 M, [G] =1.6 nM, [S]=[TuM]=[TuD]=0.0 • For each component, compute concentration change over small time interval t e.g. [TuM] = [TuM]t+t – [TuM]t = {kGS[G]t [TuT]t- kSG[S]t}x t • Solve all concentrations at t + t e.g. [TuM]t+t = 0.0 M + (1.6x106 M-1s-1 x 50 M x 1.6 nM – 220 s-1 x 0.0 M) x 1 ms = 0.128 nM • Set t to t+t, and [X]t to [X]t+t, and repeat from 2 until end condition is fulfilled e.g. t > tmaxor no more significant changes in any concentration

  44. Result

  45. Things to keep in mind • Main assumption: firing rates (J) remain constant over t • Therefore: [X] used to compute J must be negligible, so that t must be made sufficiently small • Sophisticated algorithms calculate most efficient t before each time step • Assessment of efficient t uses “accuracy” parameter • Approach invalid when J not constant over t (within set limits) • Causes: accuracy too low • Warning signs: wildly oscillating trajectories, negative concentrations, very different results for greater accuracy

  46. Algorithms • Forward (explicit) • Solve: • [G]t+t = [G]t + { kSG[S]t – kGS[G]t [TuT]t3 }x t • Examples: Euler, Runge-Kutta, Adams-Bashford • Use: for many ODE systems (easier to implement, individual steps less computationally intensive than implicit methods) • Backward (implicit) • Solve: • [G]t= [G]t +t + { kSG[S]t +t – kGS[G]t +t[TuT]t +t3 }x t • Examples: Gear, Adams-Moulton • Use: for “stiff” systems (concentrations changing on very different time scales)

  47. Stochastic Deal with finite numbers of items Cheap assessment of inherent “noise” Assessment of average population behaviour expensive Easy to form mental picture Easy to incorporate in physical models Deterministic Deal with populations of infinite size No information about inherent noise Assessment of average population behaviour computationally cheap More difficult to grasp and use Integration with other model types not trivial Stochastic vs deterministic approaches

  48. Stochastic or deterministic? • “Well-stirred containers”  • Deterministic faster, results less confusing than stochastic • “Noise” (deviations from average) important for interpretation experimental results  • Many items: deterministic with superimposed variation; approximate stochastic (Tau-leap, chemical Langevin) • Few items: accurate stochastic (Gillespie, Gibson-Bruck) • Firing probabilities depend on geometry that changes after each event (e.g. detailed description of microtubule ends)  • Accurate stochastic, needs tailor-made software

  49. Software tools • Programming languages • C/C++, C#, Delphi , Java, Python, VBasic… • High level math/statistics software • Maple, Mathematica, MATLAB, Octave, R … • Dynamical Systems software • Berkeley-Madonna, Dymola, ModelMaker, XPP… • Dedicated chemical kinetics software • Deterministic: CellDesigner, DBSolve , E-Cell, Gepasi, Jarnac/JDesigner, Promot-DIVA, PySces, Virtual Cell… • Stochastic: Dizzy, StochKit… • Both: Copasi, NetBuilder’…

  50. Acknowledgement • Peter Bayley, Justin Molloy (NIMR, London) • The SBML Forum • Jack Correia • The Wellcome Trust, EPSRC, UH

More Related