SIMULAÇÃO COMPUTACIONAL DE BIOMOLÉCULAS - DINÂMICA MOLECULAR

Física computacional aplicada à Biologia – uma nova ferramenta

Neila Cristina Fonseca Machado

A simulação computacional é uma ferramenta que possibilita analisar teoricamente a natureza o que pode reduzir custo e tempo de experimentos. Com o aumento na capacidade de armazenamento de informações disponíveis nos computadores atuais e com a melhoria das metodologias científicas, as simulações são cada vez mais realísticas, pois agregam cada vez mais os elementos caracterizadores dos sistemas estudados. A escolha de determinada ferramenta computacional depende especialmente das dimensões do objeto investigado.

No estudo de macromoléculas biológicas, a Dinâmica Molecular (DM) é uma técnica computacional das mais versáteis, que se tornou uma ferramenta padrão para investigação de biomoléculas, onde se verifica uma alta qualidade na descrição das propriedades relevantes, como, estruturas cristalinas de proteínas, difusão em membranas etc. Com ela é possível predizer diversas propriedades, contribuindo para gerar resultados rápidos, respaldando experimentos, ou seja, complementando aquilo que não é entendido (acontece, mas não se vê os detalhes). A DM permite estudos de grandes sistemas, ou seja, um grande número de moléculas interagindo, o que representa muito bem os sistemas biológicos [1-4]. A seguir é definido o conceito de dinâmica molecular.

1 Dinâmica Molecular

Fundamentada nos princípios da mecânica clássica, a DM calcula o comportamento dinâmico microscópico de um sistema composto por átomos, em função do tempo [2-4], assim, propriedades macroscópicas de interesse são obtidas a partir da aplicação das mecânicas clássica e estatística. A última tem por função calcular propriedades macroscópicas observáveis (pressão, temperatura, energia interna, volume, entropia, energia livre, etc.) baseada em modelos microscópicos de átomos e moléculas.

Na mecânica molecular, um sistema de moléculas é tratado como uma coleção de átomos, descritos por forças newtonianas, onde os átomos são mantidos unidos por forças harmônicas ou elásticas [1].

A partir de um "campo de força" [5-9] empírico (conjunto completo de potenciais de interação entre partículas) conhecido como uma função de energia potencial é possível se calcular a energia potencial total do sistema V(r), com base na estrutura tridimensional do sistema. Um típico campo de força é representado pela Equação 1:

Equação 1

	V (r) =  ∑ Vl + ∑ Vɵ + ∑ Vɸ + ∑ VvdW + ∑ Velet

 

onde, V(r) descreve a soma de vários termos de energia, incluindo os termos para os átomos ligados (comprimentos e ângulos de ligação e ângulos diedros), e os termos para os átomos não ligados (interação de van der Waals e interações de Coulomb), sendo:

Vl é a energia de estiramento da ligação em relação ao seu valor de equilíbrio,

V_Ɵ é a energia de deformação do angulo em relação a seu valor de equilíbrio,

Vɸ é a energia devido à torção em torno de uma ligação,

V vdW representa a energia das interações de van der Waals e

V elet representa as energias de atração ou repulsão eletrostática entre duas cargas. [1]

Assim como a equação anterior existem várias outras que podem ser inseridas em softwares para fazer o modelamento de interação entre pares de átomos ligados e não ligados.

1.1 Interações entre pares de átomos ligados

Para este tipo de interações, a Lei de Hooke (Equações 2 e 3) da mecânica clássica é utilizada como a forma padrão para representar os potenciais harmônicos devido às oscilações dos comprimentos de ângulos de ligação em relação aos seus valores de equilíbrio [1].

Equação 2

Equação 3

onde: l comprimento de ligação, Ɵ ângulo de ligação, l₀  e Ɵ₀ respectivos valores de equilíbrio e k_l e k_Ɵ são as constantes de força pra restituição ao respectivos valores de equilíbrio.

Já o potencial de energia para torções, é representado pela equação 4, onde:

Equação 4

Vn é a barreira de energia para a torção, n é o número de máximos ou mínimos de energia em uma torção completa e ɸ é o angulo diedro que pode gerar um ponto de mínimo ou máximo na posição ɸ=0 [5]. n geralmente não excede 3 e dependerá do tipo de torção considerada.

Existem outros campos de força que levam em consideração, um quarto potencial harmônico, denominado “potencial torcional impróprio", usado para descrever oscilações e manter a estrutura tridimensional em um conjunto de 4 átomos [1].

Existem interações entre pares não ligados que serão explicadas a seguir.

1.2 Interações entre pares de átomos não ligados

Pares de átomos que não são ligados covalentemente (ij) são descritos por outros potenciais. Estes potenciais são compostos pelos termos de van der Waals representado pelos potenciais de Lennard-Jones (potencial 6-12) na Equação 5 e pelo potencial eletrostático, representado pelos potenciais de Coulomb na Equação 6.

Equação 5

Equação 6

A equação 5, descreve o potencial de Lennard-Jones, Ɛij trata-se da profundidade do potencial entre a barreira atrativa e repulsiva. σij é a distância finita na qual o potencial inter-partícula é zero. Estes parâmetros são ajustados experimentalmente ou por cálculos teóricos [1].

Na equação 6 das interações eletrostáticas, q_ij é a magnitude das cargas pontuais de cada átomo, r_ij corresponde a distância entre as cargas e Ɛ₀ a permissividade do espaço livre e Ɛr a constante dielétrica relativa ao meio.

Desenvolvidos independentemente e acoplando parâmetros específicos, os campos de força podem incluir outros termos, porém sua confiabilidade depende, em grande parte, do sistema a ser estudado e das propriedades a ser investigadas e, deste modo, a escolha do campo de força é primordial na simulação da DM [1].

2 Simulação da Dinâmica Molecular

Para um sistema atômico simples, a simulação da DM [1-4, 10] consiste na resolução numérica, passo a passo, da Equação de movimento (Equação 7).

Equação 7

Nesta equação, Fi descreve a força que atua sobre cada partícula do sistema em dado instante de tempo t. a_i trata-se da aceleração do átomo i de massa m_i.

Definido o campo de força é possível calcular as forças atuantes sobre cada átomo, através do cálculo da derivada primeira da energia potencial (obtida do campo de força selecionado) em relação às posições desses átomos.

Então, com a Equação 8 pode-se gerar diretamente a aceleração de cada partícula. Integrando as equações do movimento, obtém-se as velocidades cuja integral, fornece a mudança de posição do átomo [1].

Equação 8

Este procedimento aplicado sucessivamente gera uma trajetória, onde temos um conjunto de posições e velocidades de cada partícula ao longo do tempo, ou seja, a definição da Dinâmica Molecular [1].

2.1 Equações de Movimento – Integração

Para a integração das equações de movimento, utilizam-se algorítimos baseados nos métodos das diferenças finitas, ou seja, a integração é dividida em pequenos intervalos de tempo Δt (passos de integração) [11] permitindo simular os movimentos de maior frequência do sistema.

Um dos algorítimos mais utilizados é o Verlet [11] que utiliza as posições e aceleração dos átomos no tempo t e as posições do passo anterior, r (t-Δt), para estabelecer as novas posições no tempo, t+Δt, com base na equação 9:

Equação 9

Neste ponto do processo, tomamos conhecimento das principais equações que regem a DM, como lidar com pares de átomos ligados e os pares de átomos não ligados e como a movimentação destes átomos passo a passo é integrada pela equação de movimento. Agora vamos conhecer um pouco mais sobre como definir características especificas para nosso sistema (pressão, temperatura, volume, etc.). O que desejamos manter constante durante a simulação, ou o que desejamos variar.  Esta definição é feita através de um arquivo denominado emsemble estatístico, detalhado abaixo.

2.2 Emsemble Estatístico

Emsemble é um conjunto de propriedades e configurações mantidas constantes durante a integração das equações de movimento, representando o estado do sistema [1].

Existem várias opções de emsembles, nas quais é possível se controlar a Temperatura (T), pressão (P) ou a energia total (E) separadamente, por exemplo:

• NVE ou micro canônico: neste emsemble o número de partículas, volume e energia total, permanecem constantes durante a simulação.
• NVT ou canônico: número de partículas, volume e temperatura, constantes durante a simulação.
• NPT ou isotérmico-isobárico: número de partículas, pressão e temperatura, constantes.
•  NPH ou isobárico-isoentapico: número de partículas, pressão e entalpia constantes.
•  µVT ou grand-canônico: potencial químico, volume e temperatura constantes.

Por traz de cada propriedade dos emsembles, há vários processos de controle, por exemplo, para a temperatura, utilizam-se termostatos, para a pressão, barostatos. A seguir um breve descrição dos processos mais utilizados.

2.3 Controle da Temperatura (T)

O controle da temperatura durante o processo da DM é obtido através de um processo conhecido como termostato. Dentre os mais utilizados, o termostato de Berendsen [12], é capaz de manter a temperatura, através do acoplamento do sistema a um banho térmico com temperatura T₀ fixa. Re-escalonando as velocidades durante os passos de integração da simulação, obtém um ajuste da energia cinética do sistema a fim de se alcançar a temperatura escolhida. [1]

2.4 Controle da Pressão (P)

Para o controle da pressão, existem 2 barostatos muito utilizados, o barostato de Berendsen [12] e o barostato Parrinello-Rahaman [13], na qual as variáveis T são substituídas por P e as velocidades substituídas pelas as respectivas coordenadas atômicas [1].

Definido os parâmetros de controle, torna-se necessário definir o ponto inicial do sistema, ou seja, definir suas coordenadas cartesianas iniciais para então otimizar a geometria das moléculas (minimizar o sistema) e iniciar sua trajetória.

2.5 Condições iniciais

Para se iniciar uma simulação de DM, o primeiro passo é especificar as posições iniciais dos átomos que compõem o sistema. Deve-se, situar as partículas como uma rede cristalina, evitando deste modo sobreposições indesejada. Com isso criamos uma caixa com diferentes geometrias (dodecahedro, triclínica, cúbica, etc.) que se adeque melhor ao sistema em questão. Utilizam-se também condições de contorno periódicas, que visa minimizar possíveis efeitos de superfície indesejáveis [4] devido à discrepância na quantidade de partículas na superfície do sistema, quando se compara sistemas simulados e sistemas macroscópicos reais. Esta técnica permite que os átomos que estão delimitados dentro de uma caixa, sejam replicados em todas as direções do espaço, assim, quando um átomo se move na caixa original, sua imagem periódica, em uma das caixas imagens, move-se da mesma maneira, caso um átomo saia da caixa original, sua imagem entra pela face oposta da caixa imagem, com a mesma velocidade, mantendo assim o número total de átomos da caixa original [1].

2.6 Minimização do sistema

Após a definição das condições iniciais, o sistema deve ser minimizado visando eliminar maus contatos entre os átomos. O termo minimização, também conhecido como otimização da geometria, é uma técnica que busca encontrar um conjunto de coordenadas que minimizam a energia potencial do sistema. Um sistema minimizado tem pequenas forças sobre seus átomos, e serve como ponto de partida para se iniciar simulações de DM [1].

Para essa minimização, caminha-se sobre a superfície de potencial em direção decrescente da energia, assim, o sistema é levado a um mínimo de energia. Ajustando-se as posições atômicas, o processo relaxa distorções nas ligações químicas, ângulos e contatos de van der Waals. Dentre os vários algorítimos de minimização existentes, utilizamos o método steepest descent que utiliza a derivada primeira para determinar a direção do mínimo de energia, sendo robusta para minimização de estruturas que estão distantes de um ponto de mínimo [1].

Realizada a minimização, o sistema é aquecido gradualmente para a temperatura determinada T_0, e através da distribuição de Maxwell-Boltzmann, atribui-se velocidades iniciais de todas as partículas.

A distribuição de Maxwell-Boltzmann parte da premissa que a temperatura de qualquer sistema físico é o resultado da movimentação e colisões das partículas que compõem esse sistema, porém uma grande fração destas partículas, em um determinado intervalo de velocidades, as taxas de variação de deslocamento são quase constantes e assim, é possível calcular a fração de partículas (átomos) que estão num dado nível de energia (temperatura) em relação velocidade destas partículas.

Uma vez alcançado o equilíbrio termodinâmico, é possível então, gerar as trajetórias da DM calculando-se diversas propriedades do sistema em foco.

3 Abordagens da Dinâmica Molecular

Na DM dois tipos de abordagens de simulação podem ser tratados: a simulação em nível atomístico, conhecida como Fine-grained (FG) [7] e a simulação por agrupamento de átomos, conhecida como Coarse-grained (CG) [8, 9, 14]. Ambas as abordagens podem ser implementadas através de pacotes computacionais, como por exemplo, o Gromacs [15-20], usando o sistema operacional Linux.

A representação FG é uma representação onde cada átomo é considerado nos processos de interações (Fig.1). Os estudos usando representações FG permitem análises mais refinadas das interações entre moléculas, porém, com uma maior demanda computacional e alcançando pequenas escalas de tempo de simulação [8,9], o que às vezes, se torna incompatível para observação de processos biológicos.

Fig.1. Molécula do Ácido Ascórbico. Representação FG onde todos os átomos são representados no sistema. Em vermelho os átomos de oxigênio, carbono representado pela cor ciano e os átomos de hidrogênio em branco. Fonte: Machado, N. 2014.

Em estudos de auto-associação de fosfolipídios em solução aquosa, para se formar bicamadas, micelas ou lipossomos, torna-se computacionalmente demorado simular modelos utilizando a abordagem FG, assim, uma abstração do sistema atomístico por uma descrição CG, viabiliza a extensão da simulação para um maior tempo computacional o que é indicado para a observação da formação desses agregados, inclusive de difusão de componentes.

3.1 Coarse-grained

Na metodologia CG é realizado um mapeamento nas moléculas de interesse, agrupando-se em média, quatro átomos pesados ​​e seus hidrogênios associados. Este agrupamento é representado por um único centro de interação, em uma única esfera (bead). Esse mapeamento é conhecido como mapeamento quatro para um, e com isso, é possível equilibrar eficiência computacional e representatividade química. A base para parametrizar e descrever os beads é um conjunto de parâmetros descritos em um arquivo chamado de campo de força, neste caso para CG, o campo de força Martini [14].

Na figura 2 temos exemplos de moléculas com o mapeamento CG.

Fig. 2. Representação utilizando o campo de força Martini (CG). Em A e B, uma única esfera (bead) representa 4 moléculas de água, em C, um fosfolipídio, em D, fragmento de um polissacarídeo. Todas as esferas em CG estão representadas como esferas ciano-transparentes que estão recobrindo a estrutura atômica. Fonte: Marrink, S. J. and Tieleman, D. P. 2013 [14].

A fim de manter o modelo simples, apenas quatro principais tipos de sítios de interação são definidos, com base no caráter mais ou menos polar do agrupamento em questão [8, 9, 14]:

1 - Polar (P): grupos de átomos que são solúveis em água

2 - Não-polar (N): são usados em grupos mistos com parte polar e parte apolar

3 - Apolar (C): hidrofóbicos

4 - Carregada (Q): grupos ionizados.

Dentro de cada tipo, existem subtipos que são diferenciados por uma letra, que indica as capacidades de ligação de hidrogênio ou por um número, indicando o grau de polaridade (a partir de 1 = baixa polaridade, até 5 = alta polaridade), num total de 18 tipos de esferas, que permitem eficiente representatividade química.

Os subtipos permitem um fino ajuste das interações com base na natureza química dos átomos. A figura 3 representa parte do processo de agrupamento dos átomos, transformando uma molécula da metodologia FG para CG [8, 9, 14].

Fig. 3. Transformação de uma molécula de FG para CG. Neste exemplo uma molécula com 67 átomos em FG é convertida para a representação CG com 8 esferas. Fonte: Machado, N. 2014.

O Campo de Força Martini tem suas bases em experimentação química e modelos atomísticos, contudo sua filosofia é diferente. Seu foco não está na reprodução acurada de detalhes estruturais de um sistema em particular, o foco está em oferecer uma ampla gama de aplicações (Fig. 4), sem a necessidade de re-parametrização dos modelos das moléculas a cada nova aplicação.

Inicialmente o campo de força Martini foi elaborado para tratar simulações de sistemas lipídicos, assim, no geral, as aplicações se concentram em torno de membranas lipídicas, porém, há um leque de aplicações [14]:

• Membranas lipídicas: caracterização de propriedades, formação de poros através de ligação de peptídeos ativos, fusão de membranas, compressão e expansão de monocamada;
• Interação de nanopartículas com membranas;
• Polimorfismo de lípidios;
• Lipossomos: estrutura e dinâmica;
• Proteínas: alterações conformacionais, interação proteína-lípideos, oligomerização de proteínas de membrana, auto-montagem de peptídeos e proteínas solúveis;
• Design de drogas e sistemas de entrega de genes;
• Auto-montagem de surfactantes;
• Caracterização de sistemas baseados em carboidratos;
• Polímeros: estrutura e dinâmica;

Fig. 4. Simulação em CG, interação de nanopartículas com lipídios. Em A, fulereno em membrana lipídica [21]. Em B, grafeno entre camadas lipídicas [22], em C, nanotubos de carbono envolto por surfactantes [23]. Figuras extraídas do artigo Marrink, J. and Tieleman D. P., 2013 [14].

Comparados a modelos atomísticos em FG, os modelos CG são mais rápidos em até 4 ordens de grandeza temporal e por isso estão se tornando cada vez mais populares para sistemas de lipídios e surfactantes. Apesar da natureza simplista do modelo, revela-se versátil em suas aplicações e preciso em suas previsões, combinando velocidade, versatilidade e especificidade química [8, 9, 14]. O objetivo da simulação em CG é obter uma simulação mais realista de grandes sistemas complexos, por um longo período de tempo, através da redução dos sítios de interação, eliminação de processos dinâmicos muito rápidos, mantendo assim a física essencial do sistema.

4 Bibliografia

1. NAMBA, A. M. et al. Dinâmica molecular: teoria e aplicações em planejamento de fármacos. Eclética Química, v. 33, n. 4, p. 13-24, 2008.

2. HANSSON, Tomas; OOSTENBRINK, Chris; VAN GUNSTEREN, WilfredF. Molecular dynamics simulations. Current opinion in structural biology, v. 12, n. 2, p. 190-196, 2002.

3. KARPLUS, Martin; MCCAMMON, J. Andrew. Molecular dynamics simulations of biomolecules. Nature Structural & Molecular Biology, v. 9, n. 9, p. 646-652, 2002.

4. FRENKEL, Daan; SMIT, Berend. Understanding molecular simulation: from algorithms to applications. Computational sciences series, v. 1, p. 1-638, 2002.

5. VAN GUNSTEREN, W. F.; BERENDSEN, H. J. C. Groningen molecular simulation (GROMOS) library manual. Biomos, Groningen, v. 24, n. 682704, p. 13, 1987.

6. VAN GUNSTEREN, Wilfred F. et al. Biomolecular simulation: The {GROMOS96} manual and user guide. 1996.

7. OOSTENBRINK, Chris et al. A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force‐field parameter sets 53A5 and 53A6. Journal of computational chemistry, v. 25, n. 13, p. 1656-1676, 2004.

8. MARRINK, Siewert J.; DE VRIES, Alex H.; MARK, Alan E. Coarse grained model for semiquantitative lipid simulations. The Journal of Physical Chemistry B, v. 108, n. 2, p. 750-760, 2004.

9. MARRINK, Siewert J. et al. The MARTINI force field: coarse grained model for biomolecular simulations. The Journal of Physical Chemistry B, v. 111, n. 27, p. 7812-7824, 2007.

10. VAN GUNSTEREN, Wilfred F.; BERENDSEN, Herman JC. Computer simulation of molecular dynamics: Methodology, applications, and perspectives in chemistry. Angewandte Chemie International Edition in English, v. 29, n. 9, p. 992-1023, 1990.

11. VERLET, Loup. Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Physical review, v. 159, n. 1, p. 98, 1967.

12. BERENDSEN, Herman JC et al. Molecular dynamics with coupling to an external bath. The Journal of chemical physics, v. 81, n. 8, p. 3684-3690, 1984.

13. PARRINELLO, Michele; RAHMAN, Aneesur. Polymorphic transitions in single crystals: A new molecular dynamics method. Journal of Applied physics, v. 52, n. 12, p. 7182-7190, 1981.

14. MARRINK, Siewert J.; TIELEMAN, D. Peter. Perspective on the Martini model. Chemical Society Reviews, v. 42, n. 16, p. 6801-6822, 2013.

15. Bekker H. et al.  Gromacs: A parallel computer for molecular dynamics simulations. In Physics computing. v. 92, p. 252-256, 1993.

16. BERENDSEN, Herman JC; VAN DER SPOEL, David; VAN DRUNEN, Rudi. GROMACS: a message-passing parallel molecular dynamics implementation. Computer Physics Communications, v. 91, n. 1, p. 43-56, 1995.

17. LINDAHL, Erik; HESS, Berk; VAN DER SPOEL, David. GROMACS 3.0: a package for molecular simulation and trajectory analysis. Molecular modeling annual, v. 7, n. 8, p. 306-317, 2001.

18. VAN DER SPOEL, David et al. GROMACS: fast, flexible, and free. Journal of computational chemistry, v. 26, n. 16, p. 1701-1718, 2005.

19. HESS, Berk et al. GROMACS 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation. Journal of chemical theory and computation, v. 4, n. 3, p. 435-447, 2008.

20. PRONK, Sander et al. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics, v. 29, n. 7, p. 845-854, 2013.

21. WONG-EKKABUT, Jirasak et al. Computer simulation study of fullerene translocation through lipid membranes. Nature Nanotechnology, v. 3, n. 6, p. 363-368, 2008.

22. TITOV, Alexey V.; KRÁL, Petr; PEARSON, Ryan. Sandwiched graphene− membrane superstructures. ACS nano, v. 4, n. 1, p. 229-234, 2009.

23. WALLACE, E. Jayne; SANSOM, Mark SP. Carbon nanotube/detergent interactions via coarse-grained molecular dynamics. Nano letters, v. 7, n. 7, p. 1923-1928, 2007.