Monte karlo metode I primene u bioinformatici master rad

5. MONTE KARLO METOD RAZMENE KOPIJA - PRIMENA NA PROBLEM UVIJANJA PROTEINA

5.1 Uvijanje proteina

Protein je najvažniji molekul u klasi bioloških makromolekula. Oni učestvuju u svim ćelijskim i među-ćelijskim procesima i obavljaju širok spektar funkcija unutar živih organizama. Proteini su u suštini polimeri amino-kiselina koje su poređane u linearne lance i spojene međusobno peptidnim vezama. Struktura proteina je određena redosledom amino-kiselina u polipeptidnom lancu i od nje direktno zavisi funkcija proteina. Sekvenca amino-kiselina u proteinu definisana je u genima i sadržana u genetskom kodu. Genetski kod određuju dvadeset "osnovnih" amino-kiselina.

Zbog specifičnog vezivanja lanaca amino-kiselina, proteini imaju četiri strukturna nivoa koji određuju njihov izgled u prostoru (konformaciju). U zavisnosti od količine unakrsnih povezivanja struktura proteina može biti: primarna, sekundarna, tercijarna i kvaternarna [27].

• 
  Primarna struktura
  

Ovo je najjednostavnija struktura sa minimalnim brojem unakrsnih povezivanja. Primarna struktura predstavlja redosled amino-kiselina u polipeptidnom lancu i određena je genom koji odgovara proteinu. Specifičan niz nukleotida u DNK je prepisan na iRNK, koji „čitaju“ ribozomi u procesu translacije - biosinteze proteina. Niz amino-kiselina je jedinstven za taj protein i definiše njegovu strukturu i funkciju.

• 
  Sekundarna struktura
  

Sekundarne strukture predstavljaju izgled polipeptidnih lanaca u prostoru. Sastoje se od dva ili više polipeptidnih lanaca. Obrazuje se uglavnom na bazi stvaranja vodoničnih veza između atoma koji se nalaze u sastavu peptidne veze i relativno blizu u polipeptidnom lancu. Najčešće sekundarne strukture su:

• 
  -heliks: tip sekundarne strukture u kojem se deo polipeptidnog lanca uvija u spiralu, najčešće u desnu zavojnicu.
  
• 
  -ravan: tip sekundarne strukture u kojem se dva polipeptidna lanca ili delovi istog polipeptidnog lanca svrstavaju paralelno jedan prema drugom, lanci mogu biti paralelni ili antiparalelni.
  

• 
  Tercijarna struktura
  

Polipeptidni lanac, koji u sebi već sadrži delove sa sekundarnom strukturom, sposoban je da se u celini izuvija u prostoru i zauzme položaj koji odgovara najstabilnijoj konformaciji, pri datim uslovima, koji se naziva tercijarna struktura. Ova struktura se formira kada se sekundarne strukture uvežu i formiraju trodimenzionalnu kompleksnu strukturu prvenstveno kroz hidrofobne interakcije, ali i vodonične veze, jonske interakcije i disulfidne veze. Tercijarna struktura zaokružuje sve nekovalentne interakcije koje ne razmatra sekundarna struktura, definiše sveukupno savijanje u proteinu i utiče na funkciju koju će protein obavljati.

• 
  Kvaternarna struktura
  

Kvaternarna struktura je nivo organizacije proteina koji je sačinjen od više polipeptida, tako dobijene strukture se nazivaju kompleksi ili agregati. Svaki polipeptid (sa uređenom tercijarnom strukturom) predstavlja jednu podjedinicu. Kvaternarna struktura je prostorni raspored podjedinica u složenoj celini.

Proces kojim se formiraju više strukture se zove uvijanje proteina i predstavlja jedan od najznačajnijih problema savremene biologije. Korektna trodimenziona struktura je esencijalna za funkciju, mada neki delovi funkcionalnih proteina mogu da ostanu neuvijeni. Prirodno stanje (eng. native state) proteina predstavlja oblik proteina koji poseduje najmanje slobodne energije ⁶ [36]. Posledica neuspešnog uvijanja u prirodnu strukturu su neaktivni proteini koji su obično toksični. Iako svaki jedinstveni polipeptid može imati više od jedne stabilne uvijene konformacije, svaka konformacija ima svoju biološku aktivnost i samo jedna konformacija se smatra aktivnom. Funkcija je povezana sa prirodnim stanjem proteina i kada se ova struktura prekine protein nije u stanju da izvršava svoju specifičnu funkciju. Međutim, postoje proteini koji nisu uvijeni u prirodno stanje (tzv. neuređeni proteini) koji imaju značajnu ulogu u ćelijskim procesima, skladištenju i obradi informacija, ali dovode i do nekih ozbiljnih kardiovaskularnih oboljenja, karcinoma i neurodegenerativnih bolesti.

Mehanizam uvijanja proteina nije u potpunosti razjašnjen. Uspešan metod za predviđanje strukture proteina bi imao dalekosežnu primenu i u drugim naučnim oblastima uključujući i genetiku i medicinu. Trenutne metode koje se koriste u labaratorijama su uglavnom skupe i zahtevaju dosta vremena kako bi se dobili rezultati. U današnjoj eri savremenih tehnologija nije prihvatljivo oslanjati se samo na tehnike predviđanja strukture, kao što su rendgenska kristalografija i nuklearna magnetna rezonanca, već se javlja potreba za otkrivanjem novog efektivnog i efikasnog algoritma za predviđanje strukture proteina. Međutim, čak i za pojednostavljen model proteina, koji će biti objašnjen u nastavku, ovaj problem se pokazao kao NP kompletan problem za čije rešavanje ne postoji algoritam u polinomijalnom vremenu [45]. U ovom radu će biti objašnjena Monte Karlo metoda razmene kopija, koja se veoma uspešno primenjuje na složene proteinske modele i druge probleme optimizacije, pomoću koje ćemo predvideti strukturu proteina pronalaženjem one sa najmanjom energijom koja predstavlja prirodno stanje.

5.2 Monte Karlo metoda razmene kopija

5.2.1. Uvod

Monte Karlo razmena kopija, skraćeno MKRK (eng. replica exchange Monte Carlo), poznata još i kao metoda paralelnog kaljenja (eng. Parallel tempering), predstavlja metodu simulacije čiji je glavni cilj poboljšanje dinamičkih osobina Monte Karlo metode uzorkovanja u fizičkim sistemima, a koja se zasniva na posmatranju M odvojenih sistema (kopija) na različitim temperaturama. Ova Monte Karlo simulacija, koja koristi Metropolis-Hestingov algoritam za ažuriranje stohastičkih procesa, ocenjuje energiju sistema i prihvata ili odbacuje posmatranu konformaciju u zavisnosti od temperature sistema T. Kod sistema u kojima su uzorci u velikoj korelaciji posmatrane konformacije se uglavnom odbacuju, pa se za njih kaže da se kritično usporavaju [28].

Skup konformacija za datu sekvencu predstavlja skup svih mogućih strukturnih stanja bez raskidanja kovalentnih veza, pa je moguće odrediti finalnu konformaciju korišćenjem grube sile, tj. računanjem svake konformacije i pronalaženjem one sa najmanjom energijom. Kako broj konformacija eksponencijalno raste sa brojem amino-kiselina u sekvenci, ovu metodu je moguće koristiti samo na kraćim sekvencama. Stoga je potrebno uprostiti prostor pretraživanja na uštrb tačnosti, kako bi se uštedelo na vremenu potrebnom za izračunavanje. Pažljivim odabirom temperatura i broja kopija sistema mogu se poboljšati svojstva Monte Karlo simulacije.

Razumevanje mehanizma uvijanja proteina predstavlja jedna od najizazovnijih problema molekularne biologije. U ovom odeljku će biti predstavljena metoda Monte Karlo razmene kopija koja za generisanje konformacije sistema koristi metodu molekulske dinamike (skraćeno MD), kao i njenu primenu na problem uvijanja penta-peptida⁷ Met-enkefalina u gasnom stanju.

5.2.2 Molekulska dinamika

Molekulska dinamika (MD) je oblik računarske simulacije, gde atomi i molekuli mogu da interaguju u određenom vremenskom intervalu, pokoravajući se poznatim zakonima fizike. Ovaj metod su prvi predstavili Alder i Vajnvrajt (eng. Alder and Wainwright) [29], nakon čega on postaje jedan od korišćenijih alata za istraživanje kompleksnih fizičkih sistema. Molekulska dinamika predstavlja determinističku proceduru za integrisanje jednačina kretanja (Hamiltonove jednačine) na osnovu klasičnih principa mehanike. Prvi korak ove simulacije je podešavanje kvantitativnog sistema (modela) pod datim uslovima (npr. tačan broj čestica ili ukupna energija sistema), nakon čega se na osnovu Njutnovih zakona kretanja generiše konfiguracija sistema, kao funkcija vremena. Podaci dobijeni ovom simulacijom predstavljaju snimke pozicija i brzine čestica na osnovu kojih se može odrediti "tipična karakteristika" sistema.

Glavni zadatak molekulske dinamike je integracija jednačina kretanja u datom vremenskom periodu, kao i proučavanje fizičkih i hemijskih osobina sistema u tom vremenskom periodu (npr. uticaj vode pri uvijanju proteina). Neka je -dimenzioni vektor pozicije čestice u vremenu , = () vektor mase, a brzina čestice, tada možemo definisati impuls i kinetičku energiju sistema kao:




Neka je potencijalna energija sistema, tada je ukupna energija čestica sistema:



a Hamiltonove jednačine su sledećeg oblika:





Kako kapacitet kompjutera omogućava samo diskretne operacije, što znači da nije moguće neprekidno izračunavanje jednačina kretanja čestica, u praksi se koristi aproksimacija Hamiltonovih jednačina u vidu Tejlorovog razvoja:





Jedan od najjednostavnijih i najčešće korišćenih algoritama za integraciju jednačina kretanja je Verletov algoritam [30], koji se zasniva na zapažanju da je:

Za izabrani vremenski korak , pozicija i impuls čestice imaju sledeće vrednosti:

Još jedna MD metoda koja je često u upotrebi je metoda preskakanja (eng. leap frog) [31]. Glavna karakteristika ove metode je to da se impuls izračunava na polovini vremenskog intervala, pa su jednačine sledećeg oblika:

Molekulski sistemi se generalno sastoje od velikog broja čestica, stoga je nemoguće analitički pronaći osobine tako složenih sistema, ali molekulska dinamika prevazilazi ovaj problem koristeći numeričke metode. Dizajn MD simulacija određen je mogućim kapacitetom kompjutera. Veličine simulacije (N = broj čestica), vremenski korak i ukupno vreme trajanja simulacije treba odabrati tako da se proračun završi u nekom razumnom vremenskom periodu. Ipak, simulacija treba da bude dovoljno duga da bi verno predstavila prirodni proces koji proučavamo. Većina naučnih radova na temu dinamike proteina i DNK koriste simulaciju za procese koji u prirodi traju od nekoliko nanosekundi do mikrosekundi, a da bi se ove simulacije posmatranog procesa izvršile potreban je vremenski period od nekoliko dana do nekoliko godina. Glavna prednost MD simulacije u fizičkim sistemima je njena zasnovanost na osnovnim fizičkim principima (npr. Njutnovim zakonima), ali jedan od glavnih nedostataka je to što ona zahteva veoma mali vremenski korak koji se određuje na osnovu dužine trajanja najbržeg pokreta tokom posmatranog procesa. Prema tome vremenski korak mora biti adekvatno izabran kako bi se posmatrani proces pravilno simulirao, jer previše veliki vremenski korak dovodi do veće greške koja se javlja pri izračunavanju jednačina kretanja, a suviše mali vremenski korak utiče na efikasnost simulacije jer zahteva veći broj izračunavanja. Tipičan vremenski korak u klasičnoj MD je reda veličine femtosekunde (fs), pa bi na primer, MD simulacija uvijanja proteina koja u prirodi traje oko 1 milisekunde, trajala oko 10⁶ dana.

U standardnim Monte Karlo simulacijama Metropolis-tipa predložena raspodela se ne može jednostavno uklopiti u lokalnu dinamiku ciljane raspodele. Na primer, ukoliko se posmatrani sistem sastoji od gusto zbijenih čestica, tada će se izračunati pomeraj čestice uglavnom odbacivati jer će nova pozicija biti delimično zauzeta od strane drugih čestica. Kako bi se prevazišao ovaj problem Čarls Gajer (eng. Charls Geyer) [33] je predstavio metod MKRK koja kombinuje osnovnu ideju MD sa Metropolis pravilima prihvatanja radi određivanja uzoraka željene raspodele.

5.2.3 Softverski paketi za simulaciju molekulske dinamike

Molekulska dinamika nam omogućava izučavanje dinamike velikih makromolekula, uključujući i biološke sisteme kao što su protein, DNK, RNK, membrane i dr. Danas je u širokoj upotrebi u farmaceutskoj industriji pri izradi lekova radi testiranja osobina molekula bez potrebe njihove sinteze koja je jako skupa. Dinamički događaji mogu imati ključnu ulogu u kontrolnim procesima koji utiču na funkciju biomolekula. Postoji nekoliko softverskih paketa koji se koriste za simulaciju molekulske dinamike. Svaki od njih ima različite karakteristike i pravila upotrebe, a ovde ćemo predstaviti tri najpopularnija paketa: AMBER, CHARm i Gromacs [39].

AMBER (eng. Assisted Model Building and Energy Refinement) predstavlja uopšten naziv za skup programa koji korisnicima omogućava da sprovedu i analiziraju MD simulaciju, posebno za proteine, amino kiseline i karbohidrate. Ovi programi se mogu podeliti u tri grupe: programi za pripremu, programi za simulaciju i programi za analizu. Osnovni programi za pripremu su Antechamber i LEaP. Antechamber automatizuje proces razvoja deskriptora polja sile⁸ za većinu organskih molekula, koji počinje sa određenom strukturom (obično u PDB⁹ formatu) i generiše fajl koji se kasnije može koristiti u LEaP radi modeliranja molekula. Deskriptor polja sile je dizajniran tako da bude kompatibilan sa standardnim AMBER poljima sile za protein i amino kiseline. LEaP je program koji obezbeđuje osnovnu izgradnju modela i kreiranje AMBER koordinata i parametara ulaznog fajla. On u sebi sadrži editor koji omogućava izgradnju ostataka i manipulaciju molekula. Glavni program za simulaciju molekulske dinamike predstavlja Sander koji se takođe koristi i kod metode razmene kopija, termodinamičke integracije i dr. Ptraj je program koji pripada grupi programa za analizu i koristi se za analiziranje MD trajektorija, hidrogenih veza i sl. U Tabeli 1 su predstavljene prednosti i mane AMBER programa za simulaciju [39].

Prednosti	Nedostaci
Obezbeđuje podršku za simuliranje karbohidrata, proteina i manjih organskih molekula.	Nije moguće simulirati samo jedan deo sistema, npr. samo aktivna strana enzima.
Računa slobodnu energiju koristeći termodinamičku integraciju ili kišobran tehniku uzorkovanja (eng. umbrella sampling)	Programske komponente nemaju korisnički interfejs.
Radi ubrzanja konvergencije moguće je koristiti metodu zamene kopija.	Nedostatak Monte Karlo uzorkovanja, dinamike torzionih uglova, izračunavanje “dualne topologije” slobodne energije.
Omogućava analizu trajektorija i fleksibilna ograničenja koja mogu biti zasnovana na podacima spektroskopije nuklearne magnetne resonance.	Kod je pisan od strane različitih autora tokom razvoja pa sadrži delove koji su teški za razumevanje i modifikaciju.
Ima veliku i aktivnu zajednicu korisnika, kao i tutorijale i uputstva za nove korisnike. Kod je prenosiv i moguće su dopune i modifikacije.	Korisnici moraju sami da kompaliraju program.

Tabela 1: Prednosti i nedostaci AMBER programa.
CHARMM je ime grupe široko korišćenih polja sila za molekulsku dinamiku, a isto tako i ime softverskog paketa za molekulsko dinamičke simulacije i analizu. CHARMM istraživački projekat uvrstava mrežu programera širom sveta koji rade na razvoju i održavanju programa. Licence za ovaj softver su dostupne besplatno individuama i grupama koji se bave akademskim istraživanjem. Accelrys distribuira komercijalnu CHARMM verziju, koja se zove CHARMm. CHARMM program omogućava izvođenje i analizu širokog kruga molekularnih simulacija. Najosnovnije vrste simulacija su minimizacija date strukture i računanje trajektorije molekulske dinamike. CHARMM je jedan od najstarijih programa molekulske dinamike i kao što je slučaj sa AMBER programa i ovaj program je pisan od strane različitih pojedinaca i grupa pa je kod težak za razumevanje i modifikaciju [40].

GROMACS (skraćeno od   Groningen mašina za hemijske simulacije (eng. GROningen MAchine for Chemical Simulations)) je molekulsko dinamički simulacioni paket originalno razvijen na Groningenskom univerzitetu. On se u današnje vreme održava i proširuje i na drugim mestima, uključujući Upsalski univerzitet, Stokholmski univerzitet i Maks Plank institut za istraživanje polimera. GROMACS je softver otvorenog koda pod GNU generalnom javnom licencom. GROMACS projekat je originalno započet da bi se konstruisao namenski paralelni računarski sistem za molekularne simulacije, koji je baziran na prsten strukturi. Izvorni kod specifičan za molekulsku dinamiku je prerađen u C programskom jeziku iz Fortran77-baziranog programa GROMOS [41].

Program je napisan za Unix operativne sisteme, ali se on može koristiti na Microsoft Windows mašinama koristeći Cigvin (engl. Cygwin) Unix sloj. GROMACS sadrži skript za konvertovanje molekulskih koordinata iz PDB fajla u formate koje program interno koristi. Nakon što je konfiguracioni fajl za simulaciju nekoliko molekula (po mogućnosti uključujući rastvarač) kreiran, izvršavanje simulacija proizvodi fajl sa trajektorijama koji opisuje kretanje atoma u toku vremena. Taj trajektorijski fajl se može analizirati ili prikazati brojnim alatima.

Mnogi specifični elementi su dodati u toku tranzicije GROMOS-a u GROMACS. Najznačajniji među njima su:

• 
  Generička reprezentacija svih mogućih tipova periodičnih kutija
  
• 
  Optimizovano rukovanje listom suseda putem smeštanja translacionih vektora ka najbližim susedima u periodičnom sistemu
  
• 
  Specijalizovane rutine za računanje inverznog kvadratnog korena
  
• 
  Korišćenje kubne splajn interpolacije iz tabeliranih vrednosti za evoluiranje sile/energije
  
• 
  Brza na-rešetki-zasnovana pretraga suseda
  

Visoko optimizovani kod čini GROMACS jednim of najbržih programa za molekulske simulacije. Dodatno, podrška za različita polja sila daje GROMACS-u veliku fleksibilnost.

5.2.4 Algoritam Monte Karlo Razmena Kopija

Posmatrajmo sistem od N atoma sa zadatom masom i pridruženim koordinatnim vektorom i impulsom p = . Ukupna energija sistema tj. Hamiltonijan je tada jednak:

gde je kinetička energija:

U kanonskom ansamblu¹⁰ na temperaturi T, svako stanje x = (q, p) sa Hamiltonijanom H(q, p) meri se Bolcmanovim faktorom:

gde je a Bolcmanova konstanta. Prosečna kinetička energija na temperaturi T je tada data kao:

Algoritam Monte Karlo metoda razmena kopija pretražuje opšti ansambl¹¹ od M neinteragujućih¹² kopija ili replika koje predstavljaju potencijalno rešenje problema. Svakoj od ovih kopija se pridružuje jedinstvena temperaturna vrednost Na taj način imamo 1-1 korespodenciju između kopija i temperature. Oznaka za kopije je permutacija oznake za temperaturu i obrnuto:

gde je funkcija permutacije od m, a njen inverz.

Neka X = predstavlja stanje sistema, gde je stanje i-te kopije na temperaturi Pošto su replike neinteragujuće, težinski faktor za stanje X u ovom sistemu predstavlja proizvod Bolcmanovih faktora za svaku repliku (ili na svakoj temperaturi):

Pretpostavimo da razmenjujemo par kopija i i j na temperaturama i , tim redom:

U ovom zapisu je korišćena još jedna funkcija permutacije:

pa se proces razmeme kopija može zapisati kao:

gde se za i mogu uzeti sledeće vrednosti, za koje se smatra da predstavljaju najjednostavniji i najprirodniji izbor:

Uočavamo da je ovaj proces ekvivalentan zameni para temperatura i za odgovarajuće replike i i j, kao što sledi:

Da bi proces zamene kopija konvergirao uravnoteženoj raspodeli dovoljno je nametnuti uslov detaljne uravnoteženosti verovatnoće prelaza

Predložena razmena kopija se prihvata ili odbacuje u zavisnosti od Metropolisovog kriterijuma:

gde je:

Simulacija MKRK se tada realizuje naizmeničnim izvođenjem sledeća dva koraka:

1. 
   Svaka replika u kanonskom ansamblu na fiksnoj temperaturi se simulira simultano i nezavisno za određeni MD korak.
   
2. 
   Par replika na susednim temperaturama, recimo i se zamenjuje sa verovatnoćom w() u jednačini Metropolisovog kriterijuma, prema tome algoritam MKRK je oblika:
   

U ovom pristupu koristimo MD algoritam u prvom koraku, dok se u drugom koraku zamenjuju samo replike koje odgovaraju susednim temperaturama, zato što se opseg prihvatljivosti zamene (eng. acceptance ratio) umanjuje eksponencijalno sa razlikom dveju β.

Glavna prednost MKRK u odnosu na druge metode jeste to da je težinski faktor apriori poznat, dok u drugim algoritmima određivanje težinskih faktora može biti prilično monotono i dugotrajno. Ipak, za optimalne performanse MKRK potrebno je izabrati prikladnu temperaturnu raspodelu.