Stochastic Simulation Methods for Chemical Reactions |

1. Stochastic Simulation Hauptseminar 11.07.2007 Ludwig Geistlinger

2. Gliederung • 1. Problemstellung • 2. Deterministische Methode • 3. Stochastische Methode • 4. Effizienzverbesserung • 5. Anwendung

3. 1. Problemstellung

4. 1. Problemstellung Gegeben: → räumlich homogenes Volumen V → Nchemisch aktive Spezies Si ( i = 1, …, N ) Xi := aktuelle Anzahl von Molekülen der Art Siin V → Mchemische Reaktionen Rμ( μ= 1, …, M ) → Rμdurch Reaktionskonstante kμcharakterisiert

5. 1. Problemstellung C B + N = 4 D A A B Si = { A, B, C, D } D + A Xi = { 8, 6, 5, 5} C A A C D Rμ= { A ↔ B, 2A ↔ C, A+D ↔ 2A } D C B C B A + D A B kμ= { k1 , k2, …, k6 } A B

6. 1. Problemstellung Zielstellung: Bestimmung bzw. Simulation der Zeitevolution aller Xi Xi(t) Konzentration c cEq Zeit t

7. 2. Deterministische Methode 2.1 Exkurs Differentialgleichungen 2.2 Zwei Beispiele 2.3 Probleme dieser Methode

8. 2. Deterministische Methode 2.1 Exkurs Differentialgleichungen → Ableitung als Änderungsrate (hier: Änderung der Konzentration pro Zeit) → Anfangswertprobleme mit Diff.glg. erster Ordnung: gegeben:(A) Diff.glg.: → Allgemeine Lösung: (B) Anfangsbedingung: x(t0) = x0 → ergibt spezielle Lösung

9. 2. Deterministische Methode 2.2 Zwei Beispiele (1) mit X(0) = X0 X(t) … kontinuierliche Fkt., die die Anzahl der Moleküle der Species X in V zum Zeitpunkt t repräsentiert (A) Aufstellen der Differentialgleichung: (B) Lösen der Differentialgleichung: (C) Einsetzen der Anfangsbedingung:

10. 2. Deterministische Methode (2) → Aufstellen der Differentialgleichungen:

11. 2. Deterministische Methode 2.3 Probleme der deterministischen Methode → Keine Berücksichtigung von Korrelationen → Keine Berücksichtigung von Fluktuationen → Problembereich Differentialgleichungen (Berechnung, Eindeutigkeit, Existenz) … Alternative(n) ?

12. 3. Stochastische Methode 3.1 Die Idee 3.2 Der Algorithmus 3.3 Simulation mit COPASI

13. 3. Stochastische Methode 3.1 Die Idee → det. Methode: Reaktionskonstanten ↔ Reaktionsraten N gekoppelte Differentialgleichungen → stochastische Methode: Reaktionskonstanten ↔ Reaktionswahrscheinlichkeiten Random Walk auf einer N-dimensionalen Markov-Kette Simulation über eineMaster-Equation

14. 3. Stochastische Methode Master Equation Bsp.:

15. 3. Stochastische Methode Fundamentalhypothese: kμ … Wahrscheinlichkeit, dass eine bestimmte Kombination von Reaktanden im nächsten Zeitintervall gemäß Rμ reagiert hμ … Anzahl der unterschiedlichen Kombinationen von Reaktanden für Rμ zum Zeitpunkt t in V → hμkμ… Wahrscheinlichkeit, dass eine Rμ Reaktion im nächsten Zeitintervall in V auftritt

16. 3. Stochastische Methode Dichtefunktion der Reaktionswahrscheinlichkeit P(τ, μ) … Wahrscheinlichkeit zum Zeitpunkt t, dass die nächste Reaktion in V im Zeitintervall [ t , t + τ] auftritt und eine Reaktion des Typs Rμ ist P(„Rμtritt auf“) P(„Keine Reaktion“)

17. 3. Stochastische Methode 3.2 Der Algorithmus Schritt 1 (Initialisierung) t = 0 : Xi = Xi(0) kμ → hμkμ Schritt 2 (Monte Carlo) → generiere Zufallspaar (τ, μ) aus P(τ, μ) Schritt 3 (Evolution) t = t + τ : Xi = Xi +∆Xi→ Berechne hμkμneu Schritt 4 (Terminierung) t > tstop oder alle hμ = 0

18. 3. Stochastische Methode 3.2 Der Algorithmus Schritt 1 (Initialisierung) t = 0 : Xi = Xi(0) kμ= kμ(0) → hμkμ Schritt 2 (Monte Carlo) → generiere Zufallspaar (τ, μ) aus P(τ, μ) Schritt 3 (Evolution) t = t + τ : Xi = Xi -∆Xi→ Berechne hμkμneu Schritt 4 (Terminierung) t > tstop oder alle hμ = 0

19. 3. Stochastische Methode Generierung eines Zufallpaares (τ, μ) aus P(τ, μ) „Konditionieren“: → τist stetig (0 ≤ τ<∞) und μist diskret (μ= 1, 2, …, M) ! Substitution aμ= hμkμ

20. 3. Stochastische Methode Generierung eines Zufallpaares (τ, μ) aus P(τ, μ)

21. 3. Stochastische Methode 3.2 Der Algorithmus Schritt 1 (Initialisierung) t = 0 : Xi = Xi(0) kμ= kμ(0) → hμkμ Schritt 2 (Monte Carlo) → generiere Zufallspaar (τ, μ) aus P(τ, μ) Schritt 3 (Evolution) t = t + τ : Xi = Xi -∆Xi→ Berechne hμkμneu Schritt 4 (Terminierung) t > tstop oder alle hμ = 0

22. 3. Stochastische Methode Evolution → Verwendung des generierten Zufallpaares (τ,μ): (1) setze: t = t + τ (2) simuliere Rμdurch Änderung aller involvierten Xi Beispiel: Rμ: S1 + S2 → 2S3 →X1 = X1 – 1 →X2 = X2– 1 → X3 = X3+ 2 →Berechne alle aμneu, in denen die involvierten Xi als Reaktanden auftreten

23. 3. Stochastische Methode 3.2 Der Algorithmus Schritt 1 (Initialisierung) t = 0 : Xi = Xi(0) kμ= kμ(0) → hμkμ Schritt 2 (Monte Carlo) → generiere Zufallspaar (τ, μ) aus P(τ, μ) Schritt 3 (Evolution) t = t + τ : Xi = Xi -∆Xi→ Berechne hμkμneu Schritt 4 (Terminierung) t > tstop oder alle hμ = 0

24. 3. Stochastische Methode 3.3 Simulation mit COPASI COPASI … COmplex PAthway SImulator Y(t) X(t)

25. 3. Stochastische Methode • DM (Direct Method) • … direkte Generierung von μundτmit Inversionsmethode → benötigt 2 Zufallszahlen pro Iteration • FDM (First Reaction Method) • Generiere für jede Reaktion τμ • Durchlaufe Reaktion mit niedrigstem τμ • → benötigt M Zufallszahlen pro Iteration

26. 4. Effizienzverbesserung 4.1 Verbesserungen im Überblick 4.2 LDM 4.3 Tau-Leaping

27. 4. Effizienzverbesserung • 4.1 Verbesserungen im Überblick • NRM (Next Reaction Method): Erstellung eines dependency graph Verwendung einer indizierten priority queue Wiederverwendung von Zufallszahlen • ODM (Optimized Direct Method): Umordnung der Reaktionen nach Häufigkeit Bedarf einer Prä-Simulation • SDM(Sorting Direct Method): Dynamisches ordnen der Reaktionen • 4. LDM (Logarithmic Direct Method)

28. 4. Effizienzverbesserung → gibt an, welcheaμverändert werden müssen (geg. Rμ)

29. 4. Effizienzverbesserung • 4.1 Verbesserungen im Überblick • NRM (Next Reaction Method): Erstellung eines dependency graph Verwendung einer indizierten priority queue Wiederverwendung von Zufallszahlen • ODM (Optimized Direct Method): Umordnung der Reaktionen nach Häufigkeit Bedarf einer Prä-Simulation • SDM(Sorting Direct Method): Dynamisches ordnen der Reaktionen • 4. LDM (Logarithmic Direct Method)

30. 4. Effizienzverbesserung Label Index Zeit Zeiger Indexstruktur → Effiziente Verwaltung der absoluten Reaktionszeiten

31. 4. Effizienzverbesserung • 4.1 Verbesserungen im Überblick • NRM (Next Reaction Method): Erstellung eines dependency graph Verwendung einer indizierten priority queue Wiederverwendung von Zufallszahlen • ODM (Optimized Direct Method): Umordnung der Reaktionen nach Häufigkeit Bedarf einer Prä-Simulation • SDM(Sorting Direct Method): Dynamisches ordnen der Reaktionen • 4. LDM (Logarithmic Direct Method)

32. 4. Effizienzverbesserung 4.2 LDM (Logarithmic Direct Method) → Kostenintensivster Schritt der Methode ist das Auffinden der nächsten ablaufenden Reaktion im System 1. 2. → Propensities werden beinahe zweimal aufsummiert → Suchtiefe im zweiten Schritt liegt in O(M)

33. 4. Effizienzverbesserung Idee „Binäre Suche“ → Speichern aller M partiellen Summen in einem Array → Binäre Suche: Wiederholtes Aufteilen des Suchintervalles in 2 Hälften Test:subtotal[μ–1] ≤ key < subtotal[μ] mit key := r2a → Suchtiefe für binäre Suche liegt bei O(logM)

34. 4. Effizienzverbesserung 4.3 Tau-Leaping Approximation →Überspringen von Reaktionsereignissen Dμ(τ) … Anzahl der Läufe durch Rμ im Intervall [t, t + τ] → Wähleτ: - größtmöglich - jedoch klein genug, dass sich aμ nicht signifikant ändern (leap condition) Approx.: Update:

35. 5. Anwendung

36. 5. Anwendung 24 h Simulation einer Säugerzelle → zelluläre „Uhr“ durch oszillierende Protein Feedback Loops → Einteilung der Reaktionen in Reaktionsklassen (Bedarf einer zusätzlichen Zufallszahl !) → Analyse der Stabilität der „Uhr“ gegenüber: (a) molekularem Rauschen (b) Knockout von Schlüsselgenen / -Proteinen

37. Zusammenfassung

38. Zusammenfassung Bestimmung / Simulation der Evolution von N chemischen Spezies in Reaktionsnetzwerk von M Reaktionswegen Deterministische Methode: Lösung N gekoppelter Differentialgleichungen Stochastische Methode: 4-Schritte-Simulation Dichtefunktion der Reaktionswahrscheinlichkeit Generierung zweier Zufallszahlen Verschiedene Möglichkeiten der Effizienzverbesserung Simulation komplexer Netzwerke

39. Literatur

40. Literatur Differentialgleichungen: [1] D. Jordan, P. Smith: „Mathematical Techniques“, Oxford Univ. Press, 1994 [2] O. Levenspiel: „Chemical Reaction Engineering“, John Wiley & Sons, 1999 Stochastische Methode:[3] D. Gillespie: „A General Method for Numerically Simulating the Stochastic TimeEvolution of CoupledChemical Reactions“, Jrn. of Comp. Phys. 22, 403-434, 1976 [4]D. Gillespie: „Exact Stochastic Simulation of Coupled Chemical Reactions“, Jrn. Of Phys. Chem. 81, 2340-2361, 1977 Effizienzverbesserung: [5] M. Gibson, J. Bruck: „Efficient exact stochastic simulation of chemical systems with many species and channels“, Jrn. of Phys. Chem. 104, 1876-1889, 1999 [6] Y. Cao, H. Li, L. Petzold: „Efficient formulation of the stochastic simulation algorithm for chemically reacting systems“, J. Chem. Phys. 121, 4059-4067, 2004 [7] H. Li, L. Petzold : „Logarithmic Direct Method for Discrete Stochastic Simulation of Chemically Reacting Systems“, Jrn. of Chem. Phys., 2006 [8] Y. Cao, D. Gillespie, L. Petzold: „Efficient Step Size Selection for the Tau-Leaping Simulation Method“, Jrn. Of. Chem. Phys. 124, 044109, 2006 Anwendung [9] D. Forger, C. Peskin: „Stochastic Simulation of the Mammalian Circadian Clock“, Proceedings of the National Academy of Sciences, 2005

41. ENDE