DINAMICA MOLECULARA

DINAMICA MOLECULARA

In sensul cel mai larg dinamica moleculara se ocupa cu miscarea moleculara inerenta tuturor proceselor chimice; vibratia atomilor in molecula, reactiile chimice, interactia substrat-receptor, si multe alte procese complexe sunt asociate cu diferitele forme de miscare intra- si intermoleculare. 222c27c

In timp ce mecanismul unui proces chimic este descris de cinetica chimica, energetica acestor procese reprezinta obiectul termodinamicii. Cu alte cuvinte cinetica chimica este aceea care prezinta succesiunea si viteza evenimentelor care apar la transformarea sau trecerea moleculelor dintr-o stare in alta, iar termodinamica dicteaza relatiile energetice intre diferitele stari moleculare posibile.

III.1 Bazele dinamicii moleculare

Un sistem molecular, aflat la o anumita temperatura, este caracterizat cu acuratete de miscarea sa. Dinamica moleculara simuleaza aceasta miscare in spatiul fazelor (spatiul definit prin pozitia si vitezele atomilor), prin integrarea numerica a ecuatiilor de miscare date de Newton. Mecanica moleculara modifica gradele de libertate intramoleculare intr-o maniera asemanatoare cu cea folosita in minimizarea energiei: etapele individuale in minimizarea energiei constau in cautarea directiei de coborare catre un minim; etapele individuale in simularile dinamicii moleculare constau in calcularea noilor pozitii r_i+1 si viteze v_i+1 ale atomilor, pe baza celor curente (r_i, v_i) si a modificarilor in timp.

Viteza si directia de deplasare este determinata de fortele pe care atomii sistemului le exercita unul asupra celuilalt, forte descrise de ecuatiile lui Newton.

In practica, atomilor li se asociaza viteze initiale, corespunzatoare energiei cinetice pe care sistemul o poseda, energie determinata la randul ei de temperatura dorita pentru simularea respectiva. Aceasta temperatura este atinsa printr-o incalzire lenta a sistemului, plecand de la zero absolut, urmata de o perioada de atingere a echilibrului.

Intr-o simulare de dinamica moleculara mai intai se calculeaza fortele care actioneaza asupra fiecarui atom, prin variatia energiei potentiale ca urmare a variatiei pozitiei cu o cantitate infinitezimala dr_i; cu alte cuvinte, forta care se modifica in timp, este egala cu gradientul cu semn schimbat al energiei potentiale:

(III.1)

unde E reprezinta energia potentiala, iar r_i pozitia atomului i.

Energia poate fi calculata fie folosind mecanica moleculara, fie metodele mecanicii cuantice. Mecanica moleculara este limitata la acele procese care nu afecteaza drastic structura electronica a sistemului, cum ar fi formarea sau ruperea unei legaturi chimice.

Cunoscand fortele care actioneaza asupra atomilor si masa acestora m_i, se pot determina pozitiile lor de-a lungul unei serii de intervale de timp extrem de scurte de ordinul femtosecundelor (10^-15 s). Schimbarile structurale rezultate reprezinta asa numita traiectorie in timp a sistemului molecular.

Ecuatia (III.1) poate fi scrisa sub forma:

(III.2)

In practica traiectoria nu se obtine prin rezolvarea ecuatiei (III.2) datorita lipsei unei solutii analitice. De aceea, cunoscand fortele si masele se determina acceleratiile a_i. Vitezele si pozitiile atomilor se obtin prin integrarea acceleratiilor, respectiv a vitezelor in raport cu timpul, folosind relatiile:

(III.3)

Algoritmul folosit la integrare este cel al sariturii de broasca":

traiectoria intre doua stari este impartita in intervale de timp extrem de mici

odata calculata acceleratia la momentul t, a_i(t), se determina viteza la jumatatea intervalului:

(III.4)

noua pozitie a atomilor dupa intervalul , se calculeaza folosind viteza la mijlocul intervalului:

(III.5)

In felul acesta in timp ce pozitiile si acceleratiile sunt determinate dupa intervale de timp intregi , vitezele se determina la jumatati de intervale, sarindu-se peste valorile intregi. Din aceasta cauza valorile energiei totale sunt aproximative.

Utilizand intervale de timp extrem de scurte de 0,5 - 1 femtosecunde (fs), in cazul folosirii tuturor atomilor si de 1 - 2 fs pentru modelul atomului unit, se realizeaza o integrare adecvata pentru miscarile de inalta frecventa ale sistemului molecular (de exemplu vibratiile de valenta presupun perioade de ordinul picosecundelor).

Aceasta constituie, in acelasi timp si un dezavantaj al simularilor dinamicii moleculare, deoarece lungimea traiectoriei, calculata intr-un timp rezonabil, este cu cateva ordine de marime mai scurta decat orice proces chimic si decat foarte multe procese fizice, a caror durata poate fi de ordinul nanosecundelor sau chiar mai mare. De aceea prin dinamica moleculara pot fi studiate practic numai acele procese si proprietati, care dureaza mai putin decat durata simularii, cum ar fi fluctuatiile in energie sau schimbarea pozitiei atomilor, dar nu procese de durata, ca de exemplu plierea proteinelor.

In lipsa unui termostat neexistand nici o sursa externa , sau rezervor de energie, energia totala a sistemului este constanta. Intr-un ansamblu microcanonic (numarul particulelor, volumul si energia constante), teoretic, schimbarile in energia cinetica ar trebui sa fie de aceeasi marime, dar de semn opus cu cele ale energiei potentiale, astfel incat o simulare de dinamica moleculara sa poata fi facuta in conditiile de conservare a energiei; practic conservarea energiei este dificil de atins datorita limitarilor atat de capacitate cat si de viteza ale tehnicii de calcul. Aceste limitari impun o precizie finita si un calcul aproximativ al energiei, datorita in primul rand trunchierii la o distanta finita a interactiilor atomilor nelegati.

O alta posibilitate este de a efectua simularea dinamicii moleculare in conditii de temperatura constanta prin adaugarea unui termen cu ajutorul caruia se ajusteaza vitezele atomilor mentinand sistemul la temperatura dorita; in timpul unei astfel de simulari vitezele sunt ajustate prin scalare cu un factor dupa fiecare interval de timp ; sistemul este cuplat cu o baie de incalzire la temperatura de simulare T₀ si cu un timp de relaxare . In aceste conditii factorul de scalare este:

(III.6)

Pentru o traiectorie stabila se utilizeaza o constanta de relaxare mai mare decat 0,1 ps; o constanta de 0,01 ps este prea mica si creeaza distorsiuni in simulare.

III.2 Etapele simularii unei dinamici moleculare

O simulare a dinamicii moleculare are in general trei perioade distincte:

incalzirea, simularea si racirea. Pentru a masura proprietatile sistemului molecular la echilibru perioada de simulare se imparte in perioada de atingere a echilibrului si perioada de colectare de date.

III.2.1 Conditiile initiale si perioada de incalzire

Pentru a incepe o simulare de dinamica moleculara este necesara o configuratie initiala a sistemului molecular. Aceasta trebuie sa fie o structura moleculara optimizata cu un criteriu de convergenta adecvat: pentru temperatura camerei se efectueaza o minimizare folosind o valoare mica a gradientului ; pentru temperaturi mai mici gradientul trebuie sa fie si mai mic. Daca molecula, sau anumiti atomi sunt supusi la tensiuni mari, in timpul simularii se incearca o eliberare a acestora chiar din primii pasi, conducand la o explozie" a moleculei, sau la o distorsiune a structurii sale.

De cele mai multe ori in cazul simularilor efectuate pe biomolecule, structurile de plecare sunt cele obtinute prin difractia de raze X sau NMR furnizate de bancile de date precum Brookhaven Protein Databank (https: www.resb.org/pdb

Urmeaza un proces de atribuire a unor valori pentru vitezele atomilor care sa corespunda temperaturii de plecare; acestea sunt obtinute printr-o distributie Gaussiana a vitezelor, provenita de la un generator de numere aleatorii, tocmai pentru a evita introducerea unor miscari corelate inca de la inceputul simularii. Aceste viteze sunt ajustate (scalate) pentru fiecare atom astfel incat energia cinetica totala sa fie egala cu , unde k este constanta lui Boltzmann, iar T este temperatura initiala; dupa 100 - 300 de pasi de simulare (intervale de timp), distributia vitezelor devine o distributie Boltzmann.

Dupa aceasta initializare urmeaza incalzirea pana la temperatura de simulare; aceasta incalzire trebuie facuta lent, cu pasi mici, pentru a permite sistemului sa se apropie de echilibru in timpul fiecarei trepte de incalzire, astfel incat sa atinga rapid starea de echilibrul la atingerea temperaturii de simulare.

III.2.2 Echilibrul si colectarea de date

Atingerea starii de echilibru este necesara pentru a evita introducerea vreunui artefact in timpul etapei de incalzire si a asigura faptul ca traiectoria simulata va conduce la proprietatile starii de echilibru.

Perioada necesara stabilirii echilibrului depinde atat de proprietatile dorite cat si de natura sistemului molecular. Sunt unele proprietati care necesita perioade in jur de 100 ps, in timp ce altele devin stabile dupa 10 - 20 ps.

In timpul acestei perioade sunt corectate si vitezele atomilor, care in perioada de incalzire nu simuleaza in mod necesar miscarea reala a sistemului molecular, vitezele fiind date de o distributie arbitrara a valorilor corespunzatoare unei anumite temperaturi si fortelor care actioneaza asupra atomilor intr-o structura partial minimizata energetic. Pentru generarea unor viteze caracteristice echilibrului la temperatura dorita, in aceasta etapa, atomii sunt lasati sa interactioneze intre ei prin intermediul ecuatiilor de miscare.

Pentru o simulare izoterma, vitezele sunt scalate cu factorul , dat de ecuatia (III.6), in fiecare pas al simularii.

Pentru a stabili daca sistemul a atins starea de echilibru, se monitorizeaza fluctuatiile anumitor proprietati cum ar fi temperatura, energia potentiala, energia totala, numarul de legaturi de hidrogen sau frecventa de aparitie a anumitor conformeri; atingerea echilibrului este indicata de constanta in timp a acestor proprietati, atata timp cat nu au loc schimbari conformationale majore.

O grija deosebita trebuie avuta in cazul simularilor in solutie, cand apar fluctuatii chiar dupa zeci de picosecunde. Astfel solventul poate prezenta zone "fierbinti", cu un continut ridicat in energie, atingerea starii de echilibru pentru un sistem molecular complex, cum ar fi cel al unei proteine solvatate, fiind insotita de schimbarea conformatiei acesteia. Pentru a evita astfel de situatii este necesara o echilibrare initiala numai a moleculelor de apa (prin selectarea lor), ceea ce conduce la o relaxare a acestora, lasand neschimbata structura proteinei. Numai dupa aceea se echilibreaza intreg sistemul.

O alta situatie intalnita in special in cazul moleculelor mici de solvent, cum ar fi cele de apa, este aceea a cresterii temperaturii solventului in timp ce temperatura solutului scade. Aceasta se intampla datorita achizitiei de catre moleculele mici a unor cantitati mari de energie sub forma de energie de translatie si rotatie, in timp ce o macromolecula retine cea mai mare parte a energiei cinetice sub forma numeroaselor moduri de vibratie. Desi semnalata, o astfel de problema nu a putut fi solutionata in cadrul dinamicii moleculare.

Dupa atingerea echilibrului traiectoria poate fi stabila pe durata a mii de intervale , ceea ce permite achizitia de date sub forma de instantanee video (snapshots) sau fisiere de date privind conformatia (lungimi de legatura, unghiuri de valenta si unghiuri diedre) sau proprietati moleculare ce pot fi utilizate in calcule termodinamice.

Astfel, calculele dinamicii moleculare din pachetul de programe HyperChem mediaza si salveaza automat intr-un fisier cu extensia .csv urmatoarele valori:

Energia cinetica (E_C)

Energia potentiala (E_P)

Energia totala (E_T)

Temperatura (T)

Orice coordonata interna definita si denumita prin selectie: lungime

de legatura, unghi de valenta, unghi diedru

Abaterea medie patratica D_X (RMS-ul abaterilor):

(III.7)

unde X poate fi E_C, E_P, E_T,T, sau orice coordonata interna, n numarul

inregistrarilor, iar si reprezinta mediile marimi X si .

In general cu cat perioada de simulare este mai mare, cu atat datele privind proprietatile moleculare vor fi mai de incredere.

Se poate face medierea rezultatelor obtinute prin simulari repetate asupra aceluiasi sistem molecular, folosind conditii initiale diferite (durata etapelor), sau schimband vitezele si distributia lor. Schimbarea distributiei se poate face prin modificarea numarului si a marimii pasilor de incalzire, in timp ce schimbarea directiilor vitezelor initiale se face prin modificarea coordonatele atomilor in urma rotirii moleculei in sistemul de coordonate molecular.

III.2.3 Racirea lenta (simulated annealing)

Racirea lenta a sistemului molecular urmeaza etapei de incalzire si de atingere a starii de echilibru la temperatura aleasa pentru simulare, avand rolul de a indeparta tensiuni datorita incalzirii; acest proces este asemanator celui folosit pentru indepartare tensiunilor din sticla incalzita la o temperatura ridicata in timpul prelucrarii. Se foloseste in special cand se fac simulari la temperaturi inalte, de exemplu 600 K, urmate de racirea lenta la 300 K

Acest proces fie conduce starile conformationale de energie inalta catre conformatii stabile, fie ca prin depasirea barierelor de energie potentiala forteaza molecula sa treaca intr-o conformatie de energie mai coborata decat cea care s-ar obtine doar prin optimizarea geometriei prin mecanica moleculara.