Simulating the Dynamics of Complex Biological Systems

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