PA054: Formální modely v systémové biologii David Šafránek 30.4.2010 Obsah Nelineární systém v rovině • f (x(t)) = f1(x(t)), f2(x(t)) T • předpokládáme ∀i. fi (x(t)) obecná funkce spojitá a diferencovatelná na X • soustava ODEs: d dt x1 x2 = f1(x1, x2) f2(x1, x2) Linearizace Mějme systém: dx1 dt = f1(x1, x2) dx2 dt = f2(x1, x2) Nechť ¯x1, ¯x2 t.ž. f1(¯x1, ¯x2) = 0 ∧ f2(¯x1, ¯x2) = 0. Uvažujme funkce z1(t), z2(t) s obory hodnot v malém okolí 0 t.ž. x1(t) = ¯x1 + z1(t) x2(t) = ¯x2 + z2(t) Linearizace Označme z(t) = z1(t), z2(t) T . Uvažujme Taylorův rozvoj funkce f blízkém okolí ¯x = ¯x1, ¯x2 : f (¯x + z(t)) ˙= f (¯x) + Jf (¯x) · z(t) Kde Jf (¯x) je Jacobiho matice funkce f vyhodnocená v bodě ¯x1, ¯x2 . Jf (¯x1, ¯x2) = ∂f1 ∂x1 ∂f1 ∂x2 ∂f2 ∂x1 ∂f2 ∂x2 x1,x2 = ¯x1,¯x2 Linearizace Platí: dx1 dt = d(¯x1+z1(t)) dt = dz1 dt dx2 dt = d(¯x2+z2(t)) dt = dz2 dt Připomeňme f (¯x) = 0, z čehož dostáváme: f (z(t)) = Jf (¯x) · z(t) d dt z1 z2 = Jf (¯x) · z1 z2 Máme tedy lineární systém charakterizující dynamiku v malém okolí stabilního bodu ¯x. Linearizace Řekneme, že bod ¯x1, ¯x2 je hyperbolické ekvilibrium pokud f (¯x1, ¯x2) = 0 a navíc všechna vlastní čísla Jf (¯x1, ¯x2) mají nenulovou reálnou složku. Hartman-Grobmanova věta: Nechť ¯x1, ¯x2 hyperbolické ekvilibrium. Pak v blízkém okolí ¯x1, ¯x2 je vektorové pole funkce x(t) topologicky ekvivalentní vektorovému poli funkce z(t). Analýza stability obecných spojitých modelů Uvažujme soustavu: dx dt = f (x) kde: • x = x1, ..., xn T • f (x) = f1(x), ..., fn(x) T • n > 2 Analýza stability obecných spojitých modelů Pro daný stabilní bod ¯x (f (¯x) = 0) charakterizujeme typ stability následujícím postupem: 1. určíme všechna vlastní čísla matice Jf (¯x) ⇒ λ1, ..., λn 2. definujeme následující množiny: • L+ = {λi | λi > 0} • L− = {λi | λi < 0} • L0 = {λi | λi = 0} • Lc+ = {λi | ∃λj . λi = ¯λj ∧ Re(λi ) > 0} • Lc− = {λi | ∃λj . λi = ¯λj ∧ Re(λi ) < 0} • Lc0 = {λi | ∃λj . λi = ¯λj ∧ Re(λi ) = 0} 3. analyzujeme vlastní podprostory Jakobiho matice: • vl. čísla L+ ∪ Lc+ určují tzv. nestabilní manifold • vl. čísla L− ∪ Lc− určují tzv. stabilní manifold • podprostor charakterizovaný L0 ∪ Lc0 určují tzv. centrální manifold, chování v něm nelze linearizací analyzovat Obsah Bayramovův model – reakční síť Z r4 r5 r3 YX r2 r1 2 Y + X k1 Y + Z 2Y Z k2 X Y k4 k5 k3 Bayramovův model – simulace Bayramovův model – stoichiometrická analýza Formálně uvažujeme model M = {X, Y , Z}, {r1, r2, r3, r4, r5}, recnet, ∅, map, rates , kde: map : M =   0 −1 0 1 0 1 −1 −1 0 1 −1 1 0 0 0   ∀i ∈ {1, ..., n}. rates(ri ) = ki recnet = { Y + Z, 2Y , Y + X, X , 0, X , 0, Y , Y , 0 } • konzervační vztahy mezi koncentracemi substrátů neexistují: h(M) = 3 • třída komplexnosti systému κ: complexes(M) = 6, λM = 2 κ = 6 − 2 − 3 = 1 Bayramovův model – stoichiometrická analýza • elementární módy: M · v = 0 : v = 1, 1, 0, 1, 0 ∨ v = 0, 0, 1, 0, 0, 1 r5 r3 YX r2 r1 2 r4 Z Z r4 r5 r3 YX r2 r1 2 Bayramovův model – kinetika dX dt = k4 − k2yx dY dt = k5 − k3y − k2yx + k1yz dZ dt = k2yx − k1yz Bayramovův model – kinetika dX dt = k4 − k2yx dY dt = k5 − k3y − k2yx + k1yz dZ dt = k2yx − k1yz k1 = 300 k2 = 250 k3 = 0.1 k4 = 0.0005 k5 = 0.0001 Bayramovův model – analýza stability • stabilní stav ¯x ˙= 0.002, 0.001, 0.00167 • určení typu stability: Jf (¯x) ˙=   −0.25 −0.5 0 −0.25 −0.1 0.3 0.25 0 −0.3   • vlastní čísla: λ1 ˙= 0.0048 − 0.11i λ2 ˙= 0.0048 + 0.11i λ3 ˙= − 0.66 • nestabilní manifold: spirálová rostoucí oscilace (λ1, λ2) • stabilní manifold: λ3 (exponenciální rozpad) • ¯x je nestabilní ekvilibrium Bayramovův model – fázový diagram (simulace) Bayramovův model – fázový diagram v okolí ekvilibria Obsah Model regulační sítě nutrienty enzymy regulatory produkty metabolismu proteiny signaly zpetne vazbyzpetne vazby Model regulační sítě A1 r1 A2 A3 A4 A5 r2 r3 r4 r5 2 3 Model regulační sítě A1 r1 A2 A3 A4 A5 r2 r3 r4 r5 2 3 Model regulační sítě • uvažujeme modely tvaru M = S, R, reanet, regnet, map, rates , kde regnet = ∅ • regnet ⊆ S × R × {inh, act} • každou regulaci lze vyjádřit reakční sítí • důvod pro zavádění regulací je abstrakce ⇒ snížení dimenzionality ⇒ zjednodušení modelu • regulace lze kvantizovat parametry • omezující podmínky specifikující aktivní konfiguraci regulátoru • např. prahová koncentrace regulátoru • příkladem je regulace proteosyntézy (transkripční regulace) Distribuce časových škál v buňce Experimentálně zjištěný parametr E.Coli Vazba molekuly signálu na transkripční faktor vedoucí ke změně aktivity faktoru ∼ 1msec Vazby aktivního faktoru na operon DNA ∼ 1sec Transkripce + translace jednoho genu ∼ 5min Životnost mRNA ∼ 2 − 5min 50% změna koncentrace stabilního proteinu ∼ 1h Transkripce v prokaryotické buňce Transkripce v eukaryotické buňce Jak modelovat dynamiku transkripční regulace? • chceme modelovat změny koncentrace proteinů se zaměřením na jejich syntézu a rozpad v buňce • časová škála životnosti proteinů • např. v E. coli ∼ desítky minut - hodiny • vs. doba jednoho buněčného cyklu (u E. coli min. 30 minut) • doba odezvy syntézy jedné molekuly proteinu ∼ jednotky minut ⇒ dynamika proteinů zajímavá v časové škále desítek minut Transkripční kinetika dle mass action r1 Y r2 • uvažujme dynamiku (syntéza,rozpad) proteinu Y syntéza: → Y rozpad: Y → • obě reakce zahrnují mnoho dílčích jevů r1 tvorba mRNA, činnost transkripčních faktorů, tvorba tRNA, činnost RNA polymerázy, . . . r2 rozpad proteinu v buňce, snížení koncentrace vlivem růstu, . . . Transkripční kinetika dle mass action • uvažujme dynamiku (syntéza,rozpad) proteinu Y syntéza: β −→ Y β ... produkční koeficient [M−1 · s−1] rozpad: Y γ −→ γ ... degradační koeficient [s−1] Transkripční kinetika dle mass action • uvažujme dynamiku (syntéza,rozpad) proteinu Y syntéza: β −→ Y β ... produkční koeficient [M−1 · sec−1] rozpad: Y γ −→ γ ... degradační koeficient [sec−1] • dle mass action kinetics dostáváme: dY dt = β − γY Transkripční kinetika dle mass action • dle mass action kinetics máme: dY dt = β − γY • β ∼ β · p β ... rychlost transkripce (inflow mRNA) [#mRNA · sec−1] p ... množství molekul Y vyrobených translací z jedné mRNA • γ ... rychlost rozpadu proteinu + roztahování buňky Negativní regulace (represe) transkripce • transkripční faktor X degraduje transkripci proteinu Y • řízeno snížením řídícího signálu S (uvolnění X z vazby X : S) r1 Y X r2 S XS r3 r4 Negativní regulace (represe) transkripce • transkripční faktor X degraduje transkripci proteinu Y • řízeno snížením řídícího signálu S (uvolnění X z vazby X : S) Y Y gen Y X dY dt = β − γYr1 Y r2 X Negativní regulace (represe) transkripce • transkripční faktor X degraduje transkripci proteinu Y • řízeno snížením řídícího signálu S (uvolnění X z vazby X : S) Y Y gen Y X*X signal Sx dY dt = β − γYr1 Y r2 X Negativní regulace (represe) transkripce • transkripční faktor X degraduje transkripci proteinu Y • řízeno snížením řídícího signálu S (uvolnění X z vazby X : S) gen Y X*X signal Sx X* dY dt = −γYr1 Y r2 X Vazba represoru na promotor gen Y DNA P X • reakční síť: XPX + P • totální koncentrace promotoru (PT ) konzervována mezi volné promotory (volné molekuly substrátu P) a vázané promotory (komplexy [XP]) PT = [XP] + P Dynamika vazby represoru na promotor • dynamika komplexu XP (mass action): d[XP] dt = konXP − koff [XP] • kin. konstanta mass action (zde kon) je obecně limitována fyzikálními možnostmi kolize proteinových molekul: kon ∼ 108 − 109 M−1 sec−1 • uvažujeme-li kolize molekuly proteinu X s DNA, je limit vyšší (fyzikální tendence “1-dimenzionálního pohybu” po DNA): kon ∼ 1010 − 1011 M−1 sec−1 • koff [sec−1] determinováno vlastnostmi chemické vazby X : P Dynamika vazby represoru na promotor • reakční toky v síti XPX + P se ustálí: d[XP] dt = 0 ⇔ konXP − koff [XP] = 0 ⇔ koff kon = XP [XP] • Kd = koff kon [M] nazýváme disociačním prahem • čím větší je Kd tím slabší je vazba X : P Dynamika vazby represoru na promotor • připomeňme PT = [XP] + P • v buňce je mnoho kopií DNA, tedy mnoho promotorů • pro procento všech volných (aktivních) promotorů platí: P PT = 1 1 + X Kd • většina vazeb represor-promotor splňuje koff > 1sec−1 • časová škála [XP] ∼ [sec] (vs. transkripce ∼ [min]) • průměr přes mnoho událostí syntézy a rozpadu [XP] v čase ⇒ P PT ≈ pravděpodobnost jevu, že P je v daném okamžiku volný P PT = 1 ⇔ X = 0 P PT = 1 2 ⇔ X = Kd Aktivita promotoru • v nepřítomnosti represoru (P volný) se váže RNA-polymeráza • tato situace určuje maximální transkripční koeficient βmax βmax ∼ 10−4 − 100 [#mRNA · sec−1 ] • βmax je dána mnoha biofyzikálními aspekty • např. pozice vazebného místa [DNA]:[RNApolymeráza] • aktivita promotoru P regulovaného represorem X: α− P (X) = βmax 1 + X Kd • tzv. Hillova funkce Aktivita promotoru v závislosti na regulujícím represoru [X] Kd βmax tc α− (X) Exprese genu v závislosti na regulujícím represoru X Y P gen y r1 Y r2 X reakční síť: r1 : → Y ; X r2 : Y → model kinetiky: dY dt = α− P (X) − γY Exprese genu v závislosti na regulujícím represoru X Y P gen y K r1 Y r2 X reakční síť: r1 : → Y ; X− (K) r2 : Y → model kinetiky: dY dt = βmax 1 1 + X K − γY Aktivní vs. pasivní forma represoru • represor X je účinný pouze ve volné (aktivní) formě (X∗) • vazba se signální molekulou S • neaktivní forma represoru je zcela v komplexu [XS] Signal S X Y P gen y S:X Aktivní vs. pasivní forma represoru Signal S X Y P gen y S:X r1 Y X r2 S XS r3 r4 reakční síť: r3, r4 : X + S ↔ [XS] r1 : → Y ; X− (Kd ) r2 : Y → Aktivní vs. pasivní forma represoru Signal S X Y P gen y S:X r1 Y X r2 S XS r3 r4 reakční síť: r3, r4 : X + S ↔ [XS] r1 : → Y ; X− (Kd ) r2 : Y → rates(r1), rates(r2) << rates(r3), rates(r4) Aktivace represoru • reakční síť: X + S XS • totální koncentrace represoru v buňce (XT ) konzervována mezi X (volné výskyty) a [XS] (vázané výskyty) XT = [XS] + X • dynamika komplexu [XS] (mass action): d[XS] dt = konXS − koff [XS] Aktivace represoru • ve stabilním stavu platí: koff kon [XS] = XS • Ks = koff kon [M] je tzv. disociační konstanta [XS] = XT S S + Ks • koncentrace aktivního represoru: X = XT − [XS] X = XT 1 + S Ks Represní vstupní funkce promotoru • předpokládáme promotor P a represor X aktivovaný signálem S • dána aktivitou P regulovanou aktivním represorem X: f − P (X) = α− P (X) = βmax 1 + X Kd • závislost na vstupním signálu (konst. koncentrace XT ): f − P (S) = βmax 1 + XT Kd + Kd S Ks • závislost na vstupním signálu a koncentraci aktivního represoru: f − P (X, S) = βmax 1 + X+[XS] Kd + Kd S Ks Model regulace represoru Signal S X Y P gen y S:X r1 Y X r2 S XS r3 r4 reakční síť: r3, r4 : X + S ↔ [XS] r1 : → Y ; X− (Kd ) r2 : Y → model dynamiky: dY dt = f − (S, X) − γY Model regulace represoru Signal S X Y P gen y S:X r1 Y X r2 S XS r3 r4 reakční síť: r1 : → Y ; X− (Kd ) r2 : Y → zjednodušený model dynamiky (pro S konstantní): dY dt = f − (X) − γY Pozitivní regulace (aktivace) transkripce • transkripční faktor X aktivuje transkripci proteinu Y • řízeno zvýšením řídícího signálu S (X vázáno v X : S) gen Y X dY dt = −γY Pozitivní regulace (aktivace) transkripce • transkripční faktor X aktivuje transkripci proteinu Y • řízeno zvýšením řídícího signálu S (X vázáno v X : S) gen Y X*X signal Sx dY dt = −γY Pozitivní regulace (aktivace) transkripce • transkripční faktor X aktivuje transkripci proteinu Y • řízeno zvýšením řídícího signálu S (X vázáno v X : S) Y Y Y Y gen Y X*X signal Sx X* transkripce RNAP probiha dY dt = β − γY Aktivační vstupní funkce promotoru • předpokládáme promotor P a aktivátor X aktivovaný signálem S • aktivita promotoru P regulovaného aktivátorem X: α+ P (X) = βmax X Kd + X • dána aktivitou P regulovanou aktivním aktivátorem X: f + P (X∗ ) = α+ P (X∗ ) = βmax X∗ Kd + X∗ Aktivita promotoru v závislosti na regulujícím aktivátoru [X] βmax tc Kd α− P (X) Aktivace aktivační vstupní funkce • aktivátor X je účinný pouze ve vázané (aktivní) formě (X∗) • vazba se signálním komplexem S • aktivní forma represoru je zcela v komplexu [XS] S:...:S:X X Signal S Y P gen y Aktivace aktivační vstupní funkce • aktivátor X je účinný pouze ve vázané (aktivní) formě (X∗) • vazba se signálním komplexem S • aktivní forma represoru je zcela v komplexu [XS] f + (S) = βmax X∗ Kd + X∗ X∗ = [XS] = XT S Kx + S Aproximace pro konstantní S: f + P (X) ≈ βmax X Kd + X Příklad modelu regulace P2 Crp Fis P1 crp Crp:S Signal S Příklad modelu regulace P2 Crp Fis P1 crp Crp:S Signal S • reakční síť: S + Crp ↔ [Crp : S] → Crp; Fis, Crp Crp → Příklad modelu regulace P2 Crp Fis P1 crp Crp:S Signal S • reakční síť: S + Crp ↔ [Crp : S] → Crp; Fis, Crp Crp → • model: d[Crp] dt = βmax P1 f − (Fis) + βmax P2 f − (Fis)f + (Crp : S) − γ[Crp] Příklad modelu regulace P2 Crp Fis P1 crp Crp:S Signal S • reakční síť: → Crp; Fis, Crp Crp → • zjednodušený model (S konstantní): d[Crp] dt = βmax P1 f − (Fis) + βmax P2 f − (Fis)f + (Crp) − γ[Crp] Vícerozměrné vstupní funkce (AND,OR) Vícerozměrné vstupní funkce Obsah Aproximace vstupních funkcí • vzhledem k S-charakteru vstupní funkce lze uplatnit její aproximaci pomocí schodové funkce Aproximace vstupních funkcí • vzhledem k S-charakteru vstupní funkce lze uplatnit její aproximaci pomocí schodové funkce f + (X) ∼ βs+ (X, K) Aproximace vstupních funkcí • vzhledem k S-charakteru vstupní funkce lze uplatnit její aproximaci pomocí schodové funkce f + (X) ∼ βs+ (X, K) f − (X) ∼ βs− (X, K) Aproximace vstupních funkcí • vzhledem k S-charakteru vstupní funkce lze uplatnit její aproximaci pomocí schodové funkce f + (X) ∼ βs+ (X, K) f − (X) ∼ βs− (X, K) s+ (X, K) = 1, if X > K, 0, if X < K, Aproximace vstupních funkcí • vzhledem k S-charakteru vstupní funkce lze uplatnit její aproximaci pomocí schodové funkce f + (X) ∼ βs+ (X, K) f − (X) ∼ βs− (X, K) s+ (X, K) = 1, if X > K, 0, if X < K, s− (X, K) = 1 − s+ (X, K) Aproximace vstupních funkcí • vzhledem k S-charakteru vstupní funkce lze uplatnit její aproximaci pomocí schodové funkce f + (X) ∼ βs+ (X, K) f − (X) ∼ βs− (X, K) s+ (X, K) = 1, if X > K, 0, if X < K, s− (X, K) = 1 − s+ (X, K) • aproximace odpovídá zavedení tzv. “kinetické logiky” Diskretizace vstupní funkce (aktivátor) Step Function Diskrétní aproximace in silico modelů – přehled boolean networks kinetic logic continuous logic differential equations discrete continuous