The black sheep

As simula¸c˜oes de DM s˜ao t´ecnicas para determinar o movimento das part´ıculas de um sistemas, para o qual o potencial de intera¸c˜ao entre as part´ıculas (fun¸c˜ao de energia potencial) j´a ´e conhecido. Esse de descri¸c˜ao das posi¸c˜oes instantˆaneas e velocidades instantˆaneas das part´ıculas, pode ser dividido nas seguintes est´agios: estabelecer a confi- gura¸c˜ao inicial para o sistema, a computa¸c˜ao das for¸cas exercidas sobre cada part´ıcula em raz˜ao das intera¸c˜oes intermoleculares, movimenta¸c˜ao das part´ıculas ou equilibra¸c˜ao de fase, ensembles e controles de simula¸c˜ao, e finalmente, an´alise das configura¸c˜oes (trajet´orias) obtidas.

A configura¸c˜ao inicial do sistema, tem como frequente metodologia a cria¸c˜ao de uma ”caixa de simula¸c˜ao”, que geralmente ´e uma caixa c´ubica, onde todas as mol´eculas do sistema se encontram. Uma boa pr´atica das condi¸c˜oes iniciais ´e que o sistema n˜ao contenha intera¸c˜oes com altas energias, que podem por sua vez, causar instabilidades na simula¸c˜ao. Se o sistema possuir um arranjo experimental dispon´ıvel (por exemplo, estrutura de raio X) esta ser´a a configura¸c˜ao inicial suficiente para o in´ıcio da simula¸c˜ao. A estrutura cristalogr´afica mais comum ´e a rede c´ubica de face centrada, mas podem ser escolhidas outras geometrias adequadas. As dimens˜oes da caixa de simula¸c˜ao devem ser escolhidas de maneira que reproduza a densidade num´erica, N/V do sistema real e suas condi¸c˜oes termodinˆamicas. A caracteriza¸c˜ao do sistema considera a reprodu¸c˜ao das condi¸c˜oes fisiol´ogicas que o aproximam do modelo real, e os sistemas s˜ao tratados, de forma que, para a maior parte das simula¸c˜oes o ambiente escolhido ´e um solvente aquoso, podendo ser uma solvata¸c˜ao impl´ıcita ou expl´ıcita. Para o caso de solvente expl´ıcito, existem alguns modelos expl´ıcitos para a ´agua, que incluem os modelos TIP3P, TIP4P (JORGENSEN et al., 1983), SCP e SCPC/E (ADCOCK; MCCAMMON, 2006), e os modelos mais conhecidos s˜ao consistentes com a aproxima¸c˜ao SHAKE (ANDERSEN, 1983) para manuten¸c˜ao da geometria molecular, e ainda para diminuir os efeitos de fronteiras, condi¸c˜oes de contorno s˜ao utilizadas (LEACH, 2001).

Para iniciar a simula¸c˜ao de DM, a solu¸c˜ao da equa¸c˜ao de movimento, ´e necess´ario fornecer n˜ao somente as posi¸c˜oes das part´ıculas, mas tamb´em as suas velocidades iniciais, j´a que as equa¸c˜oes diferencias s˜ao de segunda ordem.

Durante a simula¸c˜ao as for¸cas entre as part´ıculas decaem com a distˆancia entre estas, intera¸c˜oes Coulombianas, a medida mais aplicada ´e na caixa de simula¸c˜ao, substituir a for¸ca de intera¸c˜ao entre duas part´ıculas por uma for¸ca efetiva, soma das intera¸c˜oes entre todas as imagens peri´odicas dessas part´ıculas, pelo m´etodo de soma infinita. O m´etodo de soma de Edwald (YORK; YANG, 1994) oferece uma abordagem rigorosa para a avalia¸c˜ao das intera¸c˜oes eletrost´aticas, al´em do Edwald, existem m´etodos que s˜ao variantes deste, como o m´etodo Particle Mesh Ewald (PME), que analisa as for¸cas diferenciando analiticamente, diminuindo significantemente a demanda de mem´oria computacional. Pela eficiˆencia em velocidade de processamento e rigor, o PME est´a implementado em muitos pacotes de simula¸c˜ao (TOUKMAJI; BOARD, 1996; WANG; GAO; FANG, 2016).

2.5.2.1 Ensemble Estat´ıstico

Um sistema cl´assico com N part´ıculas interagindo em equil´ıbrio, a uma temperatura constante T, no ensemble canˆonico ´e reproduzido pela distribui¸c˜ao de Boltzmann. A m´edia da quantidade G(~rN_{), ´e unicamente uma fun¸c˜ao da posi¸c˜ao ~r de todas as N part´ıculas e ´e} expressada como uma integral de espa¸co de fase incluindo a energia potencial U (~rN_):

hGi = R G(~r

N_{)exp(−βU(~r}N_))d~rN

R exp(−βU(~rN_))d~rN (75)

na equa¸c˜ao 75acima β = _kBT1 , onde kB ´e a constante de Boltzmann.

Se, depois de um tempo suficientemente grande, sob certas condi¸c˜oes, o sistema em equil´ıbrio explora todo o espa¸co de fase de microestados com a mesma energia, onde a densidade de probabilidade ´e n˜ao nula, ou seja todos os microestados acess´ıveis s˜ao igualmente prov´aveis nesse longo per´ıodo de tempo, esse sistema ´e erg´odico. Essa ´e a hip´otese erg´odica no qual a m´edia de uma amostragem ao longo do tempo, equivale a m´edia no ensemble que ´e definida como:

hGit = lim M →∞hGiM = limM →∞ 1 M M X µ=1 G(~rN(tµ)) (76)

onde hGiM ´e a m´edia sobre as medi¸c˜oes M, da quantidade G em todos os tempos {tµ}Mµ=1. Para uma boa amostragem, ou seja, um grande n´_{umero de medi¸c˜oes M, a m´edia hGi}M ´e considerada uma aproxima¸c˜ao conveniente para a m´edia do ensemble em equil´ıbrio (SALINAS, 1997).

No decorrer da simula¸c˜ao, integra¸c˜ao das equa¸c˜oes de Newton, s˜ao mantidas constantes uma cole¸c˜ao de configura¸c˜oes e propriedades que caracteriza o estado do sistema. Essa cole¸c˜ao de todos os prov´aveis sistemas que possuem estados microsc´opicos diferentes, mas participam de um ´unico estado macrosc´opico (termodinˆamico), chamamos de ensemble estat´ıstico. Existem v´arios ensemble com caracter´ısticas diferentes, nos quais distintos parˆametros permanecem constantes simultaneamente como o ensemble canˆonico (NVT): conjunto de sistemas nos quais o estado termodinˆamico ´e caracterizado por um n´umero fixo de ´atomos, N, fixado o volume, V, e fixada a temperatura, T; o ensemble isob´arico-isoent´alpico (NPH- n´umero de ´atomos, press˜ao e entalpia constantes), o ensemble gran-canˆonico (µVT- potencial qu´ımico, volume e temperatura constantes), o ensemble isob´arico-isot´ermico (NPT- n´umero de ´atomos, press˜ao e temperatura constantes), dentre outros(CHANDLER, 1987).

Grande parte das simula¸c˜oes de DM tem suas condi¸c˜oes inicias representadas pelo ensemble microcanˆonico (NVE), onde a energia total E ´e uma constante de movimento, e sua m´edia temporal s˜ao iguais a sua m´edia no ensemble. Contudo as experiˆencias laboratoriais (por exemplo, cristalografia de raios X) s˜ao realizadas com volume e temperatura contantes (ensemble canˆonico) ou `a press˜ao e temperatura constantes (ensemble isob´arico-isot´ermico), sendo uma pr´atica na DM reproduzir, sob condi¸c˜oes fisiol´ogicas, esses estados (SCHLICK, 2010). Inicialmente `a energia constante, a temperatura flutua como consequˆencia da interconvers˜ao espontˆanea das componentes cin´etica e potencial da energia total. A energia interna total pode ser escrita como a soma das contribui¸c˜oes cin´eticas (K ) e potencial (U ) envolvendo a soma sobre todos os pares de mol´eculas:

E = hK i + hU i (77)

A temperatura e press˜ao podem ser calculados via teorema do virial. Para o caso atˆomico, somando sobre o 3N termos, a energia cin´etica tem a forma:

h N X

i=1

|pi|2/mii = 2hκi = 3NkBT, (78) onde kB´e a constante de Boltzmann, mi e pi s˜ao, respectivamente, massa e momento do ´atomo I e N ´e o n´umero total de ´atomos. A equa¸c˜ao78acima fundamenta-se no princ´ıpio

da equiparti¸c˜ao e a fun¸c˜ao de temperatura cin´etica instantˆanea pode ser determinada como: T = 2K /3N kB = 1 3N kB N X i=1 |pi|2/mi. (79)

Sendo a temperatura T (t) no instante t, e as velocidades alteradas por um fator λ, as varia¸c˜oes de temperatura podem ser calculadas da seguinte maneira:

∆T = 1 2 X i=1 2mi(λυi) 2 NdfkB − 1 2 X i=1 2 miυ 2 i NdfkB (80) ∆T = (λ2− 1)T (t) λ =pT0/T (t)

a temperatura tem como controle simples, multiplicar as velocidade atˆomicas instantˆaneas pelo fator λ. Este m´etodo tem como variante o termostato de Berendsen (H ¨UNENBERGER, 2005), que modera a raz˜ao com a qual a temperatura do sistema chega ao valor determinado. Assumindo o ensemble NVT, para o qual a energia total n˜ao ´e conservada, o inicio da simula¸c˜ao equipara-se de fato a um banho t´ermico, fazendo desse ensemble indispens´avel para o estagio de termaliza¸c˜ao, quando tratamentos t´ermicos como o Simulated Annealing s˜ao utilizados (KARPLUS; PETSKO et al., 1990; ERPENBECK; WOOD, 1984).

No ensemble NPT o volume do sistema pode flutuar alterando as dimens˜oes da caixa de simula¸c˜ao e tamb´em pelo reescalonamento das posi¸c˜oes,o fator γ = pp(t)/p3

0 ´e assim, respons´avel por manter a taxa de flutua¸c˜ao equilibrada, de maneira similar como feito para a temperatura. A press˜ao pode ser calculada com a contribui¸c˜ao do teorema do virial de acordo com a equa¸c˜ao:

p(t) = N kBT V − 1 3V N X i=1 N X j>1 rij(t) · Fij(t) (81) onde rij(t) = rj(t)−ri(t) ´e o vetor posi¸c˜ao das mol´eculas i e j, e Fij(t) a for¸ca intermolecular entre i e j, no instante t. Existem alguns m´etodos em simula¸c˜oes de DM que fazem com que a press˜ao fique corretamente constante durante a simula¸c˜ao, como por exemplo, o m´etodo dinˆamica de Langevin (FELLER et al., 1995).

2.5.2.2 An´alise e armazenamento trajet´orias geradas na DM

O est´agio inicial da dinˆamica a partir da configura¸c˜ao inicial ´e a termaliza¸c˜ao, e deve ser inspecionado com muito cuidado at´e que a configura¸c˜ao inicial n˜ao mais persistir. Neste est´agio, como mencionado acima, a temperatura e a press˜ao devem ser escalonadas, e o banho t´ermico sistematicamente observado, at´e que o sistema se encontre no equil´ıbrio. A credibilidade dos resultados no est´agio seguinte, produ¸c˜ao, depende do monitoramento adequado do est´agio de equilibra¸c˜ao.

O est´agio de produ¸c˜ao ´e a reuni˜ao das trajet´orias (coordenadas ou velocidades atˆomicas) no tempo de simula¸c˜ao que sucede aquele do equil´ıbrio. A partir dessas trajet´orias ´e que ser˜ao calculadas propriedades f´ısico-qu´ımicas do sistema e m´edias sobre os ensembles das trajet´orias (SCHLICK, 2010).

Dentre v´arias propriedades que podem ser calculadas a partir das trajet´orias est˜ao as propriedades estruturais, como a fun¸c˜ao de distribui¸c˜ao de pares, g(r), que produz a probabilidade de se encontrar uma part´ıcula β em uma camada esf´erica de espessura dr, e sua distˆancia at´e a part´ıcula α vale r posicionada na origem. Considerando as coordenadas cartesianas retiradas das trajet´orias de simula¸c˜ao, a fun¸c˜ao radial de pares pode ser dada por: g(r) = ρ−2 * X i X j6=i δ(ri)δ(rj− r) + = V N2 * X i X j6=i δ(r − rij) + (82) Assim sendo, o n´umero de part´ıculas N (r) que localiza-se na faixa estreita ∆r a uma distˆancia r da origem, remete a distribui¸c˜ao radial dada pela seguinte equa¸c˜ao:

g(r) = N (r)

4πr2_N (83)

O deslocamento quadr´atico m´edio RMSD (do inglˆes, root mean square displacement) calcula a diferen¸ca estrutural do sistema em um determinado passo de simula¸c˜ao com rela¸c˜ao as configura¸c˜oes iniciais, ou seja, realiza uma compara¸c˜oes entre as coordenadas atˆomicas de estruturas proteicas, por exemplo, adquiridos das trajet´orias de DM, que tem como dois principais interesses:

Observar a evolu¸c˜ao temporal, dado um conjunto de ´atomos espec´ıficos ´e realizada uma m´edia de RMSD a cada passo de simula¸c˜ao, aferindo as mudan¸cas ocorridas durante a simula¸c˜ao, que s˜ao indicadas por flutua¸c˜oes em torno de uma posi¸c˜ao de equil´ıbrio. A

evolu¸c˜ao temporal de RMSD ´e geralmente visualizada atrav´es de um gr´afico de RMSD m´edio versus tempo, e o RMSD m´edio ´e calculado segundo a equa¸c˜ao:

hRMSDi = s

PN

j=1[rj(ti) − r0]2

N (84)

onde rj(ti) ´e a distˆancia do ´atomo j no passo de simula¸c˜ao i, r0 ´e a distˆancia de referˆencia e N ´e o n´umero total de ´atomos.

Analisar a m´edia do RMSD por ´atomo, que ´e encontrado gerando as m´edias dos RMSD, para cada ´atomo, em todos os passos da simula¸c˜ao, e tem como finalidade principal, por compara¸c˜ao, encontrar quais ´atomos ou regi˜oes do sistema flutuam mais ou menos. A m´edia do RMSD ´e encontrado utilizando a equa¸c˜ao:

hRMSDji = s

PM

i=1[rj(ti− r0)]2

M · ∆t (85)

na equa¸c˜ao acima M ´e o n´umero de configura¸c˜oes e ∆t o intervalo de tempo entre os passos de simula¸c˜ao, de maneira que multiplicando M por ∆t, tem-se o tempo de simula¸c˜ao total.