Dinâmica molecular: teoria e aplicações em planejamento de fármacos

Molecular dynamics: theory and applications in drug design

Dinâmica molecular: teoria e aplicações em planejamento de fármacos

Molecular dynamics: theory and applications in drug design

A. M. Namba^*; V. B. da Silva; C. H. T. P. da Silva

Departamento de Ciências Farmacêuticas, Faculdade de Ciências Farmacêuticas de Ribeirão Preto, Universidade de São Paulo, Av. do Café, s/n, Monte Alegre, 14040-903, Ribeirão Preto-SP, Brasil

RESUMO

Dinâmica Molecular (DM) é uma ferramenta computacional poderosa usada em Química Medicinal para o planejamento racional de fármacos. DM é uma extensão da Mecânica Molecular, onde o comportamento dinâmico de um sistema molecular é simulado através da integração numérica das equações de movimento. Esta técnica tem sido usada extensivamente para auxiliar e complementar o planejamento de novos ligantes de um alvo terapêutico, bem como estimar a sua potência. Este artigo enfoca a teoria básica da DM clássica e suas importantes aplicações no planejamento racional de potenciais compostos bioativos, particularmente compostos com atividade anti-HIV.

Palavras-chave: dinâmica molecular; campo de força; planejamento de fármacos; atividade anti-HIV.

ABSTRACT

Molecular dynamics (MD) is a powerful computational tool used for rational drug design in medicinal chemistry. MD is an extension of Molecular Mechanics, where the dynamic behavior of a molecular system is simulated through numerical integration of the equations of motion. This approach has been extensively used to aid and complement the design of novel ligands to a therapeutic target as well as to estimate their potency. Our review focuses on the basic theory of classical MD and their important applications in rational drug design of potential bioactive compounds, particularly anti-HIV compounds.

Keywords: molecular dynamics; force field; drug design; anti-HIV activity.

Introdução

Simulação de Dinâmica Molecular (DM) é uma das técnicas computacionais mais versáteis para o estudo de macromoléculas biológicas [1-4]. No planejamento racional de fármacos baseado em estrutura, as simulações de DM têm contribuído exaustivamente em diversos estágios do processo [4-6].

A metodologia da DM é fundamentada nos princípios da Mecânica Clássica e fornece informações sobre o comportamento dinâmico microscópico, dependente do tempo, dos átomos individuais que compõem o sistema [7,8]. Para se obter as propriedades macroscópicas de interesse a aplicação da mecânica estatística é requerida, a qual tem a função de calcular propriedades observáveis macroscópicas (pressão, energia interna, volume, temperatura, entropia, energia livre, etc), a partir de outras microscópicas.

Com base na Mecânica Molecular (MM), as moléculas são tratadas como uma coleção de átomos que pode ser descrita por forças newtonianas, ou seja, são tratadas como uma coleção de partículas mantidas unidas por forças harmônicas ou elásticas. Um conjunto completo dos potenciais de interação entre as partículas é referido como "campo de força" [9,10]. O campo de força empírico, tal como é conhecido como uma função energia potencial, permite que a energia potencial total do sistema, V(r), seja calculada como a partir da estrutura tridimensional (3D) do sistema. V(r) é descrito como a soma de vários termos de energia, incluindo os termos para átomos ligados (comprimentos e ângulos de ligação, ângulos diedros) e os termos para átomos nãoligados (interações de van der Waals e de Coulomb). Um típico campo de força é representado pela Equação 1:

onde, V_l é a energia de estiramento da ligação em relação a seu valor de equilíbrio (ou ideal), V_θé a energia de deformação do ângulo de ligação em relação a seu valor de equilíbrio, V_φ é a energia devido à torção em torno de uma ligação, V_νdW representa a energia das interações de van der Waals e V_elec representa as energias de atração ou repulsão eletrostática entre duas cargas.

A forma padrão para representar os potenciais harmônicos devido às oscilações dos comprimentos e ângulos de ligação com relação aos valores de equilíbrio é a utilização da Lei de Hooke, mostrada nas equações 2 e 3:

onde l e θ são os comprimentos e ângulos de ligação, respectivamente, l₀ e θ₀ são os correspondentes valores de equilíbrio e k_l e k_θas constantes de força para a restituição aos respectivos valores de equilíbrio. A aproximação harmônica é válida somente para pequenas distorções com relação aos valores de equilíbrio.

A forma funcional padrão para representar o potencial de energia para uma torção está mostrada na Equação 4:

onde V_n é 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, φ é o ângulo diedro e γ é o ângulo de fase (defasagem no ângulo diedro que pode gerar um ponto de mínimo ou de máximo na posição φ = 0) [11]. O valor de n dependerá do tipo de torção considerada e, geralmente, não excede o valor 3.

Em alguns campos de força, um quarto potencial harmônico é definido para manter a estrutura tridimensional em um conjunto de quatro átomos, descrevendo ainda oscilações da mesma. Esse potencial é denominado "potencial torcional impróprio".

As interações entre pares de átomos nãoligados covalentemente (i, j) são descritas por potenciais compostos pelos termos de van der Waals e eletrostático, representados respectivamente pelos potenciais de Lennard-Jones (potencial 6-12, Equação 5) e de Coulomb (Equação 6).

Na Equação 5, ε_ij é a profundidade do potencial entre a barreira atrativa e a repulsiva, e σ_ij é a distância finita na qual o potencial interpartícula é zero. Ambos são parâmetros ajustados experimentalmente ou por cálculos teóricos. No caso das interações eletrostáticas (Equação 6), q_i e q_j correspondem à magnitude das cargas pontuais de cada átomo, r_ij a distância entre as cargas, ε_o a permissividade do espaço livre e ε_ra constante dielétrica relativa do meio.

Os campos de forças existentes foram desenvolvidos de maneira independente e com todos os conjuntos de parâmetros específicos. Alguns incluem outros termos para descrever especificamente as ligações de hidrogênio ou para acoplar oscilações entre ângulos e comprimentos de ligação, com o objetivo de se obter uma melhor concordância com espectros vibracionais [9,10]. A confiabilidade dos resultados é baseada na elaboração de um campo de força com parâmetros bem definidos. A escolha do campo de força depende, em grande parte, do sistema a ser estudado e das propriedades que serão investigadas. No caso de sistemas biomoleculares, os campos de força mais utilizados são CHARMM [12], GROMOS [11,13], AMBER [14], OPLS, [15], CVFF [16], entre outros.

Dinâmica molecular clássica

A simulação de DM [17-25] consiste da solução numérica, passo a passo, da Equação de movimento, que pode ser descrita para um sistema atômico simples pela Equação 7:

onde F_i é a força que atua sobre cada partícula do sistema em um instante de tempo t, e a_i é a aceleração do átomo i de massa m_i. Uma vez definido o campo de força, é possível calcular as forças que atuam sobre cada átomo, calculando-se a derivada primeira da energia potencial, obtida do campo de força escolhido, em relação às posições desses átomos (Equação 8). A Equação 8 gera diretamente a aceleração da partícula. A partir desta, integrandose as equações de movimento, pode-se obter as velocidades, cuja integral, por sua vez, proporciona a mudança de posição do átomo. Com as novas posições e velocidades de cada partícula, obtêm-se as energias potencial e cinética do sistema. Aplicando-se sucessivamente esse procedimento, obtém-se o que se denomina de "trajetória", que nada mais é do que o conjunto de posições e velocidades de cada partícula ao longo do tempo.

Integração das equações de movimento

A integração das equações de movimento é realizada através de algoritmos baseados nos métodos das diferenças finitas, nos quais a integração é dividida em pequenos intervalos de tempo (passos de integração) [27], Δt, permitindo simular os movimentos de maior freqüência do sistema, que muitas vezes são as vibrações das ligações. Um dos métodos mais utilizados em dinâmica molecular para integrar as equações de movimento é o algoritmo de Verlet [17,27], que utiliza as posições e acelerações dos átomos no tempo t e as posições do passo anterior, r (t-Δt), para determinar as novas posições no tempo t+Δt, de acordo com a Equação 9:

Este algoritmo pode ser obtido através da expansão em série de Taylor, inicialmente para frente (Equação 10) e depois para trás (Equação 11). Somando-se as duas e isolando r (t+Δt), chega-se à Equação 9.

Além do algoritmo de Verlet, existem outros algoritmos de integração confiáveis como o algoritmo de Beeman [28] e o método leap-frog [29].

Ensemble estatístico

O conjunto de configurações e propriedades mantidas constantes durante a integração das equações de Newton e que representam o estado do sistema, é chamado de ensemble [17]. Uma das primeiras opções de ensemble é o microcanônico ou NVE ( número de partículas, volume e energia total constantes durante a simulação). Outras alternativas para esse ensemble foram desenvolvidas, nas quais se controla separadamente a temperatura, T, e a pressão, P, ao invés da energia total (E). Dentre esses encontram-se o grand-canônico ou µVT (potencial químico, volume e temperatura constantes), o canônico ou NVT (com número de partículas, volume e temperatura constantes), o isobárico-isoentapico ou NPH (com número de partículas, pressão e entalpia constantes) e o isotérmico-isobárico ou NPT (com número de partículas, pressão e temperatura constantes) [17,26].

Controle da temperatura e da pressão

Durante as simulações de DM o controle de T pode ser obtido através de um processo conhecido como termostato de Berendsen [30], que para manter a temperatura, o sistema é acoplado a um banho térmico com tempertaura T₀ fixa. Durante as simulações as velocidades são re-escalonadas a cada passo de integração, a fim de ajustar a energia cinética do sistema até a temperatura selecionada [30].

Outros métodos mais sofisticados como os termostatos de Nosé-Hoover [31,32] e Nosé-Poincaré [32] também são utilizados. Para o controle da pressão, as variáveis T são substituídas por P, e as velocidades dos átomos pelas coordenadas atômicas. Os barostatos mais utilizados são os de Berendsen [30] e Parrinello-Rahaman [33].

Condições iniciais, minimização de energia e equilibração

Para se realizar uma simulação de DM é importante especificar as posições iniciais das partículas que constituem o sistema. É conveniente situar inicialmente as partículas nas posições de uma rede cristalina, o que evita sobreposição indesejada entre elas, gerando uma caixa de simulação que pode adotar diferentes geometrias. É comprovado que os resultados das simulações são praticamente independentes do tipo de geometria utilizada para as caixas de simulação. As mais comumente usadas são a cúbica, a paralelepipédica e a octaédrica truncada.

Antes de se iniciar as simulações de DM, o sistema deve ser minimizado para eliminar maus contatos entre os átomos. A minimização de energia, também conhecida como otimização da geometria, é uma técnica que visa encontrar um conjunto de coordenadas que minimizam a energia potencial do sistema de interesse. O procedimento básico consiste em caminhar sobre a superfície de potencial na direção em que a energia decresce, de maneira que o sistema é levado a um mínimo de energia local próximo. Através de ajustes nas posições atômicas, o processo relaxa as distorções nas ligações químicas, nos ângulos entre ligações e nos contatos de van der Waals. [7,34]. O sistema minimizado possui forças pequenas sobre cada átomo e serve, portanto, como estrutura de partida para iniciar as simulações de DM.

Os algoritmos de minimização mais conhecidos são: o método steepest descent [3537], o método dos gradientes conjugados [38] e o método de Newton-Raphson [39,40]

O método steepest descent usa a derivada primeira para determinar a direção para o mínimo. Esta técnica é robusta e é usada para minimizar inicialmente uma estrutura que está longe de um ponto de mínimo. Uma minimização mais eficiente pode ser obtida utilizando-se o algoritmo dos gradientes conjugados ou o método de Newton-Raphson. A técnica de gradientes conjugados utiliza informação das prévias derivadas primeiras para determinar uma direção ótima da busca. O método de Newton-Raphson usa tanto as derivadas primeiras quanto as segundas da função. Assim, além de utilizar as informações do gradiente, utiliza a curvatura para predizer onde a função mudará a direção ao longo do gradiente.

Após a minimização de energia do sistema, este é gradualmente aquecido para a temperatura de interesse T₀, atribuindo as velocidades iniciais de todas as partículas através de uma distribuição de Maxwell-Boltzmann.

As condições iniciais de uma simulação de DM distam muito do equilíbrio termodinâmico.

Assim sendo, os primeiros passos de uma DM constituem o que se denomina de "período de equilibração" (ou inicialização). Nesse estágio, as propriedades do sistema não se mantêm constantes. O período de equilibração é variável e depende do sistema em estudo. Geralmente, dá-se por finalizado quando se alcança o equilíbrio termodinâmico. A partir desse ponto pode-se então gerar as trajetórias da DM e calcular as diferentes propriedades para o sistema de interesse.

Constrangimentos (constraints) e restrições (restraints) moleculares

Os sistemas moleculares apresentam movimentos que ocorrem em diferentes escalas temporais, considerando tanto as forças intermoleculares como as intramoleculares. Normalmente essas últimas são de ordem de magnitude superior às das primeiras, de modo que os movimentos vibracionais são muito mais rápidos que os translacionais e rotacionais. Dessa forma os passos de integração devem ser aumentados, implicando em um aumento substancial do tempo de cálculo. No entanto, se os movimentos mais rápidos forem congelados (subtraídos do cálculo das interações), constrangendo os graus de liberdade correspondentes aos movimentos de maior frequência, pode-se conseguir uma maior eficiência através da utilização de um passo temporal maior na integração das equações do movimento, reduzindo consideravelmente o custo computacional [19].

Berendsen e van Gunsteren verificaram que quando se aplicam constraints aos comprimentos das ligações em proteínas, o passo temporal pode ser aumentado quatro vezes sem que o comportamento físico do sistema seja alterado [41].

Existem vários métodos usados para realizar constrangimentos moleculares [42]. Os mais utilizados são os algoritmos SHAKE [17,43], RATTLE [44] e LINCS [45].

Condições periódicas de contorno

Em simulações de DM, os sistemas possuem normalmente 10²-10⁶ partículas, sendo insignificante em comparação com o número de átomos contidos em sistemas macroscópicos (frações do número de Avogadro). Assim, nas simulações de DM o número relativo de partículas na superfície do sistema é muito maior do que para um sistema macroscópico, causando efeitos de superfície indesejáveis. A forma mais comum de minimizar os efeitos de superfície e obter as propriedades macroscópicas a partir das trajetórias das simulações de DM é a utilização de condições periódicas de contorno (ou fronteira) [18,19]. Efeitos de superfície são eliminados e as partículas consideradas sofrem o efeito das forças como se estivessem no interior do sistema. Nesta técnica, os átomos do sistema são colocados em uma caixa (normalmente cúbica), e a caixa original é replicada em todas as direções do espaço. Durante as simulações, quando um átomo move-se na caixa original, sua imagem periódica, em uma das caixas imagem, move-se da mesma maneira. Se um átomo sai da caixa original, sua imagem entra pela face oposta da caixa imagem, com a mesma velocidade. Assim, o número total de átomos na caixa central e no sistema como um todo é conservado.

Convenção da imagem mínima e raio de corte

A implementação das condições periódicas de contorno requer a avaliação das forças exercidas sobre cada molécula por todas as demais moléculas do sistema. Se um sistema contém N partículas e o potencial se expressa como a soma das interações entre pares de partículas, então existem termos para essa soma. No entanto, as interações com as partículas das caixas imagem também têm que ser incluídas nessa soma. Com isso, a soma constará de infinitos termos e, na prática, sua avaliação requer a realização de aproximações.

Para o tratamento das interações de curto alcance, onde V ~ r⁶ (interações de van der Waals), emprega-se o procedimento denominado convenção da imagem mínima [17-19,46], que é uma conseqüência do uso das condições periódicas de contorno. Nesse caso, é empregado o truncamento do potencial em um raio de corte esférico, geralmente da ordem de L/2 (L é o comprimento da célula de simulação) para evitar a interação da partícula com sua própria imagem ou com a mesma partícula duas vezes [17].

Tratamento das interações de longo alcance

Como dito anteriormente, as interações de vdW podem ser tratadas apenas como interações de curto alcance, uma vez que o potencial de interação intermolecular decai rapidamente com o aumento da distância entre as partículas [1719]. No entanto, para o tratamento das interações de longo alcance (potencial de Coulomb), em que o potencial decai com a distância que é menor do que r^-d (d é a dimensionalidade do sistema), o truncamento abrupto do potencial produz efeitos que, em muitos casos, comprometem severamente os resultados das simulações [47]. Este problema tem sido contornado com a utilização de métodos baseados em somas de Ewald [46,48], como P³M (Particle-Particle-Particle Mesh Ewald) [49], PME (Particle Mesh Ewald) [50] e SPME (Smoth Particle Mesh Ewald) [51], que produzem trajetórias estáveis, ou por métodos que empregam aproximação do campo de reação [17,46,48].

Dinâmica molecular de sistemas solvatados

Existem dois modelos de solvatação que podem ser utilizados para realizar simulações de DM de sistemas solvatados: o modelo de solvatação explícita e o modelo de solvatação implícita. No primeiro modelo, centenas ou milhares de moléculas de solvente são incorporadas explicitamente ao sistema como um componente adicional [52]. Vários potenciais foram desenvolvidos e parametrizados para descrever as moléculas de água explicitamente. Os mais utilizados são os da família TIP (Transferable Intermolecular Potentials), desenvolvidos por Jorgensen e colaboradores [53], TIP3P, TIP4P e TIP5P, os modelos de carga pontual, SPC e SPC/E, desenvolvidos por Berendsen e colaboradores [54], e o modelo ST2 de Stillinger e Rahman [55].

Em alguns casos, tais como sistemas de proteína muito grandes, torna-se necessário a utilização de modelos de solvatação implícita, que reduzem o custo computacional das simulações.

Em geral, um tratamento implícito de solvatação considera a influência média do solvente através da estimativa direta da energia livre de solvatação, definida como o trabalho reversível necessário para transferir o soluto em uma configuração fixa do vácuo para a solução. O método de solvatação implícita não considera os graus de liberdade do solvente explicitamente, mas o trata como um meio contínuo. Métodos empíricos, tais como o modelo de área superficial acessível ao solvente (SASA) [56,57], muitas vezes constituem formas simples e rápidas de avaliação da energia de solvatação com uma precisão comparável aos modelos teóricos. Na abordagem SASA, a energia livre de solvatação do soluto é expressa como a soma das contribuições atômicas, ponderada pelas suas áreas expostas ao solvente. A contribuição de cada átomo é quantificada por um coeficiente de superfície, o qual reflete a hidrofobicidade ou hidrofilicidade do tipo de átomo em questão [56,57].

O modelo de solvatação contínuo tem se destacado para descrever a solvatação eletrostática. Nesta abordagem o soluto é considerado como uma cavidade embebida em um meio dielétrico. A correspondente energia livre de solvatação eletrostática pode ser calculada rigorosamente através das soluções da Equação de Poisson-Boltzmann (PB) [58,59], ou aproximadamente pelo uso do modelo generalizado de Born (GB) [60,61].

Em estudos de associação de macromoléculas biológicas com ligantes, o modelo de solvatação implícita de Born generalizado (GB) é o mais aplicado devido ao seu menor custo computacional [62-65]

A combinação GB/SA tem sido reconhecida como uma excelente escolha para o tratamento implícito de solvatação em simulações biomoleculares [66]. Diversas otimizações dos métodos de solvatação implícita têm sido realizadas [67-70] e encontram-se implementados em diversos pacotes de programas de modelagem molecular como CHARMM [12], AMBER [71] e XPLOR [72]. Adicionalmente, continuam os esforços de pesquisa para a parametrização destes modelos com diferentes solventes [73-78].

Aplicação da dinâmica molecular no planejamento de fármacos

Simulações de DM têm sido rotineiramente utilizadas para investigar a estrutura e função de biomoléculas e seus respectivos complexos, e interações no processo de planejamento de fármacos assistido por computador [79-83]. Esta técnica incorpora flexibilidade de ambos, ligante e receptor, melhorando suas interações e reforçando a complementaridade entre eles [4]. Além disso, a habilidade de incorporar moléculas de solvente nas simulações de sistemas receptorligante é muito importante para o entendimento do papel da água e seus efeitos sobre a estabilidade de complexos proteína-inibidor [4].

Cavallari e colaboradores [84] realizaram simulações de DM para estudar a estabilidade da enzima protease do HIV-1 (HIV-1 PR) em solução aquosa em diferentes estados de enovelamento, contendo diferentes conteúdos de hélices α e folhas β (13 e 25%, respectivamente). O campo de força CHARMM foi utilizado juntamente com o modelo de solvatação TIP3P em sistema NVT em longas trajetórias de 6 ns de simulações para cada um dos três estados da enzima. Interações de longo alcance foram calculadas usando o procedimento SPME e um banho térmico Nosé-Hover, a 300 K. As mudanças conformacionais da proteína observadas durante as simulações de DM foram descritas através dos desvios quadráticos médios (RMDS, equação 12) e as flutuações da estrutrura média da proteína foram expressas em termos de flutuações dos desvios quadráticos médios atômicos (RMFS, equação 13).

onde r_i (t) e r_i (0) são as coordenadas do i-ésimo átomo no tempo t e 0, respectivamente, e N é o número de átomos no domínio de interesse (geralmente C_α ou átomos da cadeia principal).

onde são as coordenadas do i-ésimo átomo no passo de tempo j, , suas posições médias e é o tempo de simulações, expresso como o número total de passos de tempo (time steps) coletados.

Os resultados mostraram boa estabilidade das estruturas em solução aquosa. As conformações obtidas nas simulações de DM para a protease foram complexadas com uma variedade de ligantes através de técnicas de docking, gerando importantes informações sobre os modos de ligação de pequenas moléculas nos diferentes estados de enovelamento da enzima.

Docking molecular combinado com simulações de DM e cálculos de energia livre de ligação (método LIE, Linear Interaction Energy [85-89]) foi realizado por Carlsson e colaboradores [90] para predizer os modos de ligação e estimar as energias livres de ligação de 43 inibidores nãonucleosídeos da enzima transcriptase reversa do HIV-1. A energia livre de ligação foi estimada através da Equação 14, onde os Δs denotam que as diferenças entre as energias para os estados ligados e livres têm sido consideradas. são obtidos das simulações de DM e representam os valores médios das energias de van der Waals e eletrostática, respectivamente, para o ligante solvatado.

As simulações de DM para os estados ligado e livre para cada inibidor foram realizadas pelos autores com o pacote de programas Q, utilizando modelo de solvatação explícita (TIP3P) e o campo de força OPLS. O algoritmo SHAKE foi utilizado para restringir todos os comprimentos de ligação, juntamente com um cutoff de 10 Å. Interações de longo alcance foram tratadas com aproximação do campo de reação local. Trajetórias de 500 ps de simulação foram obtidas com o sistema NVT. Os resultados indicaram boa correlação com os valores de energias livres de ligação experimentais para 39 dos 43 inibidores estudados, sugerindo que a metodologia utilizada é bastante promissora em estudos de planejamento de fármacos baseado em estrutura.

De maneira similar, Nervall e colaboradores utilizaram, em trabalho recente [91], as simulações de DM como ferramenta para cálculos de energia livre de ligação de complexos obtidos pelo método de docking flexível empregando duas estruturas cristalográficas da enzima transcriptase reversa do HIV-1 como receptor e 34 inibidores derivados de benzilpiridinonas. Trajetórias de 650 ps de simulação foram obtidas no sistema NVT, e os resultados mostraram que as energias livres de ligação preditas pelo método LIE foram hábeis em propor o correto modo de ligação, depois validado, além das energias livres relativas obtidas mostrarem boa correlação com os dados experimentais obtidos previamente.

Simulações de DM também foram realizadas por Hu e colaboradores em estudos sobre a interação do inibidor EBR28 com o domínio central da enzima integrase do HIV-1 (HIV-1 IN) nas suas formas monomérica e dimérica [82]. Quatro simulações de DM foram realizadas independentemente com um modelo de solvatação explícito (TIP3P) e campo de força AMBER para os dois sistemas: monômero - complexo monômero/EBR28 (obtido através de docking molecular), e dímero - complexo dímero/EBR28 (obtido através de docking molecular). Condições periódicas de contorno, raio de corte de 0,8 nm e cutoff de 1,2 nm para todas as interações não-ligadas foram empregadas em cada simulação. Uma longa trajetória de 1900 ps de simulação não restringida foi realizada para cada sistema. Os gráficos da energia potencial e do desvio médio quadrático (RMSD) em função do tempo de simulação de cada sistema evidenciaram uma alta estabilidade de todos os sistemas em água, a 300K. Através das simulações de DM associadas aos estudos de docking, os autores conseguiram propor os modos de ligação do inibidor EBR28 com o domínio central na forma monomérica e dimérica, identificando as regiões 'chave' de interação das proteínas.

Adicionalmente, cálculos de energia livre de ligação entre o monômero do domínio central da HIV-1 IN com uma série de peptídeos mutantes foram também realizados. Snapshots das simulações foram extraídos das trajetórias para o cálculo da energia livre de ligação segundo a equação 15:

onde ΔG_bind é a energia livre de ligação e ΔE_MM representa a soma das energias da mecânica molecular no vácuo, a qual pode ser assim repartida: eletrostática, van der Waals e interna. ΔΔG_sol é a energia livre de solvatação, que inclui a energia livre de solvatação polar calculada com o modelo de Born Generalizado (GB) e a parte não-polar obtida com o modelo de área superficial acessível ao solvente (SASA). Finalmente, TΔS é a entropia do soluto. Os valores de energia livre de ligação, assim calculados, mostraram alta correlação com resultados experimentais previamente obtidos, e um plausível mecanismo de inibição foi proposto, fornecendo valiosas informações para o planejamento de fármacos anti-HIV, principalmente pequenos peptídeos.

Cheng e colaboradores [92] realizaram simulações de DM do monômero da enzima protease do vírus HIV-1 (HIV-1 PR) em água, em diversas temperaturas. Uma vez que o sítio catalítico é formado pelas duas subunidades da enzima, Cheng e colaboradores estudaram a dinâmica do monômero com o intuito de obter insights sobre o mecanismo de dimerização da enzima. As simulações foram realizadas com o pacote de programas GROMACS, utilizando-se o campo de força GROMACS e condições periódicas de contorno. A proteína foi solvatada com o modelo de água SPC. Os comprimentos de ligação foram restringidos com o algoritmo LINCS. Um cutoff de 0,9 nm foi utilizado e as interações eletrostáticas foram calculadas utilizando o algoritmo PME. O tempo de simulação foi de 2,0 fs, a 300K. Adicionalmente, trajetórias da ordem de ns de simulações foram também obtidas, variando-se a temperatura (300, 330 e 350 K). Uma análise dos desvios médios quadráticos revelou que a estrutura secundária do monômero se manteve estável durante todas as simulações, nas diferentes temperaturas, com exceção dos N- e C- terminais, os quais mostraram flutuações, sendo essas mais intensas durante as simulações a 350 K. Durante as simulações a 300 e 330K, foi observada a formação de uma folha β entre os N- e C- terminais no final das simulações, o que refletiu em uma estabilidade nos gráficos de RMSD em função do tempo de simulação. A partir dos resultados obtidos, os autores puderam concluir que pequenas moléculas orgânicas ou pequenos peptídeos podem se ligar aos domínios terminais do monômero, estabilizando a estrutura secundária e, dessa forma, interferir no modo de dimerização da protease. Os resultados forneceram novas perspectivas sobre o mecanismo de dimerização da HIV-PR, sugerindo novos modelos de inibidores de dimerização.

Da Silva e colaboradores [93] utilizaram simulações de DM para investigar a dinâmica e a estabilidade de um potencial inibidor, análogo nucleosídeo, da enzima transcriptase reversa do vírus HIV-1 (HIV-1 TR). Trajetórias de 1,5 ns de simulações foram obtidas com o módulo discover do pacote de programas Insight II. O campo de força CVFF foi utilizado no sistema NVT. Baixos valores de RMSD (~0,08 Å) para o complexo, observados durante todo o tempo de simulação, indicaram alta estabilidade do novo potencial inibidor no interior do sítio ativo da enzima, sugerindo que a molécula proposta é um promissor inibidor em potencial para essa enzima e para a quimioterapia anti-retroviral.

Lee e colaboradores [94] realizaram simulações de DM da enzima protease do HIV-1 complexada com álcool difenílico C₆₀ (protótipo do álcool derivado do fulereno C₆₀) em três estados de protonação (desprotonado, monoprotonado e diprotonado) de dois resíduos catalíticos, Asp25(A) e Asp25(B), situados na cavidade da enzima. Trajetórias de 1200 ps de simulação de DM foram obtidas com o pacote de programas AMBER v.7.0 e campo de força AMBER (parm99). Os complexos foram solvatados em camadas de até 10 Å a partir da superfície molecular. O sistema NPT foi tratado sob condições periódicas de contorno a 1 atm de pressão e temperatura de 298 K. Um cutoff de 12 Å foi aplicado para interações não-ligadas, com o método PME. O algoritmo de Berendsen foi empregado a 298 K. Outras análises foram realizadas pelos autores, com as quais foi possível investigar a estabilidade e a análise conformacional dos complexos, utilizando-se para tanto os gráficos de desvios quadráticos médios (RMSD) vs tempo de simulação, os movimentos internos de ângulos torcionais do inibidor (através do monitoramento dos valores de ângulos torcionais específicos), a localização do inibidor dentro da cavidade, a rotação do inibidor ao redor de seu próprio eixo (através do monitoramento das distâncias de átomos específicos no sítio ativo), as ligações de hidrogênio que foram formadas entre o inibidor e os resíduos de ácido aspártico do sítio ativo da enzima (considerando como critério: distância entre doador e aceptor de hidrogênio menor que 4 Å e ângulo entre doador/aceptor/átomo H menor que 60º). Adicionalmente, sabendo-se que o mecanismo catalítico da protease é possivelmente uma combinação entre ligações favoráveis do inibidor com a enzima e o papel das moléculas de água no sítio catalítico, os autores realizaram o monitoramento das posições e orientações das moléculas de água, utilizando funções de distribuição radiais átomo-átomo, , centradas no i-ésimo átomo. Os autores concluíram que o sítio ativo da protease de HIV-1, quando complexado com o álcool difenílico C₆₀ , encontra-se no seu estado diprotonado.

Os resultados forneceram informações importantes sobre as mudanças conformacionais da enzima e do receptor, bem como da estabilidade dos complexos (RMSD variando de 1,4 a 2,2 Å), permitindo identificar o melhor estado de protonação do complexo (estado diprotonado). Além disso, os dados estruturais sobre a ligação e orientação do protótipo de inibidor C₆₀ na cavidade da enzima e os dados de solvatação obtidos das simulações de DM (estabilidade da ligação água-enzima: di- < mono- < desprotonado), forneceram valiosas informações para a descoberta de novos e mais eficientes candidatos a fármacos em terapia anti-HIV.

Conclusões

Existem muitos exemplos onde as simulações de DM têm mostrado papel importante no planejamento de novos fármacos. Em estudos sobre dinâmica e estrutura de sistemas biológicos, essa técnica aponta o efeito que acarreta pequenas mudanças conformacionais sobre o ligante ou sítio de ligação, facilitando a otimização de candidatos a fármacos ou racionalizando o efeito de uma dada mutação. Simulações incluindo enzimas-alvo do vírus HIV-1, passíveis de intervenção farmacológica, bem como seus inibidores, mostram-se bastante úteis no planejamento de candidatos a fármacos anti-retrovirais. A combinação de DM com outras técnicas, como docking molecular e cálculos de energia livre de ligação, são essenciais para o entendimento das interações ligante-receptor, direcionando o planejamento racional de novos e potencias inibidores (ou ligantes de um alvo terapêutico).

Referências

[1] T. Hansson, C. Oostenbrink, W. van Gunsteren, Curr. Opin. Struct. Biol. 12(2002)190.

[2] M. Karplus, J. A. McCammon, Nat. Struct. Biol. 9(2002)646.

[3] J. Norberg, L. Nilsson, Quart. Rev. Biophys. 36 (2003) 257.

[4] H. Alonso, A. A. Bliznyuk, J. E. Gready, Med. Res. Rev. 26 (5) (2006) 531.

[5] S. A. Adcock, J. A. McCammon, Chem. Rev. 106 (5) (2006) 1589.

[6] D. Rognan, Perspect. Drug Discov. Des. 9/10/11 (1998) 181.

[7] U. Burkert, N. L. Allinger, Molecular Mechanics, Molecular Mechanics, ACS Monograph 177, American Chemical Society, Washington, D.C., 1982, 64-72.

[8] H.-D. Höltje, W. Sippl, D. Rognan, G. Folkers, Molecular Modeling: Basic Principles and Applications, Wiley-VCH, Darmstadt, 2nd ed, 2003.

[9] C. L. Brooks III, M. Karplus, B. M. Pettitt, Proteins: A theoretical perspective of dynamics structure and Thermodynamics, Advances in chemical physics, John Wiley & Sons, New York, 1988.

[10] W. F. van Gunsteren, H. J. C. Berendsen, Angew. Chem. Int. Ed. Engl., 29 (1990) 992.

[11] W. F. van Gunsteren, H. J. C. Berendsen, Groningen Molecular Simulation (GROMOS) Library Manual, Biomos, Groningen, The Netherlands, 1987.

[12] B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, M. Karplus, J. Comput. Chem. 4 (1983) 187. A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiorkiewicz-Kuczera, D. Yin, M. Karplus, J. Phys. Chem. B 102 (1998) 3586.

[13] W. F. van Gunsteren, S. R. Billeter, A. A. Eising, P. H. Hünenberger, P. Krüger, A. E. Mark, W. R. P. Scott, I. G. Tironi, Biomolecular Simulation: The GROMOS96 Manual and User Guide, Vdf Hochschulverlang, ETH Zurich Switzerland, 1996.

[14] S. J. Weiner, K. A. Kollman, D. A. Case, U. C. Singh, C. Ghio, G. Alagona, S. Jr. Profeta, P. K. Weiner, J. Am. Chem. Soc. 106 (1984) 765. S. J. Weiner, P. A. Kollman, D. T. Nguyen, D. A. J. Case, J. Comput. Chem. 7 (1986) 230.

[15] W. L. Jorgensen, J. Tirado-Rives, J. Am. Chem. Soc. 110 (1988), 1657. W. L. Jorgensen, D. S. Maxwell, J. Tirado-Rives, J. Am. Chem. Soc. 118 (1996) 11225.

[16] A.T. Hagler, E. Huler, S. Lifson, J. Am. Chem. Soc. 96 (1974) 5319. A. T. Hagler, P. S. Stern, S. Lifson, S. Ariel, J. Am. Chem. Soc. 101 (1979) 813.

[17] M. P. Allen, D. J. Tildesley, Computer simulation of liquids, Oxford University Press, New York, 1987.

[18] W. F. van Gunsteren, H. J. C. Berendsen, Angew. Chem. 102 (1990) 1020.

[19] D. Frenkel, B. Smit, Understanding Molecular Simulations: from Algorithms to Applications, Academic Press, London, 2nd edn., 2002.

[20] S. A. Adcock, J. A. McCammon, Chem. Rev. 106 (2006) 1589.

[21] H. A. Scheraga, M. Khalili, A. Liwo, Annu. Rev. Phys. Chem. 58 (2007) 57.

[22] T. Hansson, C. Oostenbrink, W. van Gunsteren, Curr. Opin. Struct. Biol. 12 (2002) 190.

[23] M. Karplus, J. A. McCammon, Nat. Struct. Biol. 9 (2002) 646.

[24] J. Norberg, L. Nilsson, Quart. Rev. Biophys., 36 (2003) 257.

[25] N. Attig, K. Binder, H. Grubmüller, K. Kremer, Advanced simulations for hydrodynamic problems: Lattice Boltzmann and Dissipative Particle Dynamics, in Computational soft matter: From synthetic polymers to proteins, NIC Series, Jülich, Vol. 23, 2004.

[26] K. Huang, Statistical Mechanics, Wiley John & Sons, New York, 2nd edn,1990.

[27] L. Verlet, Phys. Rev.165 (1967) 201.

[28] D. Beeman, J. Comput. Phys. 20 (1976) 130.

[29] R. W. Hockney, Methods Comput. Phys. 9 (1970) 136.

[30] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, J. R. Haak, J. Chem. Phys. 81 (1984) 3684.

[31] S. Nosé, Mol. Phys. 52 (1984) 255.

[32] W. G. Hoover, Phys. Rev. A 31 (1985) 1695.

[33] M. Parrinello, A. Rahman, J. Appl. Phys. 52 (1981) 7182.

[34] P. G. Pascutti, Introdução à Modelagem e Dinâmica Molecular. In: Pedro G Pascutti. (Org.). Introdução à Modelagem e Dinâmica Molecular v. 1, p. 1-38, 2002.

[35] G. Arfken, Mathematical Methods for Physicists, Academic Press, Orlando, FL, 3rd ed., 1985.

[36] D. Menzel (Ed), Fundamental Formulas of Physics, New York, 2nd ed. 1960.

[37] P. M. Morse, H. Feshbach, Methods of Theoretical Physics, Part I, McGraw-Hill, New York 1953.

[38] M. R. Hestenes, E. Stiefel, J. Res. Nat. Bur. Stand. 49(6) (1952) 409.

[39] R. Fletcher, Practical Methods of Optimization, Wiley, Chichester, 1980.

[40] E. Kreyszig, Advanced Engineering Mathematics, 7th ed., John Wiley, New York, 1993.

[41] W. F. van Gunsteren, H. J. C. Berendsen, Mol. Phys. 34 (1977) 1311.

[42] D. van der Spoel, E. Lindahl, B. Hess, A. R. van Buuren, E. Apol, P. J. Meulenhoff, D. P. Tieleman, A. L. T. M. Sijbers, K. A. Feenstra, R. van Drunen, H. J. C. Berendsen, "Gromacs User Manual version 3.2", Groningen, The Netherlands, 2004.

[43] J. P. Ryckaert, G. Ciccotti, H. J. C. Berendsen, J. Comput. Phys. 23 (1977) 327.

[44] H. C. Andersen, J. Comput. Phys. 52 (1983) 24.

[45] B. Hess, H. Bekker, H. J. C. Berendsen, J. G. E. M. Fraaije, J. Comput. Chem. 18 (1997) 1463.

[46] A. R. Leach, Molecular Modelling Principles and Applications, Addison Wesley Longman, Essex, England, 1996.

[47] T. E. Cheatham, J. L. Miller, T. Fox, T. A. Darden, P. A. Kolman, J. Am. Chem. Soc. 117 (1995) 4193.

[48] I. G. Tironi, R. Sperb, P. E. Smith, W. F. van Gunsteren, J. Chem. Phys. 102 (1995) 6199.

[49] R.W. Hockey, J.W. Eastwood, Computer Simulation using Particles, McGraw-Hill, New York, 1981.

[50] T. Darden, D. York, L. Pedersen, J. Chem. Phys. 98 (1993) 10089.

[51] U. Essmann, L. Perera, M. L. Berkowitz, T. A. Darden, H. L. Lee, L. G. Pedersen, J. Chem. Phys. 103 (1995) 8577.

[52] A. R. Bizzarri, S. Cannistraro, J. Phys. Chem. B 106 (2002) 6617.

[53] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impy, M. L. Klein, J. Chem. Phys. 79 (1983) 926.

[54] H. J. C. Berendsen, J. R. Grigera, T. P. Straatsma, J. Phys. Chem. 91 (1987) 6269.

[55] F. H. Stillinger, A. Rahman, J. Chem. Phys. 60 (1974) 1545.

[56] B. Roux, T. Simonson, Biophys. Chem. 78(1-2) (1999) 1.

[57] P. Ferrara, J. Apostolakis, A. Caflisch, Proteins Struct. Funct. Genet. 46 (2002) 24.

[58] M. K. Gilson, M. E. Davis, B. A. Luty, J. A. Mccammon, J. Phys. Chem. 97 (1993) 3591.

[59] N. A. Baker, Curr. Opin. Struct. Biol. 15 (2005) 137.

[60] W. C. Still, A. Tempczyk, R. C. Hawley, T. Hendrickson, J. Am. Chem. Soc. 112 (1990) 6127.

[61] D. Bashford, D. A. Case, Annu. Rev. Phys. Chem. 51 (2000) 129.

[62] X. Zou, Y. Sun, I. D. Kuntz, J. Am. Chem. Soc. 121 (1999) 8033.

[63] N. Majeux, M. Scarsi, J. Apostolakis, C. Ehrhardt, A. Caflisch, Proteins 37 (1999) 88.

[64] R. D. Taylor, J. W. Essex, P. J. Jewsbury, J. Comput. Chem. 24 (2003) 1637.

[65] H. Y. Liu, X. Zou, I. D. Kuntz, J. Phys. Chem. B 108 (2004) 5453.

[66] M. Feig, C. L. Brooks, III, Curr. Opin. Struct. Biol. 14 (2004) 217.

[67] J. Chen, C. L. Brooks, J. Khandogin, Curr. Opin. Struct. Biol. 18(2) (2008) 140.

[68] M. S. am Busch, A. Lopes, A. Najette, C. Bathelt, T. Simonson, BMC Bioinformatics 9 (2008) No pp. given.

[69] J. Chen, W. Im, C. L. Brooks, III, J Am Chem Soc 128 (2006) 3728.

[70] S. Jang, E. Kim, Y. Pak, Proteins 66 (2007) 53.

[71] D. A. Case, T. E. Cheatham, III, T. Darden, H. Gohlke, R. Luo, K. M. Merz, Jr., A. Onufriev, C. Simmerling, B. Wang, R. Woods. J. Computat. Chem. 26 (2005) 1668.

[72] C. D. Schwieters, G. M. Clore, J Magn Reson 149 (2001) 239.

[73] J. Li, T. Zhu, G. D. Hawkins, P. Winget, D. A. Liotard, C. J. Cramer, D. G. Truhlar, Theor. Chem. Acc. 103 (1999) 9.

[74] J. Li, G. D. Hawkins, C. J. Cramer, D. G. Truhlar, Chem. Phys. Lett. 288 (1998) 293

[75] M. Orozco, F. J. Luque, Chem. Phys. 182 (1994) 237.

[76] G. I. Almerindo, D. W. Tondo, J. R. Pliego Jr., J. Phys. Chem. A 108 (2004) 166.

[77] F. J. Luque, M. Bachs, C. Alemán, M. Orozco, J. Comp. Chem. 17 (1996) 806.

[78] F. J. Luque, Y. Zhang, C. Alemán, M. Bachs, J. Gao, M. Orozco, J. Phys. Chem. 100 (1996) 4269.

[79] Y. Joti, H. Nakagawa, M. Kataoka, A. Kitao, J. Phys. Chem. B 112(11) (2008) 3522.

[80] N. Diaz-Vergara, A. Pineiro, J. Phys. Chem. B 112(11) (2008) 3529.

[81] N. J. English, J. Mol. Model. 13(10) (2007) 1081.

[82] J. P. Hu, X.O. Gong, J. G. Su, W. Z. Chen, C. X. Wang, Biophys. Chem. 132(2-3) (2008) 69.

[83] A. Lavecchia, S. Cosconati, E. Novellino, E. Calleri, C. Temporini, G. Massolini, G. Carbonara, G. Fracchiolla, F. Loiodice, J. Mol. Graph. Model. 25(6) (2007) 773.

[84] M. Cavallari, C. Ghio, S. Monti, M. Ferrario, A. Maritan, P. Carloni, THEOCHEM 769(1-3) (2006) 111.

[85] J. Åqvist, C. Medina, J. Samuelsson, Protein Eng. 7 (1994) 385.

[86] T. Hansson, J. Åqvist, Protein Eng. 8 (1995) 1137.

[87] J. Åqvist, J. Comput. Chem. 17 (1996) 1587.

[88] W. Wang, J. Wang, P. A. Kollman, Proteins Struct. Funct. Genet. 34 (1999) 395.

[89] Y. Y. Sham, Z. T. Chu, A. Warshel, Proteins Stuct. Funct. Genet. 39 (2000) 393.

[90] J. Carlsson, L. Boukharta, J. Aqvist, J. Med. Chem. 51(9) (2008) 2648.

[91] M. Nervall, P. Hanspers, J. Carlsson, L. Boukharta, J. Aqvist, J. Med. Chem. 51(9) (2008) 2657.

[92] M. C. Yan, Y. Sha, J. Wang, X. Q. Xiong, J. H. Ren, M. S. Cheng, Proteins 70(3) (2008) 731.

[93] C. H. T. P. da Silva, I. Carvalho, C. A. Taft, J. Theor. Comp. Chem. 5(3) (2006) 579.

[94] V. S. Lee, P. Nimmanpipug, O. Aruksakunwong, S. Promsri, P. Sompornpisut, S. Hannongbua, J. Mol. Graph. Model. 26(2) (2007) 558.

Recebido em 18 de junho de 2008

Aceito em 14 de outubro de 2008

Datas de Publicação

• Publicação nesta coleção
  12 Jan 2009
• Data do Fascículo
  Dez 2008

Histórico

• Recebido
  18 Jun 2008
• Aceito
  14 Out 2008