Modelovaní dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinis1: PB051: Výpočetní metody v bioinformatice a systémové biologii David Šafránek 27.4.2011 Tento projekt je spolufinancován Evropským sociálním fondem a státním rozpočtem České republiky. EVROPSKÁ UNIE SACÍMI INVESTICE DO ROZVOJE VZDELÁVANÍ Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinis1: 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 Determinis1: 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 Determinis1: Energetický proces chemických reakcí n e kata lyžovaná reakce ©h /, Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinis1: Energetický proces chemických reakci • 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 Determinis1: Model dynamiky chemických reakci • 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 A/i molekul látky Si a A/2 molekul látky S2, pak náhodnou proměnnou % charakterizující pravděpodobnost kolize daného počtu molekul Si a S2 v časovém okamžiku dt lze modelovat: X = (cdt)N1- N2 • předpokládáme S\ a S2 různé látky • c je stochastická konstanta charakterizující průměrnou frekvenci úspěšných kolizí molekul Si a S2 za jednotku času • dt uvažováno limitně Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinism Model dynamiky chemických reakci • závislost frekvence kolizí c na absolutní teplotě T (Arheniův zákon): c oc e rt • Ea ... aktivační energie reakce • R ... plynová konstanta rovnovážný poměr frekvence dopředně a zpětné reakce u reversi bi In ich reakcí: Keq = e rt • AG ... 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 Determinis1: Stochastický model reakční dynamiky • uvažujme systém n substancí S = {S\, ...,Sn} provázaných m reakcemi R = {R\,Rm} • uvažujeme pouze reakce 1. a 2. řádu • systém zapisujeme pomocí stoichiometrické matice M rozměru n x m: —K, je-li K ■ S; reaktantem Rj K", je-li K • Sj produktem Rj • závislé reakce: dep(Ri, Rj) <^ 3k. Mki • MkJ < 0 Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Stochastický model reakční dynamiky • počet molekul substance S; v čase ř budeme značit A/,(ř) • náhodnou proměnnou rozložení počtů molekul substancí v čase ř charakterizujeme vektorem: X(t) = (N1(t),...,Nn(t)) • vývoj tohoto rozložení X(ř) v čase charakterizujeme jako stochastický proces: {X(ř)|řGR+} Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Motivace pro spojitý Markovův řetězec • stochastický proces {X(ŕ)|ŕ e M+} • 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{Ľ > t + t\U > t} = Pľ{U > t} • tuto vlastnost má exponenciálně distribuovaná proměnná • W - Exp(A) • y ■■■ průměrná čekací doba Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinis1: Exponenciálni rozloženi X ~ Exp(A) pokud: '\e~Xx, x>0, 0, jinak. Pro distribuční funkci dostáváme: 6c(x) Střední hodnota: 0, x < 0, 1 - e~Xx, x > 0. E(X) = i Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinis1: Exponenciálni rozloženi Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinism 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) = (A/i,Nn) je doba do provedení lib. reakce R; G R charakterizována rozložením Exp(xi{X,Cj)) R 0^ * R S j —> * Xi(X, c,) = (c, • dt) ■ Nj R Sp+ S q —> * Xi(X, c,) = (c, • dt) • A/p • A/q R 2Sj * • doba do nejbližší reakce má rozložení Exp(x{X,c)), kde m x{X,c) = J2xí{X,q), c = (ci,...,cm) ;=i • pravděpodobnost provedení reakce /?,•: P(/?;) = x'(x'^ Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinism Stochastický model reakční dynamiky Ri Xi = (ci • t/ŕ) • A/i R2 C c2 J2 -> X2 = (c2 • t/ŕ) • A/2 R3 Si + S2 X3 = (c3 • t/ŕ) • A/i • A/2 Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinis1: Monte Cario simulace Gillespiho přímá metoda 1. inicializace X(0) 2. výpočet ci) V/ G {1,m} v aktuálním stavu X 3. výpočet x(X, c) = J2T=i Xi{X, q) ^. simulace doby t do následující události - sampluj t G Exp(\{X, c)) 5. t := t + t č>. výběr reakce /?,• s pravděpodobností ^x'c) 7. X(t) :=X7 + M(j) 8. pokud t < Tmax, iteruj (2) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinis1: 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 nxxp://magnet.sysxemsbiology.net/software/Dizzy/ Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinis1: 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 Determinism Interakce při expresi genů - prokaryota Sigma protein recognizes promoter Core RNA polymerase DNA complementary strand Watson / DNA coding strand Promoter CTGjTTGACA|ATTAATCATCGAACTAG[rATAATÄGTACGCA -35 box -10 box mRNA start Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Interakce při expresi genů - eukaryota Activators Basal factors Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Interakce při expresi genů 1. vytvoření slabé vazby TF-DNA (nespecifické regiony DNA) 2. lineární pohyb TF po DNA (ID 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 Determinis1: 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 Determinis1: Interakce při expresi genů • modelování na úrovni buňky • ID 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 • ID pohyb po DNA je stochastický proces Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinism Interakce při expresi genů • pravděpodobnost ID pohybu lze charakterizovat exponenciálním rozložením vzhledem k rozdílu volných energií původní a cílové pozice: _ . . AG P(m) oc e rt- • AG = G' — Gq ... rozdíl volné Gibsovy energie iniciální (Go) a cílové pozice (G1) • T ... absolutní teplota [K] • R ... molární plynová konstanta [J ■ K^1 ■ mol-1] e RT Dr~3'\ _ e RT RT +e f P{mzach AGf' AGj?, v 3 ' AG|' AGJ, 1+e +e TíTTíT ag| ag|, 1+e +e Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Modelováni formace transkripčniho 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 Modelováni formace transkripčního komplexu definujeme přechody mezi stavy: DNA + DNAA DNAa -H DNA + A DNA + B -H DNAB DNAb -H DNA + B DNAB + A k^B DNAAB DNAAb k^B DNAB + A DNAA + B k^B DNAab DNAab 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 Determinis1: Modelováni 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 a ~ DNA + DNA • KA-A + DNA • KB • B + DNA • KB • KB2AB ■ B ■ A MeKB = ^KA = ^KB2AB = ^ • 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 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 B TA RNAo RNAo RNAl steps : M RNA! RNA0ffßp + BTA RNAoffBP m RNA, steps : N m RNA -^í CmRNA kdegmRNA CmRNA -> ribks'cmA'rib AAo AAo ^ÍAAusteps : P AAi clr + Hb clr ^5 prot. steps : Q ^degprot prot —> Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinism 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 Determinis1: 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.ž. We V.£„evp«vy» = l • pro každý stav v £ V je přiřazen parametr čekací doby av G M+ • 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 Determinis1: Spojitý Markovův řetězec Příklad Modelováni' dynamiky Transkripční' regulace Varianty Gillespiho algoritmu Aproximativní' metody Determinism Simulace spojitého Markovova řetězce II inicializace počáteční distribuce X(0) t := 0 u := X(0) while ( true ) wait_time := Exp (a,,) for each s, t x} = Pr{X1 > x n X2 > x n ... n Xn > x} n n = []Pr{X > x} = Yl e"Ax = e-x^'=iA/ i=l i=l Pro parametr minimálního rozložení platí: PT{Xk = min{X1,...,Xn}} Xk Ai + • • • + Xn Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody 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(Xaw) Ae měřící časový úsek s distribucí £xp(A3E) • jakmile některé doběhnou, přesun do příslušného stavu Pľ{min{Aw,AE} = Aw} - Aai4/+Aa£ PT{min{Aw,AE}=AE} = J£fr- Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinis1: Spojitý Markovův řetězec Převod na tradiční zápis - příklad • AE ~ £xp(15), Aw ~ Exp(5) . Pi{min{Aw,AE} = AE} = ^ = .75 . Pr{min{Aw,AE} = Aw} = = .25 Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinis1: Gillespiho přímá metoda (direct) 1. inicializace X(0) 2. výpočet Xi{X, ci) € {1; m} v aktuálním stavu X 3. výpočet x(X,c) = E™iX;(X,c(-) ^. simulace doby t do následující události - sampluj t G Exp(\{X, c)) 5. t := t + t fj. výběr reakce R; s pravděpodobností ^[x'^ 7. X(t) -=XT + M(j) 8. pokud t < Tmax, iteruj (2) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Varianta Gillespi: Metoda nejbližši reakce 1. inicializace t := 0, X(t), c = (ci,cm) 2. inicializace succs-times = 0 3. :/■: |1.....m\: • výpočet X/(X, Cj) v aktuálním stavu X{t) • pro reakci R-, sam pluj t,- G Exp(xi{X, c,)) • succs-times := succs-times U {t,-} ^. výpočet j, i,- = 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 Determinism 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 / simulace m(/) náhodných čísel, kde m(/) je počet následníků / 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 Determinis1: Inverzní rozloženi Uvažujme proměnnou U s rovnoměrným (uniformním) rozložením U ~ 1/(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 0 má proměnná X = —j \n(U) rozložení X ~ £xp(A). Rozepíšeme-li distribuční a kumulativní funkci f(x), F(x): f(x) = Xe-Xx F(x) = 1 - e-Xx F(x) je invertibilní: Triviálně platí (1 — U) ~ 1/(0,1) a tedy lze položit x = —j ln(u) F-\u) = -j\n(l-u) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Skálování exponenciálního rozložení Uvažujme proměnnou X o rozložení X ~ Exp(A). Proměnná Y = aX má rozložení Y ~ Exp(^). FY{y) = P{Y t Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinis1: Gibson-Bruck vs. Přimá metoda • G-B: pouze jedna simulace času v každém kroku • klíčovou vlastností G-B je selektivní update x,-(X, q) - pouze u relevantních (závislých) reakcí • selektivní přímý výpočet X/(X, q) a x(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 Determinis1: 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 Determinism 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: e-Xt(\t)k Pr{X = k}=Ě kde A je odpovídající parametr Poissonova rozložení. Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinis1: Poissonovo rozloženi X ~ Po(A) Modelováni' dynamiky Transkripční' regulace Varianty Gillespiho algoritmu Aproximativní' metody Determinism 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)|ř G M^}. Uvažujme časové intervaly délky t. 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 + T) -X{t) ~ Po(At) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinism Souvislost Poissonova a exponenciálního rozložení Uvažujme Poissonův proces {X(ř)|ř > 0} t.ž. X ~ Po(A). 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(Xt). FT(t) = Pr{7 < t} = 1 -Pr{7 > t} = 1 - Pr{A/t = 0} _ -i _ (e-At)(At)° — 1 0! = 1 - e-At Tedy platí: T - Exp(A) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinis1: 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 Ař Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Obecné schema aproximativního algoritmu 1. inicializace t := 0, X(ř), c = (ci,cm) 2. pro každé / e {1,m}\ • výpočet Xi{X,ci) v aktuálním stavu X{t) • simulace počtu reakcí R-, v intervalu At: • #/?; € Po(Xi(X, c;)At) 3. update X := X + M ■ #/?2, #/?m>7 ^. update ř := ř + Ař 5. pokud ř < Tmax iteruj (2) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Problémy aproximativního algoritmu • volba velikosti Ař • chceme rychlý, ale přitom dostatečně přesný výpočet • pro jednu simulaci nemusí být konstantní Ař vyhovující • předpokládáme konstantní frekvence reakcí po celé Ař • možnost definovat variabilní Ař v závislosti na aktuálním stavu X(ř) a škále jednotlivých reakčních konstant c Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Gillespiho T-leap metoda • aproximativní algortimus s variabilní Ař (r) • r 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 r • charakterizováno magnitudou změny frekvencí reakcí v průběhu r • r 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 Determinism Gillespiho T-leap metoda • post-leap (přechod od X(ŕ) k X' = X(ř + r)): • V; < m ověření velikosti • pokud velikosti nejsou dostatečně malé, přepočítej s menším t • problém: směřuje k malým změnám stavů Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinism Gillespiho r-leap metoda • pre-leap (odhad přechodu X(ŕ) k X' = X{t + r)): • 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^^x+X^c^t • parametr přesnosti e: V/ G {1,m}. \Xi(X', c,) - Xi(X, c,)| < e ■ X{X, c) • lze kombinovat s exaktní metodou (v kroku kde např. • nepřesnost vzniká při nelineární dynamice reakčních frekvencí Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinisi 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 Determinisi Deterministický model reakční dynamiky • uvažujme systém n substancí S = {S\, ...,Sn} provázaných m reakcemi R = {R\,Rm} • uvažujeme pouze reakce 1. a 2. řádu • systém zapisujeme pomocí stoichiometrické matice M rozměru n x m: —K, je-li K ■ S; reaktantem Rj KJe-li K • Sj produktem Rj Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Deterministický model reakční dynamiky • uvažovány vysoké molární koncentrace látek v buňce • koncentraci substance S; v čase ř budeme značit [S/](í) • systém v čase ř charakterizujeme vektorem: X(ř) = ([S1](r),...,[Sn](r)) • vývoj X v čase: • průměrné chování lze charakterizovat exponenciální funkcí Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinisi Převod počtu molekul na molární koncentraci • molární koncentrace [M]: n m=v kde n je množství látky [mol], V je objem roztoku [/] • vyjadřuje se pomocí Avogadrovy konstanty (počet částic v 1 molu): N kde Na Avogadrova konstanta [mol 1], V objem roztoku [/] a N je počet molekul. • převodní faktor 7 = Na ■ V: N = m ■ 7 Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Deterministický model reakční dynamiky • 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 Determinisi Deterministický model reakční dynamiky = MM • jaká funkce má stejný tvar jako její derivace? • f(t) = 1 + t + t2/2! + t3/3! + t4/4! + ... f(t) = e> • platí oře* Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinisi Deterministický model reakční dynamiky d-W = k- \A](t) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Deterministický model reakční dynamiky _MM = k [A](t)^ m) = [A](0) e-k • lineární dif. rce 1. řádu • jednoznačné řešení • numericky aproximovatelné Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody 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 R-, promítnuta do okamžitého reakčního toku qr,(t) R 0^ * qRl(t) = k; ■ dt R S j —> * gR,(t) = (k-dt)-[Sj](t) R Sp + S q —> * gR,(t) = (k; ■ dt) ■ [Sp](t) ■ [S,](t) R 2Sj * QRl(t) = {k- dt) ■ [SjY Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinisi Stochastický vs. deterministický model Pro reakci /?,- G R definujeme převodní vztah mezi stoch. frekvencí Cj a det. kinetickou konstantou k,: typ reakce R\ c; kj Sj —> * Sp + Sq * 2Sj -»■ * k; = c j ki = c j • 7 K, — 2 Modelovaní dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determi Modelováni' dynamiky Transkripční' regulace Varianty Gillespiho algoritmu Aproximativní' metody Determinisi Eulerova metoda • aproximativní řešení y(ř) (Euler): y'{t) = f{t,y{t)) y(0) = yo • přesné řešení 0, tn = nAt: yn w ip(tn) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinisi Eulerova metoda Pro exaktní řešení ip(t) platí: = V(tn) + St1 f(tMt))dt Schéma numerické aproximace: yn+i = yn + o- kde ftn+l a& / f(t,(tn + ^) a t£>(tn + Ař) • kroky algoritmu: kn,i = f(tn,yn) kn,2 = f(tn + jf,yn + ^knA) knA = f{t„ + At,y„ + hkn>3) yn+i = yn + %[kn,i + 2knt2 + 2/c„)3 + knA\ Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinisi Porovnáni metod diferenciální rovnice: y'(t)=y-2t y(0) = 3 • exaktní řešení: y(t) = 2 + 2t + et • srovnání simulací: Elder Improved Euler Runge—Kut ta steps error #evals error #evals error #evals 5 2.3 x ÍO"1 5 1.6 x ÍO"2 10 3.1 x 10~s 20 50 2.7 x íir2 50 1.8 x 10-4 100 3.6 x ÍO"9 200 500 2.7 x ÍO"3 500 1.8 x 10"6 1000 3.6 x 10"13 2000 Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinisi Přesnost numerické simulace Uvažujme systém: y'(t) = f(t,y),y(0)=yo Označme exaktní řešení {tn), platí y„+i = (t„) + hf(tn,(t) pro (j>'(tn) = f(tn,{tn) + h4>'{tn) Uvažujme Taylorův rozvoj pro (t„+i) = (t„ + h): <}>(tn + h) = (tn) + '(tn)h + l-4>"{tn)h2 + ^4>'"(tn)h3 + ■■■ Pn+i = 4>{tn+i) -y„+i [4,(tn) + 4>'(tn)h+±rh2 + ±h3 + \4>"{tn)h2 + i + [4>(tn) + h'(tn)} Chyba aproximována (pro konstantu K a délku kroku h): Pn+1 « Kh2 + 0(h3) Modelováni' dynamiky Transkripční' regulace Varianty Gillespiho algoritmu Aproximativní' metody Determinisi Adaptivní metody • chceme, aby chyba v každém kroku nebyla větší než e • řešení pomocí prediktoru chyby • při n-tém kroku (výpočet yn+\ z y„): 1. použití dvou různých algoritmů (dva výsledky pro yn+i — A\ a A2) 2. aproximace lokální chyby p„+i « \A\ — A2\ 3. pokud ^ < e, pak yn+í = A2 4. jinak iteruj (1) pro h' < h Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinisi Adaptivní Eulerova metoda Uvažujme systém: y'(t) = f(t,y),y(0)=yo Označme {tn + h) + ^Kh2 + 0(h3) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinisi Adaptivní Eulerova metoda Ai-A2 = (tn + h) + Kh2 + 0(h3) - (tn + h)-\Kh2- 0(h3) = \Kh2 + 0{h3) Lze tedy aproximovat \Kh2 k> A\ — A2 a označit _ \Aľ-A2\ _ iK. p — h ~ 2r\n. Pokud p > e, postup opakujeme pro h' t.ž. \\K\h1 & ^h1 < e. Např. h' = .9-/7 P Pokud p < e, nastavíme yn+\ = 2A2 — A\. Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Další adaptivní metody • Fehlbergova metoda • v každém kroku 3 odhady: h = f{tn,yn) a _ w .irr.rl k = f{tn + h,yn + hf1) ^Zy"U\1 + h] h = f{tn + %,yn + %[f1 + f2]) A2 = yn + l[h + f2 + *h] • aproximace p ■ h update faktor h' = •9y'^/J • Kutta-Merson, LSODE, ... Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody Determinism 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 Determinism 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, ř := 0 2. vypočti X/(X, q) v čase ŕ 3. na základě výsledku nastav Ař det. simulace 4- vypočti trajektorie spojitých proměnných v intervalu [ř, ř +Ař] 5. na základě výsledku vypočti X/(X, q) a rozhodni provedení případné diskrétní události (reakce) 6. pokud žádná diskrétní událost, nastav ř := ř + Ař 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 := ti • updatuj spojité proměnné k časovému bodu t\ • updatuj diskrétní proměnné relevantní reakci Rj 8. pokud ř < Tmax, iteruj (2) Modelování dynamiky Transkripční regulace Varianty Gillespiho algoritmu Aproximativní metody 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, ř := 0 2. výpočet X/(X, q) rychlých reakcí a nastav r pro r-leaping 3. pro pomalé reakce předpokládej konst. x,(X, q), rozhodni provedení pomalé reakce v intervalu [ř, ř + Ař] 4- pokud žádná pomalá reakce k provedení, vypočti r-leap update na rychlých reakcí na t := t + r 5. pokud pomalá reakce k provedení • identifikuj tuto reakci Rj a čas tj pro její provedení • update t := tj • proved r-leap update rychlých reakcí na tj • proved update pomalých reakcí vzhledem k provedení Rj 6. pokud ř < Tmax, iteruj (2)