TUTORIAL
Realizando uma simulação de dinâmica molecular no GROMACS (versão 4.6.7)
1) Introdução
Este tutorial busca apresentar as etapas necessárias
para a realização de uma simulação de dinâmica molecular, utilizando como
referência o pacote de programas GROMACS 4.6.7. Esta versão do
software pode ser baixada e instalada em seu computador pessoal seguindo-se as
instruções contidas no site oficial.
Alternativamente, você também pode utilizá-lo por meio do cluster do Departamento Acadêmico de Física (DAFIS), da Universidade
Tecnológica Federal do Paraná (UTFPR). Basta ser um usuário do servidor, acessando-o
remotamente. Isso deverá ser feito de duas formas num sistema operacional Linux:
1) Secure
Shell (SSH): conexão
com o servidor via terminal.
Através deste, você poderá executar comandos GROMACS a partir da submissão de um job (script PBS no qual são definidos principalmente os comandos e os recursos do cluster a serem usados). Para se conectar, abra um terminal em sua área de
trabalho e digite:
$ ssh –p xxxx usuario@xxxx.edu.br
Após teclar “Enter”,
faça o login fornecendo sua senha. Você estará no ambiente “server”. Para submeter seu job,
dirija-se pela linha de comando à pasta que contém seu arquivo job.sh e
digite:
$ qsub job.sh
Agora seu trabalho entra no sistema de filas do servidor, passando a ser
executado assim que houver nós disponíveis.
OBSERVAÇÃO: GROMACS é um programa que pode exigir bastante interação do usuário, especialmente se este for um iniciante. Como será mostrado neste tutorial, uma série de parâmetros, necessários para a construção dos arquivos básicos do programa, pode ser definida por você. Apesar do uso de flags complementares tornar tais interações desnecessárias, recomenda-se (por questões de aprendizagem) que a execução da maior parte dos comandos ocorra da forma como exposto aqui num computador pessoal. Execute no servidor apenas tarefas que exijam maior tempo e recurso computacional (comandos com o programa mdrun nas duas últimas etapas do tutorial). Aqui, encontra-se disponível um modelo de script PBS para isso. Obtenha o seu arquivo job.sh procedendo da seguinte forma: copie e cole o conteúdo do link anterior num editor de texto, modifique alguns de seus termos para adaptá-los aos nomes dos seus documentos e pastas e salve o arquivo final em formato .sh.
2) SSH File Transfer Protocol (SFTP): conexão com o servidor via um programa gráfico
(Fillezila, por exemplo). A partir da interface deste, você poderá realizar
transferências de arquivos do seu computador para o servidor e vice-versa. Se
você usa o Filezilla, abra-o digitando o seguinte comando em um novo terminal:
$
filezilla
Quando a interface gráfica surgir em sua tela,
clique seguidamente em “Gerenciador de Sites”
(botão no canto superior esquerdo), “Meus
Sites” e “Novo Site”. Preencha
a aba “Geral” que aparece em sua
tela fornecendo os dados necessários para sua conexão exatamente como mostra a
Figura 1.
Figura 1.
Modelo de preenchimento de dados na aba “Geral’.
Após a realização desses dois procedimentos de
conexão, você estará devidamente conectado ao servidor.
Organize-se
Como você verá, sempre a cada etapa, arquivos são requeridos e gerados. Sendo assim, para que a realização de sua simulação ocorra de forma organizada (seja em seu computador pessoal ou no cluster), crie pastas em seu diretório de trabalho “tutorial” (nomeando-as sem acento) para executar cada uma das etapas que virão a seguir: preparação de arquivos, minimização de energia, restrição, dinâmica e análise (mkdir preparacao minimizacao restricao dinamica analise).
Como você verá, sempre a cada etapa, arquivos são requeridos e gerados. Sendo assim, para que a realização de sua simulação ocorra de forma organizada (seja em seu computador pessoal ou no cluster), crie pastas em seu diretório de trabalho “tutorial” (nomeando-as sem acento) para executar cada uma das etapas que virão a seguir: preparação de arquivos, minimização de energia, restrição, dinâmica e análise (mkdir preparacao minimizacao restricao dinamica analise).
Objeto
de trabalho
Para a realização deste tutorial, será
utilizada a estrutura resolvida pelo método 2D-NMR de TS-Kappa, proteína tóxica
produzida pelo escorpião amarelo brasileiro Tityus
serrulatus, e depositada no banco de dados Protein Data Bank (PDB) sob o
código 1TSK. Tal proteína (Figura 2)
é de tamanho relativamente pequeno (35 resíduos), além de contar com a presença
de três pontes dissulfeto, o que contribui para a redução no espaço das configurações espaciais de sua cadeia polipeptídica. Dessa forma, o tempo das simulações e o esforço computacional exigidos são bastante reduzidos, o ideal para uma simulação simples como essa. Após baixar o arquivo 1tsk.pdb, coloque-o na pasta referente à etapa de preparação de arquivos e inicie o tutorial.
Figura 2. Estrutura 3D da proteína TS-Kappa
depositada no PDB.
2) Preparação
de arquivos
Conversão
de formato
Comece os
trabalhos realizando uma conversão de extensão: converta o arquivo .pdb para um
de extensão .gro, o formato utilizável no GROMACS. Para isso, dirija-se à pasta "preparacao" através da linha de comando e digite:
pdb2gmx –f 1tsk.pdb –o
proteina01.gro –p topologia.top –ignh –inter
Onde -f indica arquivo
de entrada, -o indica arquivo de saída e –p é indicativo de que um arquivo de
topologia está entrando ou sendo gerado a partir do processamento (o caso desta
linha de comando).
Ao todo, o comando significa
que o programa pdb2gmx processa o arquivo de entrada no formato .pdb
(1tsk.pdb), gerando os correspondentes arquivo estrutural no formato .gro
(proteina01.gro) e de topologia (topologia.top). Para a construção deste,
pdb2gmx, assim como mais programas a seguir, toma como base as informações contidas
em uma série de bancos de dados.
Importante: As principais diferenças entre os
arquivos estruturais, além do formato, é que o arquivo .gro também pode incluir
dados de velocidade, porém não possui suporte para identificação de cadeias
como um arquivo .pdb possui. Ou seja, se você não precisa de dados sobre
velocidade e/ou trabalha com estruturas com mais de uma cadeia, o melhor seria
gerar um arquivo em formato .pdb a partir da opção –o.
As demais opções são
complementares e opcionais (sem elas, o programa realiza a configuração padrão):
–ignh indica que pdb2gmx deve desconsiderar os átomos de hidrogênio do arquivo
estrutural de partida durante o processamento e –inter permite que o usuário
faça escolhas, por meio do prompt, quanto a alguns parâmetros para a construção
do arquivo de topologia. São eles:
- Campo de força: conjunto de parâmetros aplicáveis numa função potencial com o objetivo de caracterizar todas as interações intra e intermoleculares do sistema. A versão 4.6.7 de GROMACS conta com 19 opções de campos de força (Figura 3), sendo quatro obsoletos (DEPRECATED), ou seja, a utilização deles não é recomendável uma vez que versões atualizadas (e, possivelmente, mais confiáveis) dos mesmos conjuntos de parâmetros encontram-se disponíveis. Para este tutorial, escolha a opção 15 (OPLS-AA/L).
Figura 3. Lista de campos de força
presentes na versão 4.6.7 de GROMACS.
- Tipo de modelo de água: O modelo padrão para representar moléculas de água é o Simple Point Charge (SPC). Por meio dele, átomos de H e de O são tratados como pontos de carga concentrada: carga positiva é atribuída a átomos de H e carga negativa ao átomo de O, de tal forma que a carga total da molécula de água seja nula. Valores para distância de ligação OH e de ângulo HOH variam muito dependendo do tipo de modelo SPC (Figura 4). O modelo recomendado, e que será utilizado neste tutorial, é o TIP4P (modelo de quatro sítios de interação). Como pode ser visto na Figura 5, tal modelo conta com um quarto sítio de interação (virtual e indicado por M), ao qual sempre é atribuída uma carga negativa com o objetivo de melhorar a distribuição eletrostática ao redor da molécula.
Figura
4. Lista de modelos
de água presentes na versão 4.6.7 de GROMACS.
Figura
5. Esquematização
dos modelos TIP3P e TIP4P de solvatação.
- Cargas em aminoácidos: atribua cargas aos aminoácidos exatamente da forma como indicado na Tabela 1.
Tabela 1. Valores de carga a serem atribuídos a cada
aminoácido.
Aminoácido
|
Quantidade
|
Carga
|
Lisina
|
4
|
+1
|
Arginina
|
3
|
+1
|
Glutamina
|
1
|
+1
|
Ácido
aspártico
|
2
|
-1
|
- Construção de pontes dissulfeto: construa todas as pontes dissulfeto, escolhendo y para as três opções.
Note que ao final também é gerado posre.itp, arquivo
que contém a constante de força necessária para restringir a movimentação dos
átomos durante o equilíbrio (útil para a dinâmica de restrição).
Caixa
de simulação
O próximo passo
consiste em criar o ambiente de simulação, ou seja, a caixa na qual serão
posicionadas a proteína e as moléculas de água de solvatação. Para isso,
digite:
editconf
–f proteina01.gro -o proteina02.gro –bt
cubic –c –d 1.2
A
linha de comando informa que o programa editconf toma como base as coordenadas
fornecidas pelo arquivo de entrada (proteina01.gro) para criar um novo arquivo
de estrutura (proteina02.gro), no qual a proteína é centralizada (opção –c) numa caixa do tipo cúbica. O tamanho
desta também é definido: cada lado passa ter o comprimento correspondente à
soma do diâmetro da molécula (maior distância entre dois átomos) mais duas
vezes o valor indicado por –d (distância entre o soluto e a caixa), 1.2 nm nesse
caso como esquematiza a Figura 6:
Figura 6. Esquema de como o diâmetro da molécula e a
distância indicada por –d são definidos.
O resultado final é a
caixa com as especificações mostradas na tela (Figura 7):
Figura
7. Dimensões da
caixa de simulação criada nesta etapa.
Solvatação
A solvatação ocorre por meio do preenchimento da caixa de simulação vazia, criada anteriormente, com
moléculas de água. Tais moléculas provém de um arquivo de coordenadas
específico para isso (tip4p.gro), o qual contém uma caixa de solvente (216
moléculas de água), já pré-equilibrada (a 300 K por 20 ps) e ocupando um volume
definido. O preenchimento dá-se pelo empilhamento das coordenadas das
moléculas de água. Eventualmente durante esse processo, algumas dessas
moléculas podem ser descartadas. Isso ocorre porque estas se encontram muito
próximas da molécula de proteína. De forma padrão, a distância entre qualquer
átomo de uma molécula de solvente e entre qualquer átomo de uma molécula de
soluto não pode ser menor que a soma dos raios de van der Waals desses átomos
(banco de dados vdwradii.dat). genbox faz tudo isso por meio do seguinte
comando:
genbox
–cp proteina02.gro -o proteina03.gro –p topologia.top
–cs tip4p.gro
Durante
o processamento, portanto, a molécula de proteína dentro da caixa vazia
(proteina02.gro) é tratada como o soluto (opção –cp), e a caixa de água do
arquivo tip4p.gro (disponível dentro do diretório do campo de força utilizado)
é definida como solvente (opção –cs). Note que essa escolha está de acordo com
o modelo de solvatação definido anteriormente por meio de pdb2gmx. Ao final,
um novo arquivo de coordenadas é produzido, contendo a proteína já solvatada, além do correspondente arquivo de topologia (topologia.top).
Neutralização
Para esta etapa, dois
programas serão necessários: GROMACS PreProcessor (grompp) e genion.
grompp é o
pré-processador de GROMACS e, como tal, atua no pré-processamento de um arquivo,
antes deste ser utilizado por Molecular Dynamics Run (mdrun). Sua função é garantir a validade do mesmo realizando uma
série de modificações em outros arquivos que fornecem informação para sua
construção (topologia e de parâmetros). Sendo assim, agora será necessário
incluir mais um arquivo em sua pasta atual, o arquivo de parâmetros de
simulação minimizacao.mdp.
Faça isso da seguinte forma: copie e cole os parâmetros do link anterior num
editor de texto, salve no formato de texto e depois mude sua extensão para
.mdp. Após transferir o arquivo, digite o seguinte comando:
grompp
–f minimizacao.mdp –o minimizacao.tpr –c
proteina03.gro –p topologia.top
Durante o processamento, grompp lê três arquivos de
entrada distintos: arquivo de parâmetros de simulação (minimizacao.mdp),
arquivo de coordenadas (proteina03.gro)
e seu correspondente arquivo de topologia (topologia.top). Durante a leitura,
modificações são feitas, se necessário, nos arquivos de topologia (aplicação de
restrições a ligações e ângulos de ligações, e/ou remoção de restrições e de
interações envolvendo sítios virtuais) e de parâmetros (para garantir
aceleração nula durante a simulação). A leitura do arquivo de coordenadas serve
apenas para garantir que os nomes atômicos encontrados nele sejam os mesmos
lidos no arquivo de topologia (um alerta é emitido em caso contrário). Como
resultado do processamento desses arquivos, é gerado o único arquivo input
binário que mdrun utilizará (minimizacao.tpr).
Importante:
Apesar dos nomes atômicos serem lidos, apenas os correspondentes tipos atômicos
são contabilizados para gerar parâmetros de interação.
Em seguida, genion é
usado para realizar a neutralização do sistema, uma exigência quando se
trabalha num ensemble microcanônico (NVE) (a presença de um excesso de carga
num sistema desse tipo pode causar desequilíbrio da topologia proteica durante
a simulação). Em sua linha de comando, digite:
O arquivo de entrada de genion (indicado por -s)
trata-se do arquivo input de simulação criado anteriormente (minimizacao.tpr).
genion utiliza as informações contidas neste arquivo para gerar novos arquivos
de coordenadas e de topologia contendo o sistema neutralizado.
Conforme
a atribuição de carga feita a alguns aminoácidos anteriormente, o sistema
recebeu uma carga total de +6 (oito cargas positivas e duas negativas). Porém originalmente a estrutura proteica já contém certos pontos de carga, o que pode
gerar uma carga líquida diferente dessa (nosso caso, já que a estrutura possui
na realidade uma carga total de +7, como mostra a Figura 8). genion realiza a neutralização adicionando ao sistema quantidades
apropriadas (implicitamente asseguradas pela opção –neutral) de cátions e
ânions monovalentes: 21 íons cloreto (CL), por meio da flag –nname, e 14 íons
sódio (NA), por meio da flag -pname. Porém a adição desses íons implica na retirada do mesmo número de moléculas de solvente. genion realiza isso randomicamente e com a participação do usuário para a escolha do tipo de molécula a ser substituída (recomenda-se a escolha por moléculas de solvente). Portanto, digite 13.
Figura 8. Número total de cargas positivas
e negativas presentes no sistema e opções de grupo de moléculas a serem
removidas.
Além disso, mais íons
são adicionados até que a concentração especificada por –conc (0.15 mol/L) seja
atingida. O cálculo da concentração toma como base o volume da caixa do
sistema, indicado no arquivo .tpr de entrada.
3)
Dinâmica de minimização de energia
OBSERVAÇÃO: para esta etapa, serão necessários os seguintes arquivos: minimizacao.mdp e, os recém criados, proteina04.gro
e topologia.top. Dirija-se à sua pasta "minimizacao" pela linha de comando e copie-os para lá.
Aqui novamente o pré-processador concatenará
arquivos de coordenadas, topologia e de parâmetros para gerar o arquivo input
binário a ser processado por mdrun. Sua utilização se faz necessária novamente
devido às novas moléculas incluídas ou excluídas do sistema (água e íons). Para
tanto, digite:
grompp
–f minimizacao.mdp –o minimizacao.tpr –c
proteina04.gro –p topologia.top
Em seguida, o principal programa do GROMACS,
e responsável pela execução das simulações de dinâmica molecular, entra em ação.
Nesta etapa, sua função será promover o reordenamento do sistema, minimizando o
conteúdo energético do mesmo. Tal procedimento é importante, pois elimina a
possibilidade de ocorrerem os chamados contatos de van der Waals aberrantes, uma
consequência da inclusão de novas moléculas ou rompimento de ligações de
hidrogênio. Por fim, apenas digite:
mdrun –s minimizacao.tpr -v
O processamento encerra-se após a realização de
três das 200 etapas de minimização de energia, quando a força máxima aplicada
ao sistema atinge um valor menor que o limite estipulado no arquivo de
parâmetros (10000.0 kJ • mol-1 • nm-1). Em
termos de energia potencial, a energia do sistema é minimizada para um valor da
ordem de aproximadamente -1,7 x 105 kJ/mol.
Importante: de acordo com o
manual oficial do GROMACS, o parâmetro emtol possui um valor padrão de 10.0
kJ • mol-1 • nm-1, um valor muito abaixo do adotado para este
tutorial. Para atingi-lo, aumente o número de etapas. Se tal valor não for
atingido mesmo assim, estipule um valor maior, porém adequado, para emtol.
4) Dinâmica de restrição
OBSERVAÇÃO: dirija-se agora pela linha de comando à pasta "restricao" e copie nela os seguintes arquivos:
restricao.mdp (seus parâmetros são encontrados aqui),
confout.gro (produzido na etapa de minimização de energia), topologia.top e
posre.itp (produzido na fase de preparação). Se você for um usuário do cluster, crie em seu ambiente "server" um diretório específico para isso (mkdir /home/usuario/tutorial/restricao) e transfira para lá (via Filezilla) apenas o arquivo gerado a partir de grompp (minizacao.tpr), além de seu job, que pode ser encontrado aqui. Ao final, transfira todos os arquivos gerados para sua correspondente pasta em seu computador pessoal.
Novamente a partir de
grompp, concatene todos os arquivos existentes na pasta digitando:
grompp
–f restricao.mdp –o restricao.tpr
–c confout.gro –p topologia.top
Acione mdrun a partir
de (usuários do cluster fazem o mesmo por meio da submissão de um job):
mdrun –s restricao.tpr –v
Durante esta dinâmica,
são realizados os procedimentos finais para que a estrutura esteja pronta para
a principal simulação: a proteína tem seus movimentos restringidos para que as
moléculas de água (livres) possam acomodar-se ao seu redor. Após certo
intervalo de tempo, chamado tempo padrão de relaxamento (~10 ps para uma
proteína pequena como 1TSK), o sistema entra em equilíbrio. Por uma questão de
segurança, um intervalo de tempo maior é utilizado na prática (pelo menos 20 ps para este
caso). Simulações que envolvam proteínas grandes devem trabalhar com valores
ainda maiores.
5) Dinâmica molecular
OBSERVAÇÃO: em sua pasta para realizar a última e principal dinâmica deste tutorial, transfira para ela os seguintes arquivos: dinamica.mdp (clique aqui para obter seus
parâmetros), confout.gro (gerado na etapa de restrição) e topologia.top. Usuários do cluster procedem de forma análoga à etapa anterior, executando mdrun a partir da submissão de um job.
Em sua linha de comando, digite seguidamente:
grompp –f
dinamica.mdp –o dinamica.tpr –c
confout.gro –p topologia.top
mdrun –s dinamica.tpr –v
De forma análoga às etapas anteriores, grompp
prepara o único arquivo input binário a ser utilizado por mdrun. Este, por sua
vez, lê o arquivo dinamica.tpr e distribui sua topologia em processadores
paralelos se necessário. Ao final da simulação, os seguintes arquivos devem ser
produzidos:
1) [.log]: arquivo de registro;
2) [.trr]: arquivo de trajetória (contém coordenadas, velocidades e, opcionalmente, forças);
3) [.xtc]: arquivo de trajetória compacta;
4) [.gro]: arquivo estrutural (contém coordenadas e velocidades da última etapa);
5) [.edr]: arquivo de energia;
6) [.cpt]: arquivo de checkpoint.
1) [.log]: arquivo de registro;
2) [.trr]: arquivo de trajetória (contém coordenadas, velocidades e, opcionalmente, forças);
3) [.xtc]: arquivo de trajetória compacta;
4) [.gro]: arquivo estrutural (contém coordenadas e velocidades da última etapa);
5) [.edr]: arquivo de energia;
6) [.cpt]: arquivo de checkpoint.
Esses arquivos
podem agora ser utilizados para que análises sejam feitas. GROMACS possui uma
série de ferramentas destinadas a isso. Na próxima etapa, você conhecerá algumas delas.
6)
Análise dos resultados
Para comparar os
resultados obtidos com os do manual de referência, iremos trabalhar com duas
ferramentas: g_rms e g_gyrate. Além delas, também conheceremos trjconv, útil na
produção de vídeos. Para isso, tenha a disposição os seguintes arquivos em sua
pasta “analise”: dinamica.tpr e
traj.xtc (ou dinamica.xtc se você o gerar pelo cluster).
g_rms –s dinamica.tpr
–f traj.xtc –o rmsd.xvg
Quando a interface
surgir pedindo-lhe a estrutura de referência para os cálculos, escolha “Backbone”.
Para melhor visualizar o resultado, abra o arquivo num programa de plotagem de gráficos
como o Xmgrace. A aparência do gráfico deve ser como esta (Figura 9):
Figura
9. Gráficos de
RMSD vs. tempo de simulação: (a) gráfico mostrado no manual de referência e (b)
gráfico obtido por meio desta simulação.
✓g_gyrate: calcula o grau de compactação da
estrutura proteica durante a simulação (parâmetro chamado raio de giro). Em sua
linha de comando, digite:
g_gyrate –s
dinamica.tpr –f traj.xtc –o raio_giracao.xvg
Faça os cálculos
escolhendo a opção “Protein”. De forma análoga ao RMSD, abra seu arquivo num
programa de plotagem de gráficos. O resultado deve ser esse (Figura 10):
Figura
10. Gráficos de
raio de giro vs. tempo de simulação: (a) gráfico mostrado no manual de
referência e (b) gráfico obtido por meio desta simulação.
Como pode ser visto nos gráficos, a dinâmica produziu os resultados
esperados de forma muito semelhante, atestando a confiabilidade do uso deste
tutorial para a realização de simulações futuras.
Performance: ~34 ns/dia (cluster do DAFIS)
Performance: ~34 ns/dia (cluster do DAFIS)
✓ trjconv: gera um filme composto pelos frames
da simulação. Para isso, digite:
trjconv –s dinamica.tpr
–f traj.xtc –o filme.pdb –ur compact –center –pbc mol
Em um programa com
suporte a múltiplos frames (VMD), visualize o filme como a seguir:
7)
Referências
1) HESS,
B.; VAN DER SPOEL, D.; LINDAHL, E. GROMACS
User Manual (Version 4.6.7). 2014. 412 p.
2) SOARES,
R. O. S. GROMACS (4.5.5): Tutorial para
iniciantes (em Português-BR). 2013. 26 p.
3) SOARES,
R. O. S. Dinâmica Molecular de
Proteínas: estabilidade e renaturação térmica. 2009. 86 p. Dissertação (Mestrado) –
Faculdade de Ciências Farmacêuticas de Ribeirão Preto, Universidade de São
Paulo, Ribeirão Preto, 2009.
