Stochastic Simulation

Hauptseminar

11.07.2007

Ludwig Geistlinger

Gliederung

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

1. Problemstellung

1. Problemstellung

Gegeben:

→ räumlich homogenes Volumen V

→ N chemisch aktive Spezies Si ( i = 1, …, N )

Xi := aktuelle Anzahl von Molekülen der Art Si in V

→ M chemische Reaktionen Rμ ( μ = 1, …, M )

→ Rμ durch Reaktionskonstante kμ charakterisiert

A

D

C D

B

B

B

C

C

C

D A

A

B

D

A

B

A B

CD

AA

A

1. Problemstellung

Xi = { 8, 6, 5, 5}

N = 4

Si = { A, B, C, D }

Rμ = { A ↔ B, 2A ↔ C, A+D ↔ 2A }

kμ = { k1 , k2, …, k6 }+

+

+

1. Problemstellung

Zielstellung:

Bestimmung bzw. Simulation der Zeitevolution aller Xi

Zeit t

Kon

zent

ratio

n c

cEq

Xi(t)

2. DeterministischeMethode

2.1 Exkurs Differentialgleichungen

2.2 Zwei Beispiele

2.3 Probleme dieser Methode

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.:0

d

d cx

t

x

→ Allgemeine Lösung: cteAtx )(

(B) Anfangsbedingung: x(t0) = x0

→ ergibt spezielle Lösung

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

YX k

kXt

X

d

d(A) Aufstellen der Differentialgleichung:

(B) Lösen der Differentialgleichung: kteAtX )(

(C) Einsetzen der Anfangsbedingung: kteXtX 0)(

2. Deterministische Methode

(2)

2 ,2 ,

5

6

3

4

1

2

XXWZXYXRk

k

k

k

k

k

→ Aufstellen der Differentialgleichungen:

d

d

t

Y 21 YkXk

t

Z

d

dZkXk 4

232

1

t

W

d

d 265 2

1XkWXk

t

X

d

d 2654

2321 2

12 XkWXkZkXkYkXk

2.3 Probleme der deterministischen Methode

2. Deterministische Methode

→ Keine Berücksichtigung von Korrelationen

→ Keine Berücksichtigung von Fluktuationen

→ Problembereich Differentialgleichungen

(Berechnung, Eindeutigkeit, Existenz)

… Alternative(n) ?

3. StochastischeMethode

3.1 Die Idee3.2 Der Algorithmus3.3 Simulation mit COPASI

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 eine Master-Equation

3. Stochastische Methode

Master Equation

Bsp.: YXk

k

1

2

)|1,1()1(

)|1,1()1()'|,(

2

1

tYXPYk

tYXPXktYXP

M

XN

XN tRXRXPTtXXP N

111 )|,...,()'|,...,( 1

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

3. Stochastische Methode

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

khkhP

PkhP

M

1

0

exp),(

)(),(

P(„Rμ tritt auf“)P(„Keine Reaktion“)

3.2 Der Algorithmus3. Stochastische Methode

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

3.2 Der Algorithmus3. Stochastische Methode

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

3. Stochastische Methode

Generierung eines Zufallpaares (τ , μ) aus P(τ , μ)

„Konditionieren“:

→ τ ist stetig (0 ≤ τ < ∞) und μ ist diskret (μ = 1, 2, …, M) !

)|()(),( PPP

M

PP1

),()(

M

P

PP

1

),(

),()|(

a

aP )|()exp()( aaP

aμ = hμ kμ

M

aa1

Substitution

Generierung eines Zufallpaares (τ , μ) aus P(τ , μ)

3. Stochastische Methode

)exp()( aaP a

aP )|(

)exp(1)( aF

1

1ln

1

ra

12

1

1

aara

]1,0[1 Ur ]1,0[2 Ur

1

)|(a

aF

3.2 Der Algorithmus3. Stochastische Methode

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

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

3. Stochastische Methode

3.2 Der Algorithmus3. Stochastische Methode

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

3.3 Simulation mit COPASI

3. Stochastische Methode

COPASI … COmplex PAthway SImulator

} 0 ,2 ,2 { 321 kkk YYYXXXR

Y(t)

X(t)

3. Stochastische Methode

DM (Direct Method)

… direkte Generierung von μ und τ mit Inversionsmethode → benötigt 2 Zufallszahlen pro Iteration

FDM (First Reaction Method)

1. Generiere für jede Reaktion τμ

2. Durchlaufe Reaktion mit niedrigstem τμ

→ benötigt M Zufallszahlen pro Iteration

4. Effizienzverbesserung

4.1 Verbesserungen im Überblick4.2 LDM4.3 Tau-Leaping

4.1 Verbesserungen im Überblick

1. NRM (Next Reaction Method):Erstellung eines dependency graph

Verwendung einer indizierten priority queue Wiederverwendung von Zufallszahlen

2. ODM (Optimized Direct Method):Umordnung der Reaktionen nach Häufigkeit

Bedarf einer Prä-Simulation

3. SDM (Sorting Direct Method):Dynamisches ordnen der Reaktionen

4. LDM (Logarithmic Direct Method)

4. Effizienzverbesserung

4. Effizienzverbesserung

→ gibt an, welche aμ verändert werden müssen (geg. Rμ)

4.1 Verbesserungen im Überblick

1. NRM (Next Reaction Method):Erstellung eines dependency graph

Verwendung einer indizierten priority queue Wiederverwendung von Zufallszahlen

2. ODM (Optimized Direct Method):Umordnung der Reaktionen nach Häufigkeit

Bedarf einer Prä-Simulation

3. SDM (Sorting Direct Method):Dynamisches ordnen der Reaktionen

4. LDM (Logarithmic Direct Method)

4. Effizienzverbesserung

4. Effizienzverbesserung

→ Effiziente Verwaltung der absoluten Reaktionszeiten

Indexstruktur

Label

Index

Zeit

Zeiger

4.1 Verbesserungen im Überblick

1. NRM (Next Reaction Method):Erstellung eines dependency graph

Verwendung einer indizierten priority queue Wiederverwendung von Zufallszahlen

2. ODM (Optimized Direct Method):Umordnung der Reaktionen nach Häufigkeit

Bedarf einer Prä-Simulation

3. SDM (Sorting Direct Method):Dynamisches ordnen der Reaktionen

4. LDM (Logarithmic Direct Method)

4. Effizienzverbesserung

4. Effizienzverbesserung

Kostenintensivster Schritt der Methode ist das Auffinden der nächsten ablaufenden Reaktion im System

12

1

1

aara

M

aa1

→

→ Propensities werden beinahe zweimal aufsummiert

→ Suchtiefe im zweiten Schritt liegt in O(M)

1.

2.

4.2 LDM (Logarithmic Direct Method)

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)

4.3 Tau-Leaping Approximation4. Effizienzverbesserung

→ Ü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:

)(~rmit )( aPoissonrD

M

iii rtXtX

1

)()(

5. Anwendung

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

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-SimulationDichtefunktion der ReaktionswahrscheinlichkeitGenerierung zweier Zufallszahlen

Verschiedene Möglichkeiten der Effizienzverbesserung

Simulation komplexer Netzwerke

Zusammenfassung

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 Time Evolution of Coupled Chemical 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

Literatur

ENDE