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