Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist PB051: Výpočetní metody v bioinformatice a systémové biologii David Šafránek 27.4.2010 Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Obsah Modelování dynamiky chemických reakcí Modelování dynamiky transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Deterministické metody Hybridní metody Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Obsah Modelování dynamiky chemických reakcí Modelování dynamiky transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Deterministické metody Hybridní metody Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Energetický proces chemických reakcí Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Energetický proces chemických reakcí • různé energetické stavy molekuly • např. komplex AB méně stabilní než individuální výskyt molekul A, B • při přechodu mezi energ. stavy dochází k výměně energie • energie požadována pro aktivaci procesu (aktivační energie) • energie uvolněna během procesu (volná energie) • pro biologický systém je zdrojem většiny energie metabolismus • absolutní teplota ovlivňuje kinetickou energii molekul • pro reakci (úspěšnou kolizi) musí byt splněno: • správná prostorová konfigurace (orientace) molekul • dostatek kinetické energie Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Model dynamiky chemických reakcí • fixujeme-li konstantní teplotu, objem a tlak, orientace a kinetická energie molekul je stochastickým jevem • uvažujeme dobře promíchané médium • lze definovat průměrnou pravděpodobnost kolize molekul v časovém okamžiku • určeno fyzikálními vlastnostmi reagujících molekul • máme-li N1 molekul látky S1 a N2 molekul látky S2, pak náhodnou proměnnou χ charakterizující pravděpodobnost kolize daného počtu molekul S1 a S2 v časovém okamžiku dt lze modelovat: χ = (c · dt) · N1 · N2 • předpokládáme S1 a S2 různé látky • c je stochastická konstanta charakterizující průměrnou frekvenci úspěšných kolizí molekul S1 a S2 za jednotku času • dt uvažováno limitně Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Model dynamiky chemických reakcí • závislost frekvence kolizí c na absolutní teplotě T (Arheniův zákon): c ∝ e −EA RT • EA ... aktivační energie reakce • R ... plynová konstanta • rovnovážný poměr frekvence dopředné a zpětné reakce u reversibilních reakcí: Keq = e− ∆G RT • ∆G ... volná energie vyměněná při reakci • např. komplex kooperujících TF má vyšší stabilitu Polach K. J., Widom J., A Model for the Cooperative Binding of Eukaryotic Regulatory Proteins to Nucleosomal Target Sites, Journal of Molecular Biology, Volume 258, Issue 5, 24 May 1996, Pages 800-812, ISSN 0022-2836, DOI: 10.1006/jmbi.1996.0288. Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Stochastický model reakční dynamiky • uvažujme systém n substancí S = {S1, ..., Sn} provázaných m reakcemi R = {R1, ..., Rm} • uvažujeme pouze reakce 1. a 2. řádu • systém zapisujeme pomocí stoichiometrické matice M rozměru n × m: Mij = −K, je-li K · Si reaktantem Rj K, je-li K · Si produktem Rj • závislé reakce: dep(Ri , Rj ) ⇔ ∃k. Mki · Mkj < 0 Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Stochastický model reakční dynamiky • počet molekul substance Si v čase t budeme značit Ni (t) • náhodnou proměnnou rozložení počtů molekul substancí v čase t charakterizujeme vektorem: X(t) = N1(t), ..., Nn(t) • vývoj tohoto rozložení X(t) v čase charakterizujeme jako stochastický proces: {X(t)|t ∈ R+ } Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Motivace pro spojitý Markovův řetězec • stochastický proces {X(t)|t ∈ R+} • spojitý čas pobytu ve stavu • lze zachytit rozložením W samplujícím „čekací dobu mezi změnami stavů • požadujeme markovskou vlastnost nezávislosti na historii: Pr{U > t + τ|U > τ} = Pr{U > t} • tuto vlastnost má exponenciálně distribuovaná proměnná • W ∼ Exp(λ) • 1 λ ... průměrná čekací doba Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Exponenciální rozložení X ∼ Exp(λ) pokud: fX (x) = λe−λx , x ≥ 0, 0, jinak. Pro distribuční funkci dostáváme: FX (x) = 0, x < 0, 1 − e−λx , x ≥ 0. Střední hodnota: E(X) = 1 λ Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Exponenciální rozložení Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Stochastický model reakční dynamiky • interleaving: při přechodu X(t) → X(t + dt) je updatována právě jedna složka X • provedení právě jedné reakce z R • provedení reakce uvažováno jako okamžitý jev (trvá nulový čas) • ve stavu X(t) = N1, ..., Nn je doba do provedení lib. reakce Ri ∈ R charakterizována rozložením Exp(χi (X, ci )) Ri ∅ → ∗ χi (X, ci ) = ci · dt Ri Sj → ∗ χi (X, ci ) = (ci · dt) · Nj Ri Sp + Sq → ∗ χi (X, ci ) = (ci · dt) · Np · Nq Ri 2Sj → ∗ χi (X, ci ) = (ci · dt) · Nj ·(Nj −1) 2 • doba do nejbližší reakce má rozložení Exp(χ(X, c)), kde χ(X, c) = m i=1 χi (X, ci ), c = c1, ..., cm • pravděpodobnost provedení reakce Ri : P(Ri ) = χi (X,ci ) χ(X,c) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Stochastický model reakční dynamiky N2 = 1 N1 = 0 N2 = 0 N1 = 1 N1 = 0 N2 = 0 N2 = 1 N1 = 1 c2 c3 c2c1 c1 R1 S1 c1 −→ χ1 = (c1 · dt) · N1 R2 S2 c2 −→ χ2 = (c2 · dt) · N2 R3 S1 + S2 c3 −→ χ3 = (c3 · dt) · N1 · N2 Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Monte Carlo simulace Gillespiho přímá metoda 1. inicializace X(0) 2. výpočet χi (X, ci ) ∀i ∈ {1, ..., m} v aktuálním stavu X 3. výpočet χ(X, c) ≡ m i=1 χi (X, ci ) 4. simulace doby τ do následující události – sampluj τ ∈ Exp(χ(X, c)) 5. t := t + τ 6. výběr reakce Ri s pravděpodobností χi (X,ci ) χ(X,c) 7. X(t) := XT + M(j) 8. pokud t < Tmax , iteruj (2) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Nástroj Dizzy • nástroj pro simulaci dynamiky sítí chemických reakcí • obsahuje stochastické i deterministické solvery • mimo přímý Gillespiho algoritmus zahrnuje další varianty stochastické simulace • podpora zobrazení modelů v Cytoscapu http://magnet.systemsbiology.net/software/Dizzy/ Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Obsah Modelování dynamiky chemických reakcí Modelování dynamiky transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Deterministické metody Hybridní metody Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Interakce při expresi genů – prokaryota Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Interakce při expresi genů – eukaryota Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Interakce při expresi genů 1. vytvoření slabé vazby TF-DNA (nespecifické regiony DNA) 2. lineární pohyb TF po DNA (1D difůze) 3. disociace vazby TF-DNA 4. prostorový pohyb – “skok” (3D difůze) 5. (1-4) iterováno dokud nenalezena regulační sekvence DNA 6. afinita TFBS a kooperativní interakce transkripčního komplexu snižují vliv difůze a stabilizují funkční vazbu TF-DNA Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Interakce při expresi genů • interference při paralelním “skenování” DNA • velké množství molekul téhož TF • molekuly různých TF • u eukaryot navíc komplikovaná lokální prostorová struktura DNA • na úrovni buňky nutno uvažovat stochasticitu (husté prostředí v okolí DNA) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Interakce při expresi genů • modelování na úrovni buňky • 1D difůze uvažována v diskrétních krocích ⇒ 1 krok ∼ 1 nukleotid za jednotku času • stupně volnosti: • pohyb ve směru 5 − 3 • pohyb ve směru 3 − 5 • zachování pozice • proteiny TF v buňce mají specifickou distribuci volné energie vůči vazbě s DNA • energie individuální molekuly fluktuuje vzhledem k náhodným kolizím s ostatními molekulami • fluktuace způsobuje dynamické změny afinity proteinů k DNA • 1D pohyb po DNA je stochastický proces Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Interakce při expresi genů • pravděpodobnost 1D pohybu lze charakterizovat exponenciálním rozložením vzhledem k rozdílu volných energií původní a cílové pozice: P(m) ∝ e− ∆G RT • ∆G = G − G0 ... rozdíl volné Gibsovy energie iniciální (G0) a cílové pozice (G ) • T ... absolutní teplota [K] • R ... molární plynová konstanta [J · K−1 · mol−1 ] P(m5 3 ) = e− ∆G5 3 RT 1+e− ∆G5 3 RT +e− ∆G3 5 RT P(m3 5 ) = e− ∆G3 5 RT 1+e− ∆G5 3 RT +e− ∆G3 5 RT P(mzach) = 1 1+e− ∆G5 3 RT +e− ∆G3 5 RT Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Modelování formace transkripčního komplexu • uvažujeme dva TF A a B kooperativně aktivující transkripci genu • zavedeme následující stavy regulatorního regionu DNA: • DNA ... volný regulatorní region DNA • DNAA ... navázán A, nikoliv B • DNAB ... navázán B, nikoliv A • DNAAB ... navázán A i B Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Modelování formace transkripčního komplexu • definujeme přechody mezi stavy: DNA DNA_A DNA_ABDNA_B DNA + A kfA −→ DNAA DNAA krA −→ DNA + A DNA + B kfB −→ DNAB DNAB krB −→ DNA + B DNAB + A kfB2AB −→ DNAAB DNAAB krB2AB −→ DNAB DNAA + B kfA2AB −→ DNAAB DNAAB krA2AB −→ DNAA + B • při kooperačním faktoru Kq lze uvažovat: kfA2AB · kfA = kfB2AB · kfB = Kq · kfA · kfB Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Modelování formace transkripčního komplexu • transkripce (tvorba mRNA) je proces o několi řádů pomalejší vzhledem k formaci transkripčního komplexu • lze jej uvažovat jako nekonečně rychlý (tj. stabilní) • okupaci promotoru pak charakterizujeme poměrem frekvence výskytu stavu DNAAB vzhledem k sumě frekvencí všech možných stavů: α = DNA · KB · KB2AB · B · A DNA + DNA · KA · A + DNA · KB · B + DNA · KB · KB2AB · B · A kde KB = kfB kbB , KA = kfA kbA , KB2AB = kfB2AB kbB2AB • ekvivalentně lze psát: α = DNA · KA · KA2AB · A · B DNA + DNA · KA · A + DNA · KB · B + DNA · KA · KA2AB · A · B Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Modelování průběhu transkripce • uvažujeme činnost komplexu RNA polymerázy • předpokládáme, že všechny faktory důležité pro transkripční komplex jsou k dispozici v libovolné míře BTA k·BTA·DNAAB −→ RNA0 RNA0 tr −→ RNA1, steps : M RNA1 tr −→ RNAoffBP + BTA RNAoffBP tr −→ mRNA, steps : N mRNA trR −→ cmRNA cmRNA kdegmRNA −→ rib ks ·cmRNA·rib −→ AA0 AA0 AAsr −→ AA1, steps : P AA1 AAsr −→ clr + rib clr AAsr −→ prot, steps : Q prot kdegprot −→ Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Obsah Modelování dynamiky chemických reakcí Modelování dynamiky transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Deterministické metody Hybridní metody Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Spojitý Markovův řetězec Spojitý Markovův řetězec lze definovat jako přechodový graf MC = V , E, p : • stavy V reprezentují prvky jevového pole • přechody E jsou ohodnoceny pravděpodobností p : E → (0, 1 t.ž. ∀v ∈ V . v ∈V p( v, v ) = 1 • pro každý stav v ∈ V je přiřazen parametr čekací doby αv ∈ R+ • indexujeme-li prvky V , pak přechodový graf lze zapsat maticí: Qij = p( vi , vj ) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Spojitý Markovův řetězec Příklad • průměrná čekací doba ve stavu Surf je 3 minuty, což je 3 60 = 1 20 hod ⇒ αS = 20 • αW = 7.5 • αE = 15 Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Simulace spojitého Markovova řetězce // inicializace počáteční distribuce X(0) t := 0 u := X(0) while ( true ) wait time := Exp(αu) for each s, t x} = Pr{X1 > x ∩ X2 > x ∩ ... ∩ Xn > x} = n i=1 Pr{Xi > x} = n i=1 e−λx = e−x Pn i=1 λi Pro parametr minimálního rozložení platí: Pr{Xk = min{X1, ..., Xn}} = λk λ1 + · · · + λn Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Spojitý Markovův řetězec Převod na tradiční zápis • rozhodovací procedura při opuštění stavu Surf • v okamžiku vstupu do stavu Surf se spustí stopky: AW měřící časový úsek s distribucí Exp(λaW ) AE měřící časový úsek s distribucí Exp(λaE ) • jakmile některé doběhnou, přesun do příslušného stavu • Pr{min{AW , AE } = AW } = λaW λaW +λaE • Pr{min{AW , AE } = AE } = λaE λaW +λaE Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Spojitý Markovův řetězec Převod na tradiční zápis – příklad • AE ∼ Exp(15), AW ∼ Exp(5) • Pr{min{AW , AE } = AE } = 15 15+5 = .75 • Pr{min{AW , AE } = AW } = 5 15+5 = .25 Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Gillespiho přímá metoda (direct) 1. inicializace X(0) 2. výpočet χi (X, ci ) ∀i ∈ {1, ..., m} v aktuálním stavu X 3. výpočet χ(X, c) ≡ m i=1 χi (X, ci ) 4. simulace doby τ do následující události – sampluj τ ∈ Exp(χ(X, c)) 5. t := t + τ 6. výběr reakce Ri s pravděpodobností χi (X,ci ) χ(X,c) 7. X(t) := XT + M(j) 8. pokud t < Tmax , iteruj (2) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Varianta Gillespi: Metoda nejbližší reakce 1. inicializace t := 0, X(t), c = (c1, ..., cm) 2. inicializace succs times = ∅ 3. ∀i ∈ {1, ..., m}: • výpočet χi (X, ci ) v aktuálním stavu X(t) • pro reakci Ri sampluj ti ∈ Exp(χi (X, ci )) • succs times := succs times ∪ {ti } 4. výpočet j, tj = min(succs times) 5. t := t + tj 6. X(t) := XT + M(j) 7. pokud t < Tmax , iteruj (2) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Přímá vs. metoda nejbližší reakce • obě metody exaktně simulují CTMC • přímá: v každém kroku simulace dvou náhodných čísel • mnr: v lib. stavu i simulace m(i) náhodných čísel, kde m(i) je počet následníků i v CTMC • přímá metoda efektivnější • mnr však poskytuje filosoficky odlišný přístup, který lze dále akcelerovat Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Inverzní rozložení Uvažujme proměnnou U s rovnoměrným (uniformním) rozložením U ∼ U(0, 1) a nechť F(·) libovolná invertibilní kumulativní distribuční funkce. Pak rozložení X = F−1 (U) je charakterizováno kumulativní distribuční funkcí F(·). P(X ≤ x) = P(F−1 (U) ≤ x) = P(U ≤ F(x)) = FU (F(x)) = F(x) Realizaci libovolného rozložení X charakterizovaného invertibilní kumulativní distribuční funkcí F(x) lze simulovat pomocí uniformního rozložení U. Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Souvislost uniformního a exponenciálního rozložení Uvažujme proměnnou U s rovnoměrným (uniformním) rozložením U ∼ U(0, 1). Pro lib. λ > 0 má proměnná X = − 1 λ ln(U) rozložení X ∼ Exp(λ). Rozepíšeme-li distribuční a kumulativní funkci f (x), F(x): f (x) = λe−λx F(x) = 1 − e−λx F(x) je invertibilní: F−1 (u) = − 1 λ ln(1 − u) Triviálně platí (1 − U) ∼ U(0, 1) a tedy lze položit x = − 1 λ ln(u). Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Škálování exponenciálního rozložení Uvažujme proměnnou X o rozložení X ∼ Exp(λ). Proměnná Y = αX má rozložení Y ∼ Exp(λ α ). FY (y) = P(Y ≤ y) = P(αX ≤ y) = P(X ≤ y α ) = FX (y α ) FX (x) = 1 − e−λx ⇒ FY (y) = FX (y α ) = 1 − e− λ α y ⇒ Y ∼ Exp(λ α ) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Algoritmus Gibson-Bruck 1. inicializace t := 0, X(t), c = (c1, ..., cm) 2. pro každé i ∈ {1, ..., m}: • výpočet χi (X, ci ) ve stavu X(0) • pro reakci Ri sampluj ti ∈ Exp(χi (X, ci )) 3. inicializace succs times = ∅ 4. pro každé i ∈ {1, ..., m}: • succs times := succs times ∪ {ti } 5. výpočet j, tj = min(succs times) 6. t := tj 7. X(t) := XT + M(j) 8. update χj (X, cj ) v novém stavu X(t) 9. sampluj tj := t + Exp(χj (X, cj )) 10. pro každé i ∈ {1, ..., m} splňující dep(Ri , Rj ): • ulož χlast i := χi (X, ci ) • update χi (X, ci ) • ti := t + χlast i (X,ci ) χi (X,ci ) (ti − t) 11. pokud t < Tmax , iteruj (3) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Algoritmus Gibson-Bruck • relativní čas (do nejbližší události) nahrazen časem absolutním (čas nejbližší události) • není nutno generovat nové časy pro všechny reakce • ale pouze pro reakce závislé na provedené reakci • umožněno díky markovovské vlastnosti • nové časy závislých reakcí se nesamplují, ale vypočítavají z předchozích • předchozí časy závislých reakcí jsou použity pro výpočet nových časů • škálování podmíněné požadavkem ti > t Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Gibson-Bruck vs. Přímá metoda • G-B: pouze jedna simulace času v každém kroku • klíčovou vlastností G-B je selektivní update χi (X, ci ) – pouze u relevantních (závislých) reakcí • selektivní přímý výpočet χi (X, ci ) a χ(X, c) může zrychlit i přímou metodu • rychlost obou metod výrazně závisí na implementaci • v G-B je obecně více operací • implementaci výpočtu minimálního reakčního času lze realizovat pomocí indexované prioritní fronty Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Obsah Modelování dynamiky chemických reakcí Modelování dynamiky transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Deterministické metody Hybridní metody Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Poissonův proces Uvažujme experiment představující četnost výskytu diskrétních událostí v časovém intervalu t. Pro zachycení tohoto měření zavedeme náhodnou proměnnou X. Pak platí: X ∼ Po(λt) Jinými slovy, pravděpodobnostní funkce určující pravděpodobnost jevu, že za dobu t nastane právě k událostí, má následující tvar: Pr{X = k} = e−λt (λt)k k! kde λ je odpovídající parametr Poissonova rozložení. Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Poissonovo rozložení X ∼ Po(λ) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Poissonův proces Nyní uvažujme náhodnou proměnnou X z předchozího slidu v závislosti na čase, tzv. stochastický proces určený množinou náhodných proměnných {X(t)|t ∈ R+ 0 }. Uvažujme časové intervaly délky τ. Předpokládejme následující podmínky: • výskyty událostí ve dvou libovolných vzájemně disjunktních časových intervalech jsou nezávislé jevy • pravděpodobnostní rozložení četnosti výskytu událostí v daném časovém intervalu závisí pouze na délce intervalu • žádné dvě události nemohou nastat současně („interleaving ) Pak pro libovolný interval (t, t + τ platí: X(t + τ) − X(t) ∼ Po(λτ) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Souvislost Poissonova a exponenciálního rozložení Uvažujme Poissonův proces {X(t)|t ≥ 0} t.ž. X ∼ Po(λ). Zavedeme náhodnou proměnnou T zachycující dobu do první nejbližší události (od počátečního okamžiku). Pro t > 0 uvažujme náhodnou proměnnou Nt zachycující počet událostí v intervalu (0, t . Z definice platí Nt ∼ Po(λt). FT (t) = Pr{T ≤ t} = 1 − Pr{T > t} = 1 − Pr{Nt = 0} = 1 − (e−λt )(λt)0 0! = 1 − e−λt Tedy platí: T ∼ Exp(λ) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Motivace pro aproximativní algoritmy • exaktní metody simulují spojitý markovův řetězec • lze aproximovat diskretizací času • rozdělení časové osy na diskrétní intervaly • počet událostí v lib. časovém intervalu charakterizován Po(λ) • lze tedy samplovat počet provedení reakcí daného typu v daném časovém intervalu ∆t Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Obecné schema aproximativního algoritmu 1. inicializace t := 0, X(t), c = (c1, ..., cm) 2. pro každé i ∈ {1, ..., m}: • výpočet χi (X, ci ) v aktuálním stavu X(t) • simulace počtu reakcí Ri v intervalu ∆t: • #Ri ∈ Po(χi (X, ci )∆t) 3. update X := X + M · #R1, #R2, ..., #Rm T 4. update t := t + ∆t 5. pokud t < Tmax iteruj (2) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Problémy aproximativního algoritmu • volba velikosti ∆t • chceme rychlý, ale přitom dostatečně přesný výpočet • pro jednu simulaci nemusí být konstantní ∆t vyhovující • předpokládáme konstantní frekvence reakcí po celé ∆t • možnost definovat variabilní ∆t v závislosti na aktuálním stavu X(t) a škále jednotlivých reakčních konstant c Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Gillespiho τ-leap metoda • aproximativní algortimus s variabilní ∆t (τ) • τ uvžováno co největší při garanci požadované přesnosti • přesnost určena mírou tolerance uvažování konstantní frekvence reakcí v rámci intervalu τ • charakterizováno magnitudou změny frekvencí reakcí v průběhu τ • τ voleno tak, aby míra změny všech reakčních frekvencí byla minimalizována Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Gillespiho τ-leap metoda • post-leap (přechod od X(t) k X = X(t + τ)): • ∀i ≤ m ověření velikosti |χi (X , ci ) − χi (X, ci )| • pokud velikosti nejsou dostatečně malé, přepočítej s menším τ • problém: směřuje k malým změnám stavů Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Gillespiho τ-leap metoda • pre-leap (odhad přechodu X(t) k X = X(t + τ)): • výpočet očekávaného nového stavu E(X ) (v t := t + τ) E(X ) = X + E(r)M kde E(r) je vektor očekávaných frekvencí reakcí: E(r)i = χi (X, ci )τ • parametr přesnosti : ∀i ∈ {1, ..., m}. |χi (X , ci ) − χi (X, ci )| ≤ · χ(X, c) • lze kombinovat s exaktní metodou (v kroku kde např. τ < 2 χ(X,c) • nepřesnost vzniká při nelineární dynamice reakčních frekvencí Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Obsah Modelování dynamiky chemických reakcí Modelování dynamiky transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Deterministické metody Hybridní metody Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Deterministický model reakční dynamiky • uvažujme systém n substancí S = {S1, ..., Sn} provázaných m reakcemi R = {R1, ..., Rm} • uvažujeme pouze reakce 1. a 2. řádu • systém zapisujeme pomocí stoichiometrické matice M rozměru n × m: Mij = −K, je-li K · Si reaktantem Rj K, je-li K · Si produktem Rj Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Deterministický model reakční dynamiky • uvažovány vysoké molární koncentrace látek v buňce • koncentraci substance Si v čase t budeme značit [Si ](t) • systém v čase t charakterizujeme vektorem: X(t) = [S1](t), ..., [Sn](t) • vývoj X v čase: dX dt = f (X) • průměrné chování lze charakterizovat exponenciální funkcí Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Převod počtu molekul na molární koncentraci • molární koncentrace [M]: m = n V kde n je množství látky [mol], V je objem roztoku [l] • vyjadřuje se pomocí Avogadrovy konstanty (počet částic v 1 molu): m = N NA · V kde NA Avogadrova konstanta [mol−1], V objem roztoku [l] a N je počet molekul. • převodní faktor γ = NA · V : N = m · γ Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Deterministický model reakční dynamiky A k −→ • předpokládejme nádobu jednotkového objemu obsahující v čase t látku A v molárním množství [A] [mol] • kolik množství látky A „odteče za jednotku času? • hodnota přímo úměrná hodnotě [A] v daném okamžiku − d[A](t) dt = k · [A](t) • koeficient úměrnosti je konstanta k [s−1 ] tzv. reakční konstanta (koeficient) - determinuje rychlost reakce rozpadu (“odtoku”) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Deterministický model reakční dynamiky [A](t) dt = k · [A](t) • jaká funkce má stejný tvar jako její derivace? • f (t) = 1 + t + t2 /2! + t3 /3! + t4 /4! + ... f (t) = et • platí det dt = et Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Deterministický model reakční dynamiky A k −→ −d[A](t) dt = k · [A](t) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Deterministický model reakční dynamiky A k −→ −d[A](t) dt = k · [A](t)⇔ [A](t) = [A](0) · e−kt • lineární dif. rce 1. řádu • jednoznačné řešení • numericky aproximovatelné Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Deterministický model reakční dynamiky • spojité chování: při přechodu X(t) → X(t + dt) jsou updatovány všechny složky X (souběžný spojitý tok reakcí) • časová informace o běhu reakce Ri promítnuta do okamžitého reakčního toku Ri (t) Ri ∅ → ∗ Ri (t) = ki · dt Ri Sj → ∗ Ri (t) = (ki · dt) · [Sj ](t) Ri Sp + Sq → ∗ Ri (t) = (ki · dt) · [Sp](t) · [Sq](t) Ri 2Sj → ∗ Ri (t) = (ki · dt) · [Sj ]2 Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Stochastický vs. deterministický model Pro reakci Ri ∈ R definujeme převodní vztah mezi stoch. frekvencí ci a det. kinetickou konstantou ki : typ reakce Ri ci → ki Sj → ∗ ki = ci Sp + Sq → ∗ ki = ci · γ 2Sj → ∗ ki = ci ·γ 2 Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Eulerova metoda Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Eulerova metoda • aproximativní řešení y(t) (Euler): y (t) = f (t, y(t)) y(0) = y0 • přesné řešení ϕ(t): ϕ (t) = f (t, ϕ(t)) ϕ(0) = y0 • pro lib. n ≥ 0, tn = n∆t: yn ≈ ϕ(tn) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Eulerova metoda Pro exaktní řešení ϕ(t) platí: ϕ(tn+1) = ϕ(tn) + tn+1 tn ϕ (t)dt = ϕ(tn) + tn+1 tn f (t, ϕ(t))dt Schema numerické aproximace: yn+1 = yn + σ kde σ ≈ tn+1 tn f (t, ϕ(t))dt Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Eulerova metoda I dy dt = f (t, y) yn+1 = yn + ∆t · f (tn, yn) 1. init t0, y0, ∆t, n; 2. for j from 1 to n do 2.1 m := f (t0, y0); 2.2 y1 := y0 + ∆tm; 2.3 t1 := t0 + ∆t; 2.4 t0 := t1; 2.5 y0 := y1; 3. end Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Eulerova metoda II • aproximace σ (obsah pod křivkou) lichoběžníkem o obsahu: ∆t 2 · [f (tn, ϕ(tn)) + f (tn+1, ϕ(tn+1))] • pro ϕ(tn+1) nutno předpočítat aproximaci (známe již yn ≈ ϕ(tn)): ϕ(tn+1) ≈ ϕ(tn) + ϕ (tn)∆t ≈ yn + f (tn, yn)∆t • aproximace σ ≈ 1 2 [f (tn, yn) + f (tn+1, yn + f (tn, yn)∆t)]∆t • celkem dostáváme: y(tn+1) ≈ yn+1 = yn + 1 2 [f (tn, yn) + f (tn+1, yn + f (tn, yn)∆t)]∆t Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Runge-Kutta • σ upřesněno Simpsonovým pravidlem pro aproximaci obsahu pod parabolou: σ = R tn+1 tn ≈ ∆t 6 [f (tn, ϕ(tn)) + 4f (tn + ∆t 2 , ϕ(tn + ∆t 2 )) + f (tn + ∆t, ϕ(tn + ∆t))] • nutno aproximovat ϕ(tn + ∆t 2 ) a ϕ(tn + ∆t) • kroky algoritmu: kn,1 = f (tn, yn) kn,2 = f (tn + ∆t 2 , yn + h 2 kn,1) kn,3 = f (tn + ∆t 2 , yn + h 2 kn,2) kn,4 = f (tn + ∆t, yn + hkn,3) yn+1 = yn + h 6 [kn,1 + 2kn,2 + 2kn,3 + kn,4] Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Porovnání metod • diferenciální rovnice: y (t) = y − 2t y(0) = 3 • exaktní řešení: y(t) = 2 + 2t + et • srovnání simulací: Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Přesnost numerické simulace Uvažujme systém: y (t) = f (t, y), y(0) = y0 Označme exaktní řešení φ(t): φ (t) = f (t, φ(t)) Uvažujme Eulerovu metodu I s krokem ∆t = h: yn+1 = yn + hf (tn, yn) (Lokální) chyba numerické simulace: ρn+1 = φ(tn+1) − yn+1 Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Přesnost numerické simulace Předpokládáme-li yn = φ(tn), platí yn+1 = φ(tn) + hf (tn, φ(tn)) a tedy z exaktnosti řešení φ(t) pro φ (tn) = f (tn, φ(tn)) dostáváme: yn+1 = φ(tn) + hφ (tn) Uvažujme Taylorův rozvoj pro φ(tn+1) = φ(tn + h): φ(tn + h) = φ(tn) + φ (tn)h + 1 2 φ (tn)h2 + 1 3! φ (tn)h3 + · · · ρn+1 = φ(tn+1) − yn+1 = [φ(tn) + φ (tn)h + 1 2 φ h2 + 1 3! h3 + · · · ] − [φ(tn) + hφ (tn)] = 1 2 φ (tn)h2 + 1 3! + · · · Chyba aproximována (pro konstantu K a délku kroku h): ρn+1 ≈ Kh2 + O(h3 ) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Adaptivní metody • chceme, aby chyba v každém kroku nebyla větší než • řešení pomocí prediktoru chyby • při n-tém kroku (výpočet yn+1 z yn): 1. použití dvou různých algoritmů (dva výsledky pro yn+1 — A1 a A2) 2. aproximace lokální chyby ρn+1 ≈ |A1 − A2| 3. pokud ρn+1 h < , pak yn+1 = A2 4. jinak iteruj (1) pro h < h Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Adaptivní Eulerova metoda Uvažujme systém: y (t) = f (t, y), y(0) = y0 Označme φ(t) exaktní řešení y = f (t, y) t.ž. φ(tn) = yn. Pro A1 uvažujme Eulerovu metodu I: A1 = yn + hf (tn, yn) A1 = φ(tn + h) + Kh2 + O(h3 ) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Adaptivní Eulerova metoda Pro A2 uvažujme Eulerovu metodu 2-step (provedení dvou kroků, každý s ∆t = h 2 ): A2 = yn + h 2 f (tn, yn) + h 2 f (tn + h 2 , yn + h 2 f (tn, yn)) Označme ymid = yn + h 2 f (tn, yn) a A2 = ymid + h 2 f (tn + h 2 , ymid ). ρmid = K( h 2 )2 + O(h3 ) ρ2nd = K( h 2 )2 + O(h3 ) Konstanta K shodná v obou případech (neplatí pro člen O(h3)). Lze tedy psát (dle Taylorovy věty): A2 = φ(tn + h) + 1 2 Kh2 + O(h3 ) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Adaptivní Eulerova metoda A1 − A2 = φ(tn + h) + Kh2 + O(h3) − φ(tn + h) − 1 2 Kh2 − O(h3) = 1 2 Kh2 + O(h3) Lze tedy aproximovat 1 2 Kh2 ≈ A1 − A2 a označit ρ = |A1−A2| h ≈ 1 2 Kh. Pokud ρ > , postup opakujeme pro h t.ž. 1 2 |K|h ≈ ρ h h < . Např. h = .9 ρ h Pokud ρ < , nastavíme yn+1 = 2A2 − A1. Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Další adaptivní metody • Fehlbergova metoda • v každém kroku 3 odhady: f1 = f (tn, yn) f2 = f (tn + h, yn + hf1) f3 = f (tn + h 2 , yn + h 4 [f1 + f2]) A1 = yn + h 2 [f1 + f2] A2 = yn + h 6 [f1 + f2 + 4f3] • aproximace ρ = |A1−A2| h ≈ Kh2 • update faktor h = .9 ρ h • Kutta-Merson, LSODE, . . . Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Obsah Modelování dynamiky chemických reakcí Modelování dynamiky transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Deterministické metody Hybridní metody Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Kombinace deterministické a stochastické simulace Předpokládáme model s kombinací spojitých a diskrétních proměnných. 1. inicializuj systém, t := 0 2. vypočti χi (X, ci ) v čase t 3. na základě výsledku nastav ∆t det. simulace 4. vypočti trajektorie spojitých proměnných v intervalu [t, t + ∆t] 5. na základě výsledku vypočti χi (X, ci ) a rozhodni provedení případné diskrétní události (reakce) 6. pokud žádná diskrétní událost, nastav t := t + ∆t a updatuj pro tento bod spojité proměnné 7. pokud diskrétní událost • najdi (nejbližší) reakci Rj a čas tj jejího provedení • t := t1 • updatuj spojité proměnné k časovému bodu t1 • updatuj diskrétní proměnné relevantní reakci Rj 8. pokud t < Tmax , iteruj (2) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinist Multi-step metody stochastické simulace Metoda Puchalka and Kierzek (2004) (STOCKS2). Předpokládá rozlišení “rychlých” a “pomalých” reakcí. 1. inicializuj systém, t := 0 2. výpočet χi (X, ci ) rychlých reakcí a nastav τ pro τ-leaping 3. pro pomalé reakce předpokládej konst. χi (X, ci ), rozhodni provedení pomalé reakce v intervalu [t, t + ∆t] 4. pokud žádná pomalá reakce k provedení, vypočti τ-leap update na rychlých reakcí na t := t + τ 5. pokud pomalá reakce k provedení • identifikuj tuto reakci Rj a čas tj pro její provedení • update t := tj • proveď τ-leap update rychlých reakcí na tj • proveď update pomalých reakcí vzhledem k provedení Rj 6. pokud t < Tmax , iteruj (2)