V12 Stochastische Effekte 8. Juli 2011 Folien von Dr. Tihamer Geyer

Übersicht Stochastische Effekte => was sind das? => warum wichtig? => wie behandeln?

Gillespie-Algorithmus direkte Implementierung

Softwarewerkzeuge SS 11 – V12

2

Stochastische Effekte

Dichte =

Softwarewerkzeuge SS 11 – V12

ununterscheidbare Teilchen Volumen

3

Dichtefluktuationen N = 10 6

0

N = 100

1.67

6

4

4

2

2

0

0

N = 1000 6

3

6

0

9

0

6

4

4

2

2

0

0

3

Softwarewerkzeuge SS 11 – V12

6

0

N = 10000

55.6

9

0

0

0

5.56

3

6

0

9

556

3

6

9 4

Poisson-Verteilung Betrachte Kontinuum w mit im Mittel g Ereignissen pro Einheitsintervall Δw Annahmen: i) Seltenheit:

) √ g(x) = exp − w/ A0 = 0.13, wA = 0.45 = 0.26, wB = 0.55 = 0.61, wC = 0.45 Softwarewerkzeuge SS 11 – V12

100

rel. Häufigkeiten

1000 Simulationsläufe, Werte sichern bei t = 7.

0 0.00 A 300

0.50 B

A0 = 1001.00

C

150 0 0.00 400 A

0.50 B

1.00

A0 = 30

C

200 0 0.00

0.50

1.00

11

Gillespie-Algorithmus

naive Integration => shortcomings

2340

Exact Stochastic Simulation of Coupled Chemical Reactions

Gillespie => Idee an der Tafel Daniel => Parallelen zu denT. Gillesple Poisson-Annahmen

Danlel T. Gillespie Research Department, NaVal Weapons Center, China Lake, California 93555 (Received May 72, 1977) Publication costs assisted by the Naval Weapons Center

There are two formalisms for mathematically describing the time behavior of a spatially homogeneous chemical system: The deterministic approach regards the time evolution as a continuous, wholly predictable process which is governed by a set of coupled, ordinary differential equations (the “reaction-rateequations”);the stochastic approach regards the time evolution as a kind of random-walk process which is governed by a single differential-difference equation (the “master equation”). Fairly simple kinetic theory arguments show that the stochastic formulation of chemical kinetics has a firmer physical basis than the deterministic formulation, but unfortunately the stochastic master equation is often mathematically intractable. There is, however, a way to make exact numerical calculations within the framework of the stochastic formulation without having to deal with the master equation directly. It is a relatively simple digital computer algorithm which uses a rigorously derived Monte Carlo procedure to numerically simulate the time evolution of the given chemical system. Like the master equation, this “stochastic simulation algorithm” correctly accounts for the inherent fluctuations and correlations that are necessarily ignored in the deterministic formulation. In addition, unlike most procedures for numerically solving the deterministic reaction-rate equations, this algorithm never approximates infinitesimal time increments dt by finite time steps At. The feasibility and utility of the simulation algorithm are demonstrated by applying it to several well-known model chemical systems, including the Lotka model, the Brusselator, and the Oregonator.

D. Gillespie, J. Phys. Chem. 81 (1977) 2340–2361 Softwarewerkzeuge SS 11 – V12

12

Das Problem Zerfallsreaktion:

A => Ø

Wahrscheinlichkeit für ein Ereignis in (t, t+Δt) bei A(t) Molekülen = A(t) k Δt Naiver Algorithmus:

Es geht, aber:

A = A0 For every timestep: get random number r ε [0, 1) if r ≤ A*k*dt: A = A-1

A(t)*k*dt meistens passiert nichts

Gillespie überlegte sich: "Wann ist das nächste Ereignis?" Softwarewerkzeuge SS 11 – V12

13

Exponentiell verteilte Wartezeiten Idee: wann findet die nächste Reaktion statt? Mit A(t) Molekülen: g(A(t), s) = Wahrscheinlichkeit, dass keine Reaktion in (t, t+s) stattfindet + Poisson-Annahmen => DGL für Verteilung der Wartezeiten: g(A(t), s) = exp[−A(t) k s] = e−s/s0

mittlere Wartezeit:

Softwarewerkzeuge SS 11 – V12

s0 =

1 k A(t)

14

Gillespie-Algorithmus: 1 Spezies g(A(t), s) = exp[−A(t) k s]

Mit r = exp[−A(t) k s]

r ∈ [0,1]

Exponentielle Verteilung der Wartezeiten:

! " ! " 1 1 1 1 s = ln = ln k A(t) r α0 r

1

0 0

life time t

=> Gillespie-Algorithmus für Reaktion A => Ø: A = A0 While(A > 0): get random number r ∈ [0, 1) t = t + s(r) A = A-1 Softwarewerkzeuge SS 11 – V12

keine Abfrage!!! einfache Sprung in der Zeit 15

Gillespie- vs. direkte Integration Naiv:

Gillespie:

"Was ist die Wahrscheinlichkeit, dass im nächsten Δt ein Ereignis stattfindet?"

"Wie lange dauert es bis zum nächsten Ereignis?"

=> kleine feste Zeitschritte

=> variable Zeitschritte

=> Näherung erster Ordnung

=> exakt

30

30

• Gillespie

• Gillespie

• naive - analytic

20

• naive - analytic

NA

NA

25

10

20 0

0

1

Softwarewerkzeuge SS 11 – V12

2

t

3

4

0.00

0.20

t

0.40 16

Gillespie – viele Spezies Für beliebige Anzahl an möglichen Reaktionen: (i) Wahrscheinlichkeiten für jede einzelne Reaktion: αi i = 1, …, N gesamte W‘keit α0 = Σ αi ! " 1 1 s = ln (ii) Wartezeit bis zum nächsten beliebigen Ereignis: α0 r1 j−1

(iii) Wähle eine der Reaktionen aus gemäß:

∑ αi

i=1

α1 0

α2 α 1 α 1+α2

α3 α 1+α2+α3

α4

j

≤ α0 r2 < ∑ αi i=1

α5 α 1+…+α4

α6 α 1+…+α5

α0

(iv) Update der Teilchenzahlen Softwarewerkzeuge SS 11 – V12

17

Ein Beispiel mit zwei Spezies Reaktionen: k1

k2

A + A => Ø

A + B => Ø

k3

Ø => A

k4

Ø => B

Ratengleichungen: dA = k3 − 2 A2 k1 − A B k2 dt

Ass =

Gleichgewichtszustand: mit

k1 = 10–3 s–1

!

k3 − k4 2k1

k2 = 10–2 s–1

=> Ass = 10, Softwarewerkzeuge SS 11 – V12

dB = k4 − A B k2 dt

Bss = k3 = 1.2 s–1

k4 k2 A k4 = 1 s–1

Bss = 10 18

(a)

(b) k1

A + A −→ ∅

25A

solution of ODEs

20

number of B molecules

number of A molecules

25

k3

15

∅ −→ A

k

2 + B −→ ∅

solution of ODEs

(2.27)

20

k

4 ∅ −→ B.

Gillespie-Algorithmus

(2.28)

15

Erban, Chapman, Maini, arXiv:0704.1908v2 [q-bio.SC]

Let us note that we are not interested in chemical species C and D. Hence, we 10 10 notation of unimportant chemical replaced them by ∅, consistent with our previous species. To simulate the system of chemical reactions (2.27)–(2.28), we perform the 5 5 = n0 , B(0) = m0 at time t = 0): following four steps at time t (starting with A(0) 0(a4) 0

Generate two random numbers r1 , r2 uniformly distributed in (0, 1). 0 20 40 60 80 100 0 20 40 60 80 100 (b4) Computetime the[sec] propensity functions of each reaction by α A(t)(A(t)− 1)k1 , time [sec] 1 = α2 = A(t)B(t)k2 , α3 = k3 and α4 = k4 . Compute α0 = α1 + α2 + α3 + α4 . Fig. 2.3. Five realizations of SSA Number of molecules of chemical A (c4) Compute the time when the(a4)–(d4). next chemical reaction takes place as t+species τ where

(left panel) and B (right panel) are plotted as functions of time as solid lines. Different colours $ correspond to different realizations. The solution of#(2.33)–(2.34) is given by the dashed line. We 1 1 −3 −1 −2 −1 −1 −1 . use A(0) = 0, B(0) = 0, k1 = 10 sec τ, k =2 = 10ln sec ., k3 = 1.2 sec and k4 = 1 sec (2.29)

α0

r1

(d4) Compute the number of molecules at time t + τ by  A(t) − 2    A(t) − 1 A(t + τ ) = A(t) + 1    A(t)

if 0 ≤ r2 < α1 /α0 ; if α1 /α0 ≤ r2 < (α1 + α2 )/α0 ; if (α1 + α2 )/α0 ≤ r2 < (α1 + α2 + α3 )/α0 ; if (α1 + α2 + α3 )/α0 ≤ r2 < 1;

(2.30)

 B(t)    B(t) − 1 B(t + τ ) =  B(t)   B(t) + 1

if 0 ≤ r2 < α1 /α0 ; if α1 /α0 ≤ r2 < (α1 + α2 )/α0 ; if (α1 + α2 )/α0 ≤ r2 < (α1 + α2 + α3 )/α0 ; if (α1 + α2 + α3 )/α0 ≤ r2 < 1;

(2.31)

Then continue with step (a4) for time t + τ.

SSA (a4)–(d4) is a direct generalisation of SSA (a3)–(d3). At each time step, we first ask the question when will the next reaction occur? The answer is given by formula (2.29) which can be justified using the same arguments as formulae (2.5) or Softwarewerkzeuge SS 11 – V12 (2.10). Then we ask the question which reaction takes place. The probability that the i-th chemical reaction occurs is given by αi /α0 . The decision which reaction takes place is given in step (d4) with the help of the second uniformly distributed random

19

Stochastisch vs. Deterministisch k1

A + A => Ø

k2

12

k4

Ø => A

Ø => B

RADEK ERBAN ET AL.

(a)

(b)

25 20 15 10 5 0 0

25

solution of ODEs number of B molecules

number of A molecules

Erban, Chapman, Maini, arXiv:0704.1908v2 [q-bio.SC]

k3

A + B => Ø

20

40 60 time [sec]

80

100

solution of ODEs

20 15 10 5 0 0

20

40 60 time [sec]

80

100

Fig. 2.3. Five realizations of SSA (a4)–(d4). Number of molecules of chemical species A (left panel) and B (right panel) are plotted as functions of time as solid lines. Different colours correspond to different realizations. The solution of (2.33)–(2.34) is given by the dashed line. We use A(0) = 0, B(0) = 0, k1 = 10−3 sec−1 , k2 = 10−2 sec−1 , k3 = 1.2 sec−1 and k4 = 1 sec−1 .

the number of molecules at time t + τ by Softwarewerkzeuge(d4) SS 11 Compute – V12  A(t) − 2    A(t) − 1

if 0 ≤ r2 < α1 /α0 ; if α1 /α0 ≤ r2 < (α1 + α2 )/α0 ;

20

Stationärer Zustand k1

k2

k3

k4

A + A => Ø

A + B => Ø

Ø => A

Ø => B

k1 = 10–3 s–1

k2 = 10–2 s–1

k3 = 1.2 s–1

k4 = 1 s–1

kontinuierliches Modell: Ass = 10, Bss = 10 14

RADEK ERBAN ET AL.

(a)

(b) !3

x 10

30

8

25

0.08 6

20 15

4

10 2

stationary distribution

number of B molecules

Erban, Chapman, Maini, arXiv:0704.1908v2

Aus langen Gillespie-Simulationen: Ass = 9.6, Bss = 12.2



0.06

0.04

0.02

5 5

10 15 20 25 number of A molecules

30

0

0

2

4

6

8 10 12 14 16 18 20 22 24 number of A molecules

Fig. 2.4. (a) Stationary distribution φ(n, m) obtained by long time simulation of (a4)–(d4) for k1 = 10−3 sec−1 , k2 = 10−2 sec−1 , k3 = 1.2 sec−1 and k4 = 1 sec−1 . (b) Stationary distribution of A obtained by (2.35).

Softwarewerkzeuge SS 11 – V12 consider that we have a system of q chemical reactions. Let αi (t) be the propensity function of the i-th reaction, i = 1, 2, . . . , q, at time t, that is, αi (t) dt is the probability that the i-th reaction occurs in the time interval [t, t + dt). Then the Gillespie SSA

21

Photosynthese ist… …die Umwandlung von Lichtenergie in chemische Energie (einer der wichtigsten Prozesse weltweit) biologischer Überblick Stochastische SImulation aus einzelnen Reaktionen => Unterschiede zu Gillespie? => Reihenfolge von Näherung und Simulation Vesiweb

P P P A

Softwarewerkzeuge SS 11 – V12

22

Photosynthese in Rb. sphaeroides ATPase

ATP ADP + Pi 0

!

H+

1

3

LH2

0

!

1

E

6 4

NRC x RRC

3

e–H+

2

e–

H+

H+

e–-H+-Paare

5

Softwarewerkzeuge SS 11 – V12

bc1

e– 5

elektronische Anregungen NLH x RLH

Q

2 c2

Lichtenergie

H+

4

2

RC/LH1

7

H+

Protonengradient

Nbc1 x Rbc1 4

chemische Energie

NATPase x RATPase

"pH

H+

6

7

ATP

ADP + P 23

Chromatophoren-Vesikel 4 nm

7

0 6 H+

3.3 nm 37

nm

H+

4

12

nm

5

2

H+

45 nm

=> einfach: 4 Proteine + 2 Transporter + H+; Kristallstrukturen u. Reaktionen => klein:

abgeschlossen: definierte Randbedingungen für die Simulation => praktisch:

Anregung mit Licht, spektroskopische Messungen Softwarewerkzeuge SS 11 – V12

Geyer, Helms, BPJ 91 (2006) 921 + 927

24

n residues in facilitating this reaction is a subject of intense st. From H'- the electron transfers to a quinone, labeled about 100 ps to form P'+HQA'-; the electron then moves out 100 ps to a second quinone, labeled QB, to form QAQgC.*The electron-transfer kinetics at room temperature mmarized in Figure 1, along with the back electron transfer determined by using RCs with subsequent electron acceptors removed or reduced. The forward rates are several orders gnitude faster than theRC: recombination ensuring a high Photonrates, => Ladungstrennung um efficiency of charge separation. ause the chromophores involved in these reactions are l, closed-shell molecules in their ground state, the interte states have two characteristics. First, charge is separated considerable distance, so the states have large electric dipole nts; second, a radical pair is formed. The large electric dipole nt makes it possible to manipulate the relative energies of ntermediates and consequently the reaction dynamics with al electric fields. The presence of unpaired electrons in the l pair allows the reaction dynamics to be altered with exmagnetic fields. Because of the wide range of experimental eoretical tools that have been used to study RCs, much is n about their structure and function. Therefore, the elecransfer reactions in the photosynthetic RC offer an excepopportunity to test and refine theories of electron transfer pecific applications to long-distance electron transfer in ns, to develop new methods to study electron transfer, and blish which parameters are important for efficient artificial synthesis. 1AIJ.pdb he following, we review some details of the properties of We then focus on three experimental issues that have been imary concerns of our laboratory: the nature of the lowest y singlet electronic excited state of P, the identity of the initial on acceptor, and the energetics of the initial charge separation Obviously, there Softwarewerkzeuge are other important experimental and SS 11 – V12 tical issues which have occupied the attention of other tories but which will not be discussed in detail; up-to-date

F

Modellierung der Proteine Boxer et al, JCP 93 (1989) 8280

_.'

I

\

Figure 3. Reaction schemes without (A) and with (B) 13 for energetics).

QA

(see Figure

H+out

summaries can be found in recent monographs.Is Q A3 Properties of the Reaction Center. The A2 R C is the smallest isolable unit which performsRCthe initial charge separation steps bsQ of photosynthesis.16 EBy far the best studied RCs are those isolated A1 QH2sphaeroides D2 Rhodobacter from two species of purple bacteria: P and Rhodopseudomonas uiridis. The arrangement of the reactive components of the RC from Rps.D4uiridis is shown in Figure 1;' bsc

(1 5) The Photosynthetic Bacterial Center-Structure and DyA4 Reaction D4 c2ox c2red namics; Breton, J., Vermeglio, A,, Eds.; Plenum: New York, 1988. Perspectiues in Photosynthesis; Jortner, J., Pullman, B., Eds.; Elsevier Press: Amsterdam, 1989. 25 (16) Feher, G.;Okamura, M. Y . In The Photosynthetic Bacteria; Clayton, R. K.,Sistrom, W. R., Eds.; Plenum Press: New York, 1978; pp 349-386.

Stochastische Reaktionen Wenn BS frei ist => Assoziation möglich: BS + X => BS:X 1) sind alle Bedingungen erfüllt? 2) chemische Reaktionskinetik: d[BS:X] = kon [BS] [X] dt

Reaktionsrate:

Bindungs-W.keit pro BS:

Pon = kon [X] ∆t

H+out A3

A2

RC bsQ

E A1

D2

P D4 bsc c2red

Q

A4

D4

QH2

for each timestep Δt: for each reaction: conditions fulfilled? determine probability: perform reaction

c2ox

Softwarewerkzeuge SS 11 – V12

Geyer, Lauck, Helms, J Biotech 129 (2007) 212

26

Bsp: Elektronentransfer im RC // R1: transfers an electron to the Quinone // using the energy from an exciton { if (bs_Q && (e_P == 1) && (e_Q == 0) && ((He_Q == 0) || (He_Q == 1))) { if (LHPoolp->take_out(LH_kon)) { e_P = 0; e_Q = 1;

Bedingungen? Wahrscheinlichkeit? Reaktion!

writeInternals();

H+out

} }

A3

}

A2

RC bsQ

E A1

D2

P

Softwarewerkzeuge SS 11 – V12

QH2

D4

Protein = {BS; Reaktionen(Zustand)} Geyer, Lauck, Helms, J Biotech 129 (2007) 212

Q

bsc c2red

A4

D4

c2ox

27

"Pools-and-Proteins"–Modell H+ outside

Light

ADP + P

QH2 LHC

E

RC

Q

bc1

ΔΦ

ATPase

ATP

titratable groups

c2red c2ox

40 aktive Proteine • unabhängig voneinander • stochastische Reaktionen mit je 1 Molekül • Anzahl wie auf dem Vesikel

H+ inside

19 passive Pools • ein Pool pro Metabolit • hier: Diffusion ist schnell

Verbindungen definieren das biologische System => Pfade als "emergent behavior" Softwarewerkzeuge SS 11 – V12

Geyer, Lauck, Helms, J Biotech 129 (2007) 212

28

Web Interface Simulationen über Konfigurationsdatei oder web-interface @ service.bioinformatik.uni-saarland.de/vesiweb

• Verstehen der Prozesse • Modell-Verifikation + Parametrisierung gegen Experimente Softwarewerkzeuge SS 11 – V12

[Florian Lauck. T.G., 2006]

29

Stochastische Effekte Oxidationszustand des Cytochrom c-Pools bei kontinuierlicher Beleuchtung H+

Light

Q 20 (LHC/RC)

5 (bc1)2

1 ATPase

QH2 10 c2 ox H+ 10 c2 red

4.4 W/m2 const.

Softwarewerkzeuge SS 11 – V12

30

Steady State Fluktuationen 60 Sek. bei konstanter Beleuchtung mit 10 RC/LHC-Dimeren und 4 bc1-Dimeren => Oxidationszustand des Cytochrom c-Pools 20

3 W/m2

8 W/m2



15 10 5 0

0

3

6

I [W/m2]

9

12

5.5 W/m2

4.5 W/m2

6.5 W/m2

=> weicher Übergang mit starken Fluktuationen Softwarewerkzeuge SS 11 – V12

31

Deterministisch vs. Stochastisch 20

20 t=1s t=3s

10

15

# c2ox

15

t=0.3s

5

5 0

10

0 0

3

6

9

I [W/m2]

Gleichungen mitteln, dann simulieren

0

3

6

I [W/m2]

9

12

Simulieren, dann Ergebnisse mitteln

scharfer Übergang



weicher Übergang auch für lange Zeiten

nur numerische Unsicherheiten



Fluktuationen ≈ Signal

Reproduzierbare Werte



nur Mittelwert reproduzierbar

Softwarewerkzeuge SS 11 – V12

32

Vesimulus: Warum nicht Gillespie? Gillespie:

Pools-and-Proteins:

hohe Anzahl von Spezies/Zuständen => Komplexität!

Belegung unabhängiger Bindungsstellen => Zustände als Filter

Jeweils eine Reaktion

Mehrere Reaktionen pro Zeitschritt

Numerisch effizient, exakt

Natürliche Modellierung

Viele unabhängige Proteine Mittelung

Gemittelte Proteine Simulation

gemittelte Belegungszahlen Softwarewerkzeuge SS 11 – V12

Simulation

?=

Belegungszahlen pro Protein Mittelung

gemittelte Belegungszahlen 33

Zusammenfassung Stochastische Effekte sind oft wichtig — wann? Modellierung I: Gillespie • numerisch effizient • numerisch exakt • biochemische Prozesse => stochastische Integration kontinuierlicher Ratengleichungen Modellierung II: Pools-and-Proteins • direkte Modellierung individueller Proteine • Netzwerktopologie bleibt erhalten (auf allen Ebenen) => mikroskopisches Verhalten jedes einzelnen Proteins

Softwarewerkzeuge SS 11 – V12

34