Second Brazilian School on Bioinformatics 2009
(EBB 2009)
Curso
Modelagem e Dinâmica Molecular
Prof. Hermes Luís Neubauer de Amorim - hlnamorim@yahoo.com.br
(Labortatório de Bioinformática Estrutural - LaBiE - ULBRA)
Prof. Paulo Augusto Netz - netz@iq.ufrgs.br
(Grupo de Química Teórica - GQT - UFRGS)
Porto Alegre, julho de 2009
Tutorial
Simulação da dinâmica molecular do complexo ternário entre uma fosfolipase A2 e a aspirina
A simulação dinâmica molecular está se tornando, nos últimos anos, uma das ferramentas mais importantes para o estudo da estrutura, da dinâmica e de outras propriedades de macromoléculas biológicas, desde proteínas até ácidos nucléicos. Estes estudos têm uma importância central na elucidação da função biológica destas moléculas, no desenho racional de fármacos e nos estudos da interação destas moléculas com diferentes meios, seja em condições fisiológicas, alteradas ou patológicos.
Um dos pacotes de programas mais usados neste tipo de estudo é o GROMACS. Trata-se de um programa de alto desempenho, especialmente adequado para a simulação de macromoléculas de interesse biológico. Ele trabalha com um campo de forças simples, harmônico, com a convenção united atom e sem termos cruzados. Um dos pontos fortes deste programa, além da sua velocidade, é a facilidade de uso.
Nesta prática, iremos usar o GROMACS 4.0.2 para simular o complexo ternário entre uma fosfolipase A2 de veneno de serpente e a aspirina. Fosfolipases A2 (FLA2) catalisam a hidrólise da posição sn-2 de glicerofosfolipídios de membrana, produzindo 1-acilfosfolipidios e ácidos graxos livres. Esta reação é particularmente importante quando o ácido graxo gerado é o ácido araquidônico, já que este último é convertido por enzimas metabólicas em vários compostos bioativos (eicosanoides), tais como as prostagrandinas e os leucotrienos. Outros produtos da reação, por exemplo ácido lisofosfatídico e lisofosfatidilcolina, também são biologicamente ativos, sendo precursores de outros potentes mediadores bioativos, tais como o fator de ativação plaquetária (PAF). Há evidências de que a ação anti-inflamatória da aspirina decorre, em parte, da inibição de enzimas da família FLA2.
![]()
Primeiro, faremos o download da estrutura do complexo enzima-fármaco a partir do banco de dados Protein Data Bank (PDB). A estrutura do complexo que iremos estudar foi determinada através de cristalografia de raio-X, a uma resolução de 1,9 Å. Para maiores detalhes, a referência primária para esta estrutura é: Singh, RK et al. “Aspirin induces its anti-inflammatory effects through its specific binding to phospholipase A2: crystal structure of the complex formed between phospholipase A2 and aspirin at 1.9 angstroms resolution” J. Drug Target 2005, 13(2):113-9 (PubMed).
Após, iremos gerar os arquivos de topologia e de entrada para o programa. A seguir, a estrutura será hidratada e procede-se a uma minimização da energia. As etapas seguintes são uma simulação com restrição de posições e uma simulação dinâmica molecular propriamente dita, para explorar as propriedades termodinâmicas e estruturais do sistema.
1.1 Recuperando a estrutura cristalográfica do complexo enzima-fármaco
Antes de mais nada, abra um terminal (com o cursor sobre a área do Desktop, clique com o botão direito do mouse - escolha a opção Terminal ou Konsole) e crie uma pasta com seu nome no diretório /home/visit/mini6/
mkdir seunome
vá para a pasta criada
cd seunome
Obtenha a estrutura da proteína complexo entre a fosfolipase A2 e a aspirina, código pdb (PDB ID) 1OXR do endereço do PDB. Para isto, digite o PDB ID na área correspondente de busca do site do Protein Data Bank. Uma vez encontrada a estrutura, faça download do arquivo PDB, da área de download files, localizada no canto superior direito da página - opção PDB File (Text). Clique com o botão direito do mouse no link para armazenar o arquivo na pasta com o seu nome (/home/visit/mini6/seunome). Mude o nome do arquivo baixado para cplx.pdb.
Use o programa VMD para visualizar a complexo:
vmd cplx.pdb
Visualização com o programa VMD
Na janela VMD Main clique em Graphics e depois em Representations
No campo Selected Atoms digite protein e depois tecle Enter
Em Drawing Method escolha NewCatoon
Em Coloring Method escolha Secondary Structure
Clique em Create Rep
No campo Selected Atoms digite resname AIN e depois tecle Enter
Em Drawing Method escolha Licorice
Em Coloring Method escolha Name
1.2 Preparando o ligante
O programa GROMACS gera a topologia e as coordenadas para a simulação de DM somente de moléculas que estão parametrizadas em seu banco de dados. Assim, a topologia e as coordenadas da aspirina devem ser obtidas a partir de um outro programa, o PRODRG.
Primeiro, precisamos extrair as coordenadas do ligante (aspirina) no arquivo cplx.pdb.
Abra o arquivo cplx.pdb com um editor de texto qualquer. Por exemplo:
gedit cplx.pdb
Role até a parte final do arquivo. Selecione e copie as coordenadas da estrutura da aspirina (linhas que começam com HETATM..........AIN), destacadas abaixo:
Utilizando um navegador de internet, vá até a página do servidor PRODRG, versão beta, e cole as coordenadas do fármaco na janela de input.
Antes de prosseguir, selecione No na opção EM (minimização de energia).
Depois, clique em Run PRODRG.
Na janela que abre, localize o campo “Coordinates” e clique em “polar/aromatic H's”
Selecione e copie as coordenadas da aspirina no formato GROMACS.
Sem fechar o navegador, abra um documento em branco no editor de texto e cole as coordenadas copiadas previamente. Salve o arquivo como ain.gro
Volte à página do PRODRG. Role para cima e, no campo “Docking / MD simulations”, clique na opção “GROMACS [D] (topology)”.
Selecione e copie a topologia da aspirina no formato GROMACS.
No editor de texto, abra um documento em branco e cole a topologia.
Salve o arquivo como ain.itp
Um problema nos arquivos de topologia obtidos a partir do PRODRG são as cargas atômicas geradas pelo programa. Compare abaixo as cargas geradas para a aspirina pelo PRODRG com as cargas geradas com o método ab initio CHELPG 6-31G**(d,p).
Estrutura | Átomo | Cargas PRODRG | Cargas CHELPG |
| O1 C7 O2 C3 C4 H4 C5 H5 C6 H6 C1 H1 C2 O3 C8 O4 C9 | -0.255 0.396 -0.255 -0.001 0.012 0.046 0.012 0.045 0.007 0.026 0.007 0.026 0.092 -0.097 0.231 -0.407 0.115 | -0.857 0.966 -0.857 -0.178 -0.062 0.117 -0.175 0.086 -0.068 0.082 -0.304 0.126 0.431 -0.558 0.995 -0.667 -0.074 |
Neste caso, é necessário substituir as cargas atômicas no arquivo ain.itp por cargas melhores.
Em função do tempo reduzido, já preparamos um arquivo de topologia do ligante com cargas CHELPG 6-31G**(d,p). Ele está identificado na pasta "tutorial" (/home/visit/mini6/tutorial/) com o nome ain_c.itp.
Assumindo que você está na pasta criada com seu nome, digite no terminal:
cd ../tutorial
Copie o arquivo ain_c.itp para sua pasta
cp ain_c.itp ../seunome
Volte para sua pasta
cd ../seunome
1.3 Gerando a topologia e as coordenadas da proteína no formato GROMACS a partir do arquivo .pdb
Já que o programa GROMACS não gera a topologia e as coordenadas para moléculas que não sejam aminoácidos ou nucleotídeos, é necessário remover as coordenadas do ligante no arquivo que contém a estrutura da proteína.
Para tanto, abra o arquivo cplx.pdb no editor de texto.
Delete as coordenadas da aspirina (linhas que começam com AIN) e as coordenadas que seguem (águas cristalográficas), até o final do arquivo. Salve o arquivo com as coordenas removidas com o nome prot.pdb.
No terminal, use o programa pdb2gmx, do pacote do GROMACS, para gerar o arquivo de topologia (extensão top) e o arquivo de entrada do programa, com as coordenadas cartesianas da proteína (extensão gro).
g_pdb2gmx –f prot.pdb –o prot.gro –p prot.top
O programa pede a escolha do campo de forças. Escolha GROMOS96 53a6 (opção 4 na versão 4.0 do GROMACS).
Examine o arquivo prot.top abrindo-o com o editor de texto.
Agora, é necessário acrescentar a topologia do ligante com as cargas corrigidas (ain_c.itp) no arquivo de topologia geral (prot.top). Também é necessário colar as coordenadas da aspirina no formato GROMACS (ain.gro) no arquivo que contém as coordenadas da fosfolipase A2 (prot.gro).
Abra o arquivo prot.top com o editor de texto. Vá até a parte final do arquivo e acrescente as linhas que incluem a topologia do ligante, conforme destacado abaixo.
Salve o arquivo como cplx.top.
Abra o arquivo ain.gro e selecione e copie as coordenadas da aspirina.
Abra o arquivo prot.gro com o editor de texto. Altere o número que consta na segunda linha do arquivo para 1206 (número de coordenadas).
Depois, role até a parte final do arquivo e cole as linhas que incluem as coordenadas do ligante, conforme destacado abaixo.
Salve o arquivo como cplx.gro.
1.4 Inserindo as condições periódicas de contorno e solvente
A seguir, definimos o sistema como centrado numa caixa cúbica, cujas faces distam no mínimo 1,0 nm da proteína. Esta caixa será então “preenchida” com moléculas de água (do tipo SPC). O arquivo de topologia é alterado e é gerado um novo arquivo de coordenadas, desta vez contendo proteína e água. Use o programa vmd para visualizar a estrutura final obtida.
g_editconf –bt cubic –f cplx.gro –o cplx.gro –c –d 1.0
g_genbox –cp cplx.gro –cs spc216.gro –o cplx_box.gro –p cplx.top
tail cplx_box.gro
Também é possível solvatar o complexo com outro solvente, como metanol, DMSO, tetracloreto de carbono, etc.
1.5 Minimização de energia
A próxima etapa é a minimização de energia. A maior parte das ações do GROMACS utiliza a combinação dos programas grompp e mdrun. O primeiro utiliza parâmetros contidos num arquivo editável de extensão mdp para produzir um arquivo binário de extensão tpr, o qual irá servir de entrada ao programa mdrun. Na minimização, usamos um arquivo em.mdp, disponibilizado na pasta tutorial.
Para obtê-lo, digite no terminal:
cd ../tutorial
Copie o arquivo em.mdp para sua pasta
cp em.mdp ../seunome
Volte para sua pasta
cd ../seunome
Analise o arquivo em.mdp:
cpp = /lib/cpp
define = -DFLEX_SPC
constraints = none
integrator = steep
nsteps = 500
;
; Energy minimizing stuff
;
emtol = 2000
emstep = 0.01
nstcomm = 1
ns_type = grid
rlist = 1
rcoulomb = 1.0
rvdw = 1.0
Tcoupl = no
Pcoupl = no
gen_vel = no
Rode em seqüência os programas grompp e mdrun, como a seguir:
g_grompp –v –f em.mdp –c cplx_box.gro –o cplx_em.tpr –p cplx.top
(este comando gera o arquivo binário "em.tpr", que contém as informações - arquivos de input - sobre as coordenadas iniciais - arquivo "cplx_box.gro" -, a topologia - aquivo "cplx.top" - e os parâmetros do cálculo - aquivo "em.mdp")
g_mdrun –v –s cplx_em.tpr –o cplx_em.trr –c cplx_em.gro –g emlog
(este comando roda a minimização de energia; o flag –v ativa o método verbose; o flag –s aponta para o arquivo criado pelo grompp, com extensão .tpr; o flag –c indica no nome do arquivo final com as coordenadas minimizadas; o flag –o indica no nome do arquivo de saída com a trajetória da minimização de energia, o qual não é importante para cálculos deste tipo)
O cálculo termina quando a minimização converge ou quando o número de passos previamente estabelecido atinge o limite.
1.6 Dinâmica molecular com restrição de posições
A próxima etapa é a dinâmica molecular com restrição de posições (position restraints) na qual a macromolécula é restrita, mas as moléculas de solvente têm liberdade de movimento. Uma dinâmica de 1,0 ps com estas condições leva a uma relaxação das interações desfavoráveis soluto-solvente. Crie, com os parâmetros abaixo relacionados, um arquivo de nome pr.mdp.
Analise o arquivo pr.mdp disponibilizado para o tutorial:
title = Yo
cpp = /lib/cpp
define = -DPOSRES
constraints = all-bonds
integrator = md
dt = 0.002 ; ps !
nsteps = 500 ; total 1 ps.
nstcomm = 1
nstxout = 250
nstvout = 1000
nstfout = 0
nstlog = 100
nstenergy = 100
nstlist = 10
ns_type = grid
rlist = 1.0
coulombtype = PME
rcoulomb = 1.0
vdwtype = cut-off
rvdw = 1.4
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
; Berendsen temperature coupling is on in three groups
Tcoupl = berendsen
tau_t = 0.1 0.1 0.1
tc_grps = Protein SOL AIN
ref_t = 100 100 100
; Pressure coupling is not on
Pcoupl = parrinello-rahman
pcoupltype = isotropic
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Generate velocites is on at 100 K.
gen_vel = yes
gen_temp = 100.0
gen_seed = 173529
A seguir, rode os programas grompp e mdrun:
g_grompp –f pr.mdp –o cplx_pr.tpr –c cplx_em.gro –r cplx_em.gro –p cplx.top
(este comando gera o arquivo binário "pr.tpr", que contém as informações sobre as coordenadas iniciais, obtidas no final da minização de energia - arquivo "cplx_em.gro" -, a topologia - aquivo "cplx.top" - e os parâmetros do cálculo - aquivo "pr.mdp")
g_mdrun –v –s cplx_pr.tpr –e cplx_pr.edr –o cplx_pr.trr –c cplx_pr.gro –g prlog >& pr.job &
(este comando roda a DM com restrição de posições; o flag –c indica no nome do arquivo final da simulação; o flag –o indica no nome do arquivo de saída com a trajetória da simulação; o flag –e indica no nome do arquivo de saída com a trajetória das energias)
tail –f pr.job
para matar o comando tail, tecle simultaneamente Ctrl C
Analise a energia da simulação
g_energy -f cplx_pr.edr -o out -w
9 0
O resultado desta simulação pode ser visualizado como programa VMD
vmd cplx_pr.gro
Visualização com o programa VMD
Na janela VMD Main clique em Graphics e depois em Representations
No campo Selected Atoms digite protein e depois tecle Enter
Em Drawing Method escolha NewCatoon
Em Coloring Method escolha Secondary Structure
Clique em Create Rep
No campo Selected Atoms digite resname AIN e depois tecle Enter
Em Drawing Method escolha Licorice
Em Coloring Method escolha Name
1.7 Dinâmica molecular sem restrições
A última etapa da simulação consiste na dinâmica molecular propriamente dita, na qual todos os graus de liberdade são amostrados. O arquivo de parâmetros correspondente tem o nome 100ps.mdp.
Para obtê-lo, digite no terminal
cd ../tutorial
Copie o arquivo 100ps.mdp para sua pasta
cp 100ps.mdp ../seunome
Volte para sua pasta
cd ../seunome
Analise o arquivo 100ps.mdp disponibilizado para o tutorial:
title = fws
cpp = /usr/bin/cpp
constraints = all-bonds
integrator = md
tinit = 0
dt = 0.002 ; ps !
nsteps = 50000 ; total 100 ps.
nstcomm = 1
nstxout = 500
nstvout = 0
nstfout = 0
nstlog = 100
nstenergy = 100
nstlist = 10
ns_type = grid
rlist = 1.0
coulombtype = PME
rcoulomb = 1.0
vdwtype = cut-off
rvdw = 1.4
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
; Berendsen temperature coupling is on in three groups
Tcoupl = berendsen
tau_t = 0.1 0.1 0.1
tc_grps = Protein SOL AIN
ref_t = 300 300 300
; Pressure coupling is on
Pcoupl = parrinello-rahman
pcoupltype = isotropic
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Generate velocites is on at 300 K.
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529
O uso seqüencial dos programas grompp e mdrun completa a simulação:
g_grompp –f 100ps.mdp –o cplx_100ps.tpr –c cplx _pr.gro –p cplx.top
g_mdrun –v –s cplx_100ps.tpr –e cplx_100ps.edr –o cplx_100ps.trr –c cplx_100ps.gro –g 100pslog >& 100ps.job &
tail –f 100ps.job
para encerrar o programa tail, tecle simultaneamente Ctrl C
Atenção! feche o terminal digitando exit (não utilize o mouse)
1.8 Análise das simulações
Há diferentes âmbitos de resultados da simulação:
Visualização da dinâmica:
vmd cplx_100ps.gro
Visualização da estrutura com o programa VMD
Proceda como anteriormente
Visualização de trajetória com o programa VMD
Na janela VMD Main clique sobre cplx_100ps.gro de forma a linha fique amarela
Em File escolha Load File Into Molecule...
No campo File Name use Browse para selecionar o arquivo cplx_100ps.trr
OK
Load
Na janela Graphical Representations, clique no botão Play localizado no canto inferior direito
Análise da energia.
g_energy –f cplx_100ps –o cplx100ps_ene (com a escolha dos termos adequados)
Obs: os arquivos .xvg podem ser visualizados com o programa xmgrace:
xmgrace cplx100ps_ene.xvg (e similar para os outros arquivos)
Análise da evolução da estrutura média através do desvio médio quadrático (RMSD) das posições atômicas:
g_rms –s cplx_pr –f cplx_100ps –o cplx100ps_rmsd (escolha o grupo carbono-alfa, 3)
O RMSD convergiu no tempo simulado? Caso contrário, o que isto significa?
Análise do raio de giro:
g_gyrate –s cplx_pr –f cplx_100ps –o cplx100ps_rg (ecolha o grupo proteína, 1)
O raio de giro mudou durante a simulação?
Análise da mobilidade relativa de segmentos:
g_rmsf –s cplx_pr –f cplx_100ps –res –o cplx100ps_rmsf (ecolha o grupo proteína, 1)
Quais segmentos são mais flexíveis durante a simulação?
Análise do gráfico de Ramachandran:
g_rama –s cplx_pr –f cplx_100ps –o cplx100ps_rama
Análise da superfície acessível ao solvente (SAS):
g_sas –s cplx_pr –f cplx_100ps –o cplx100ps_sas (ecolha o grupo proteína, 1)
Análise das ligações de hidrogênio:
g_hbond –s cplx_100ps –f cplx_100ps –num cplx100ps_hbintra (ecolha o grupo proteína, 1, duas vezes)
g_hbond –s cplx_100ps –f cplx_100ps –num cplx100ps_hbprotein-drug (proteína-ligante; ecolha o grupo proteína, 1, e depois ligante, 12)
Postagens
