Szimulációk egyszerősített fehérjemodellekkel Szilágyi András
Szimulációs módszerek alkalmazhatósági tartományai
Egyszerősített modellek Három típusát mutatjuk be: • „Játék” rácsmodellek • Realisztikusabb rácsmodellek • Diszkrét molekuladinamika
Kitérı: A fehérjék felgombolyodása DNSRNSpolipeptidláncmőködıképes fehérjemolekula transzkripció
transzláció
felgombolyodás
Anfinsen kísérletei (1961): Az aminosavsorrend tartalmazza a háromdimenziós szerkezet kialakulásához szükséges összes információt. A natív szerkezet a termodinamikailag legstabilisabb állapot. A fehérjék stabil konformációs állapotai: • A legombolyodott (denaturált) állapot. Jele: U (unfolded) vagy D (denatured) • A felgombolyodott (natív, általában biológiailag aktív) állapot. Jele: F (folded) vagy N (natív). • A köztes állapotok (köztitermékek, intermedierek). Jelük: I (intermediate)
Állandó nyomáson és hımérsékleten egy rendszer egyensúlyi állapotát a szabadentalpia (G) minimuma adja
Az entalpiacsökkentı folyamatok (pl. kölcsönhatások létrejötte) és az entrópianövelı folyamatok (pl. egy láncmolekula mozgásszabadságának megnövekedése) egyaránt csökkentik a szabadentalpiát. A hımérséklet határozza meg, hogy melyik típusú folyamat érvényesülhet jobban (a magasabb hımérséklet az entrópianövelı folyamatoknak kedvez).
A felgombolyodás termodinamikája
Alacsony hımérsékleten a natív állapot (sok kölcsönhatás, alacsony entrópia), magas hımérsékleten a denaturált állapot (kevés kölcsönhatás, magasabb entrópia) a stabilabb. Egy bizonyos hımérsékleten (Tm „olvadáspont”) kooperatív átmenet történik a két állapot között (jellegzetes csúcs a hıkapacitásgörbében). (A hımérsékletet eléggé lecsökkentve általában megint a denaturált állapot lesz a stabilabb hidegdenaturáció)
Rácsmodellek
Rácsmodell: "önelkerülõ" gyöngysor, rácspontokon
2D HP modell
3D HP modell Ken A. Dill (UCSF)
Perturbált heteropolimer modell, kétbetûs kóddal Eugene Shakhnovich (Harvard)
A rácsmodellek típusai HP modell (Hidrofób-Poláros) [Ken A. Dill és mások] • Kétféle aminosav: Hidrofób és Poláros (kétbetős kód) • 2D és 3D modellek. Eredményeik hasonlóak; a 2D modellek könnyebben számíthatóak. Hidrofób maggal rendelkezõ modell 2D-ben 16, 3D-ben 64 aminosavból készíthetõ. • Egyparaméteres: ε kontaktenergia (a HH kontaktus energiája) • A kontaktusok energiája:
(EHH, EHP, EPP) = (−ε, 0, 0) • Perturbált heteropolimer modellek [Eugene I. Shakhnovich és mások] • Erıs vonzás a monomerek között, gyenge variációval • Többnyire 27-mer, alapállapota 3x3x3-as kockán (103 346 lehetséges állapot) • Két változat: • Független kontaktenergiák: E(i,j) •
= Bij Kétbetős kód: A,B. Ált. EAB > EAA, EBB. (Ált. nem mag-felszín szeparációhoz, hanem jobb-bal elkülönüléshez vezet.)
A rácsmodellek elınyei és hátrányai Elınyök • Szekvenciatér és konformációtér teljességében felderíthetı • Egzaktul kiszámítható mennyiségek! (Jobb az egyszerő, de egzaktul számolható modell, mint a bonyolult, de csak közelítıleg számítható.) • Nem csak fehérjékre vonatkoztatható • Használható analitikus elméletek tesztelésére (mint spinüveg-modell, átlagtér-elméletek) • Használható konformációkeresı algoritmusok tesztelésére Hátrányok • Az atomi felbontás hiánya • A szerkezeti és energetikai részletek nem pontosak • A lánc sokszor irreálisan rövid
A rácsmodellek elemzése • Enumeráció: a modell összes lehetséges állapotának elıállítása, megszámlálása és jellemzése. Kisebb modelleknél járható. • Monte Carlo mintavétel, ill. szimuláció: Monte Carlo módszerrel (véletlenszerő generálás) véletlen állapotok elıállítása, ebbıl statisztikai számítás. Nagyobb modelleknél. (Hibalehetıségek vannak, pl. rossz becslésekbıl.) • A natív állapot megkeresése: Kisebb modelleknél az összes lehetséges állapot végigpróbálásával, nagyobbaknál becslés speciális algoritmusokkal. A HP modellek alapállapotának megkeresése ún. NP-nehéz feladat
Egyszerő rácsmodellek: 1. Négytagú játékmodell
A: energiaszintek, B: állapotok, C: h kontaktushoz tartozó állapotok száma, g(h) • Legegyszerőbb 2D HP rácsmodell, szekvencia: HPPH • h db kontaktusnál az energia: −hε
A: Natív populáció részaránya, B: ∆G felgombolyodási szabadentalpia. Tm olvadáspontnál a natívhányad 50%, a szabadentalpia nulla. A szabadentalpiafüggvény lineáris (valódi fehérjéknél görbült)
A: energia, B: hõkapacitás. A hõkapacitásgörbe jellegzetes átmeneti csúcsot mutat
Már ez a legegyszerőbb modell is mutatja a fehérjék egyes termodinamikai tulajdonságait
2. Hattagú játékmodell
Hattagú 2D HP modell, szekvencia: HPHPPH. Három energiaszint: három állapot
Populációeloszlások a hõmérséklet függvényében: A populáció a tisztán natív állapotból indulva a köztes állapoton át folyamatosan áttolódik a denaturált állapotra
Felgombolyodási szabadentalpia a hõmérséklet függvényében: Az intermedier állapot miatt már nem lineáris, hanem görbül, mint a valódi fehérjéknél. Ha epszilont hõmérsékletfüggõvé tesszük (a hõmérséklet csökkenésével csökken, akárcsak a hidrofób kölcsönhatás), a szaggatott vonallal jelzett görbét kapjuk. Ez mutatja a valódi fehérjékre jellemzõ, parabolaszerû görbét, a hidegdenaturációval.
Egy egyszerû hattagú modell már a fehérjék szabadentalpia-görbéjét is jól modellezi
3. Húsztagú modellek A szekvenciától függõen sokféle viselkedés: A: natív részarány, B: hõkapacitás, epszilon/kT függvényében, C: a modellek natív állapotai (legmélyebb energiaszint)
ε a hidrofób kölcsönhatás erısségét jellemzi, ε/kT a denaturálószer-koncentrációnak feleltethetı meg (i): két csúcsú széles átmenet: kétdoménes fehérje (ii), (iii): egyetlen éles csúcs
(ii) populációeloszlása a hımérséklet függvényében: T=0,3 és T=0,4-nél az eloszlásnak minimuma van: kooperatív, kétállapotú rendszer
Egy húsztagú modell már kétállapotú kooperativitást képes mutatni, akárcsak a valódi fehérjék A modell natív állapota egy β-lemezt és két hélixet tartalmaz, poláros felszínnel
Rácsmodell felgombolyodásának szimulációja Monte Carlo módszerrel (ld. elızı elıadás) lehetséges, ehhez definiálni kell egy mozgáskészletet (move set). Példa:
Egy 13 tagú HP modell felgombolyodásának szimulációja
Két útvonal, egy gyors és egy lassú. A lassú útvonalon kinetikus csapda (lokális energiaminimum). Ez az egyszerő modell is komplex felgombolyodási kinetikát mutat (a lizozimhoz hasonló).
Egy 13 tagú HP modell felgombolyodásának szimulációja (folytatás) A populációk idıfüggése: (A számok az adott állapotban meglévı HH kontaktusok számát jelentik. Szaggatott vonal: a hidrofób kontaktusok száma a natív szerkezethez viszonyítva)
Gyors hidrofób kollapszus, majd 5 nagyságrenddel hosszabb idıskálán lassú átrendezıdés a natív állapotba. Jól mutatja a valódi fehérjéknél megfigyelt sajátosságokat.
Szabadenergia-felületek rácsmodellekbıl
Rácsmodellek Monte Carlo szimulációiból számítható a szabadenergia:
F = E − TS ahol • E az energia (a kontaktusok számából) • S az entrópia (az adott energiájú állapotok számából számítható). A szabadenergiát megfelelõ mennyiségek mint koordináták függvényében ábrázolva szemléletes felületek adódnak. Ilyen mennyiségek pl.: • C: az összes kontaktusok száma • Q0: a natív kontaktusok száma (olyan kontaktus, ami a natív szerkezetben is megvan)
Egy 27 elemő modell szabadenergia-felülete: A szimuláció menete ezen a felületen szemléltethetõ: • A szimuláció random szerkezetbõl indul (kevés kontaktus): hozzáférhetõ állapotok száma kb. 1016 • Gyorsan összeesik egy kompakt szerkezetté, mely a lehetséges kontaktusok 60%-át tartalmazza, de a 60%-nak csak 25%-a natív kontaktus. A hozzáférhetõ állapotok száma kb. 1010. • Itt egy széles minimum van a szabadenergiában • Ezután lassú folyamat következik: a lánc keresi a natív állapot felé vezetõ átmeneti állapotot • Átmeneti állapotok száma 103 • Az átmeneti állapotból gyorsan betalál a natívba (1 állapot) • Zöld, piros: lehetséges, különbözõ útvonalak. Sárga: átlag.
125 elemő modell (5x5x5-ös kocka a natív állapot) szabadenergia-felülete Két mennyiség mint koordináta: Qc: a belsõ magban lévõ natív kontaktusok (34 darab), Qs: olyan kontaktusok, amelyek a legfontosabb intermedierben nincsenek benne (15 darab) A felgombolyodás menete: • Gyors összeesés, majd lassú keresés • A molekulák 15%-ában gyorsan kialakul a natív mag, majd közvetlenül a natív állapotba jut (sárga útvonal) • A molekulák 40%-a különféle nemnatív kontaktusokat létesít, így egy hosszú élettartamú intermedier alakul ki (piros útvonal), melyben két szubdomén van, egyikben a natív 5x5x5-ös kocka felsõ három síkja megvan, a másik a többi részbõl áll. Mindkét szubdomén eléggé felgombolyodott, de a köztük lévõ kapcsolatok nem jöttek létre. • A maradék 45%-ban szintén rosszul kezd felgombolyodni, de gyorsabban áthidalja a dolgot (ezt nem jelöltük)
A rácsmodellek következtetései • Bizonyos szekvenciájú rácsmodellek stabil, kompakt állapotokba esnek össze. Ezek számos tekintetben igen hasonlóak a fehérjékhez: • Kompakt, unikális natív szerkezet • Másodlagos és harmadlagos szerkezet • Apoláros mag, poláros felszín • Éles, szigmoidális, kooperatív átmenet • Kétállapotú kooperativitás • Kompakt, komplex denaturált állapotok • Többfázisú felgombolyodási kinetika • Szabadenergiafelületek, kinetikus komplexitás jól modellezhetõ
Realisztikusabb rácsmodellek • Elég sőrő 3D rácson, ha nem kötjük ki, hogy a kötések rácséleken feküdjenek, valóságos fehérjeszerkezet is reprezentálható • Minden mozgás a rácsra korlátozódik, ezért Monte Carlo módszerrel hatékonyan szimulálható • Tudás alapú (empirikus) potenciálfüggvényekkel
CABS modell (ZhangKolinski-Skolnick 2003) Atomok egy 0,87 angström cellamérető rácson Csak Cα, Cβ és az oldallánc tömegközéppontja
Monte Carlo mozgáskészlet
Empirikus potenciálfüggvény
A paramétereket az ismert szerkezetek adatbázisából vezették le (ezért tudásalapú)
A potenciálfüggvény 19 szabad paramétert (súlyt) tartalmaz. Ezeket optimalizálták egy decoy (csalétek) halmaz segítségével oly módon, hogy az energia és az RMSD közötti korreláció maximális legyen. „Tölcséresített potenciálfüggvény”: a szerkezetet beirányítja a globális minimumba. CASP2004-en a CABS modellen alapuló szerkezetpredikciós módszer kiválóan szerepelt. 100 aminosav alatti fehérjéket általában a helyes szerkezetbe gombolyít fel (NMR szerkezetnek megfelelı pontossággal).
Diszkrét molekuladinamika (DMD) • A hagyományos molekuladinamikában a Newton-féle mozgásegyenleteket integráljuk. Az idıbeli lépésköznek elegendıen kicsinek kell lennie, a molekula leggyorsabb mozgásaihoz igazodva. • DMD: a dinamikát ütközési események egymásutánjaként fogjuk fel. Két ütközés között az atomok egyenes vonalú, egyenletes mozgást végeznek, ezt nem szükséges követni, a számítás rögtön a következı ütközésre léptetheti a rendszert. • Ezáltal több nagyságrenddel gyorsabban végezhetı a szimuláció! • Ennek ára: a potenciálfüggvényt diszkretizálni kell.
A potenciálfüggvény diszkretizálása A párpotenciált lépcsızetes függvénnyel közelítjük
DMD szimuláció 1
Ballisztikus szakasz
2
Ütközés
v′j
v′i vi
vj
r = rj − ri v = vj − vi
|r + vτ| = σ
b = rv r = |r| v = |v|
|r| = σ
A következı ütközés ideje az alábbi egyenlet két megoldása közül a kisebb pozitív:
−b ± b 2 − v 2 (r 2 − σ 2 ) τ= v2
Az ütközések feldolgozása Általánosított ütközések: potenciálfallal
i
1
mi∆vi = −mj∆vj
lendület
2
Ek′ + ∆U = Ek
energia
j kilépés
s =1
belépés
s =1
∆vi = mjΦr ∆vj = −miΦr
potenciálgödör
megoldható m=
mi + m j 2
visszapattanás
s = −1
visszapattanás
s = −1
Φ=
U ∆U
r
2 4 σ m∆U −b + s b 2 − mi m j
2σ 2 m
Eseményfa Az ütközési eseményeket bináris keresıfában tároljuk root
t = 23.7 A=4 B=3 : t = 19.2 A = 9 B =1 :
t = 25.1 A=6 B =8 :
t = 24.2 A=3 B =8 :
t = 23.9 A=9 B=2 :
t = 24.7 A = 4 B = 10 :
t = 32.1 A =1 B = 7 :
t = 31.6 A=2 B=4 :
t = 43.8 A=4 B=5 :
t = 47.8 A=6 B=9 :
Egyszerősített fehérjemodell N
Négygolyós modell
C
Négyféle egyesített atom
Cα R = Cβ
R O O
H
H Cα
C N
C
N
Cα
C
H H
O
R alanin: R = Cβ
Peptidgerinc
kovalens kötés effektív kötés Cβi
A kötésszög-kényszereket 1-3 effektív kötésekkel valósítjuk meg
Ni
Φi
Ci+1
Ni+1
Cαi
Ψi
Ci
Cα, i+1
Kovalens kötés potenciálja:
U rid
r
Cβ, végtelen mély derékszögő gödör ±2% fluktuációt engedünk meg az ideális kötéshossz körül
i+1
Hidrogénkötés
Cα,
Cαi
i-1
Segédkölcsönhatásokat vezetünk be az orientációfüggés modellezésére
U
Ni
4.0Å – 4.2Å
ε HB / 2 ε HB / 2 rmin
r2
0,
rmax > r > r2
Cj
Nj+1
Cαj
rmax
r
r1
ε HB / 2 r2 > r > r1 V (r ) = ε HB r1 > r > rmin +∞
otherwise
Hımérsékletszabályozás Andersen-termosztát:
Véletlenszerő ütközések történnek a hıfürdı „szellemrészecskéivel” Két ilyen ütközés között eltelı idı valószínőségeloszlása:
P(t ) = ν e −ν t
ν
az ütközésgyakoriság, ezt úgy választjuk meg, hogy az események 1%-a legyen szellemütközés
A kinetikus energia fluktuációi jó egyezést mutatnak az elméleti értékekkel
Hexadeka-alanin
A hidrogénkötések számának fluktuációja 330 Kelvinen, 1 mikroszekundum idıtartam szimulációja során
Másodlagosszerkezet-tartalom a hımérséklet függvényében
A DMD perspektívái • Mostanában újra felfedezett módszer • Korlátja: bonyolultabb potenciálfüggvények diszkretizálása nehezen megoldható • Jól alkalmazható: peptidek és kisebb fehérjék konformációs terének feltérképezésére, asszociáció, aggregáció modellezésére
