Resumo do Tutorial do Autodock

Esse texto tem como objetivo auxiliar pessoas que estejam fazendo ou que almeijam fazer o tutorial do Autodock, disponível em: http://autodock.scripps.edu/faqs-help/tutorial/using-autodock-4-with-autodocktools


Resumo - Tutorial Autodock

A docagem molecular é o processo de "encaixar" uma molécula na outra,
a molécula que "encaixa" (ligante) tem várias conformações possíveis
para o processo, o software simula automaticamente cada uma delas
(em função das configurações do seu algorítimo) e nos apresenta as
propriedades e resultados de cada conformação.

O objetivo do Autodock é simular a docagem molecular, para melhor
visualização do processo foi usada uma interface gráfica chamada
AutoDockTools (ADT).

O tutorial começa ensinando a setar o diretório onde serão manipulados
todos os arquivos de entrada e saída que serão manipulados durante o
processo, após isso ele apresenta a interface gráfica do AutoDock e
seus elementos. Antes de formatar a molécula para o AutoDock, é
necessário resolver diversos problemas potenciais, o primeiro deles é
a remoção das águas cristalográficas, resolvido simplesmente
selecionado todos os HOH da molécula e deletando.

Macromolécula aberta - HIV Protease - Código pdb: 1hsg

 

Moléculas de água selecionadas

 

Moléculas de água deletadas

Hidrogênios adicionados

 

Após esse primeiro contato o tutorial apresenta diversos elementos de
interação com a interface do ADT.

Agora chegamos ao passo 2, que consiste em preparar o ligante para
o AutoDock, aqui carregamos o ligante no software, o tutorial segue
ensinando a encontrar a raiz da árvore de torsão do ligante, a qual
será destacada através de uma esfera verde, em seguida o tutorial
ensina a mostrar e modificar a quantidade de ângulos de torsão, logo
após essa formatação salvamos o arquivo no formato .pdbqt, que é o
formato de arquivo certo para o ligante.

Ligante aberto no software - No caso, o inibidor o qual será docado na macromolécula com a finalidade de quebrar sua função

 

Zoom do ligante dentro da macromolécula

 

Ligante dentro da macromolécula

Ligante (em destaque) dentro da macromolécula

 

 

 Raiz do ligante (zoom)

No passo 3 prepararemos a macromolécula, selecionando ela dentro
do software e salvando-a como um .pdbqt.

Após isso chegamos no passo 4, aqui o tutorial nos ensina a criar
uma caixa onde definimos um escopo para a busca do AutoDock, ele segue
mostrando como setar os parâmetros da caixa.
 
 Grade definida no software

No passo 5 aprendemos a preparar o arquivo de parâmetro AutoGrid,
o AutoDock não usa o receptor diretamente, em vez disso, ele usa um
conjunto de “mapas’ pré-calculados produzidos pelo AutoGrid, então nós
selecionamos o ligante e salvamos um arquivo de saída com o parâmetro
.gpf (Grid Parameter File).

O passo 6 simplesmente mostra como executar o AutoGrid. Também é possível executá-lo através do terminal usando o mesmo comando citado no tutorial:



Após isso chegamos no passo 7, onde é ensinado a preparar o
arquivo de parâmetro do AutoDock. Nesse passo é ensinado a configurar
o algorítimo genético o qual é usado na simulação da docagem, aqui
será gerado um arquivo de saída .dpf (Docking Parameter File) o qual
conterá as informações configuradas para o algorítimo genético.

Exemplo de arquivo .dpf gerado

É importante notar a linha em destaque, rmstol. Esse valor é o que irá definir a quantidade de clusters que o Autodock irá resultar (será explicado posteriormente), o valor padrão (2.0) é muito grande para um ligante tão pequeno, então é recomendado ajustar para um valor mais baixo, no caso usei 0.5.

Enfim chegamos na parte de rodar o AutoDock, no passo oito é
explicado quais arquivos de entrada e qual comando é usado para a
execução do AutoDock. Também é possível executar o AutoDock através do terminal, usando o mesmo comando citado no tutorial:





O arquivo de saída final é um .dlg, que contém todas as iterações do
algorítimo genético, o tutorial também ensina a visualizar o resultado terminar
dentro do software o qual pode ser analisado por diversos ângulos e
formas diferentes.


É possível abrir o arquivo .dlg com um editor de texto, assim, podendo
fazer uma análise mais específica, após as iterações nos é apresentado
a tabela de histograma de cluster organizada pelo AutoDocking de acordo
com as conformações que ele encontrou, para simplificar o resultado, o
software organiza conformações parecidas em clusters (considerando essas
conformações como apenas um valor), na última coluna ele apresenta um
histograma feito de hashtags representando graficamente a quantidade de
conformações dentro de cada cluster:

(Para facilitar a localização dessa representação, é possível apertar CTRL+F no documento e digitar "####")

Abaixo a representação gráfica dos clusters por um histograma.



A próxima tabela é a RMSD, que é um valor médio para o desvio médio de
átomos de uma estrutura X, relativamente a uma estrutura Y.

A etapa final é a escolha da melhor conformação possível do ligante. Para isso é necessário levar em conta diversos fatores, primeiramente é fundamental que tenham sido rodadas, pelo menos, 50 vezes o algorítimo e também que o rmsd esteja de acordo com sua molécula (como dito anteriormente). A próxima questão é: o número de avaliações de energia combinam com a dimensionalidade do  problema de busca? Geralmente a resposta convergirá para o cluster com menor energia de ligação. Entretanto, é necessário também levar em consideração o cluster mais populoso.

 Programas para análise de pockets, cavidades, túneis/canais e poros em proteínas

Para garantir maior proteção ao sítio ativo ou seletividade com relação ao substrato, muitas enzimas possuem suas cavidades catalíticas profundamente escondidas, acessíveis apenas por meio de túneis que conectam tais cavidades ao meio externo. Outras proteínas, como as proteínas membranares, possuem como função realizar o transporte de pequenas moléculas ou íons através de membranas, conduzindo-os até seu destino final por meio de canais. O estudo desses e outros tipos de espaços vazios podem ajudar na elucidação dos mecanismos de atividade biológica ou acessibilidade às proteínas. Softwares como Caver Analyst e Mole permitem que forma, tamanho, propriedades físico-químicas e até a dinâmica de pockets, cavidades, túneis/canais, e poros (figura 1) possam ser analisados.

Figura 1. Tipos de espaços vazios existentes em proteínas: A - pockets, B - cavidades, C - túneis/canais, D - poros.

Definição dos tipos de espaços vazios

- Pocket: depressão na superfície da proteína que serve como sítio de ligação para ligantes ou outras biomoléculas.

- Cavidade: espaço vazio dentro de uma proteína. Pode ocorrer de forma totalmente isolada ou não, possuindo comunicação com o meio externo por meio de túneis/canais. Geralmente constitui-se como o sítio ativo de uma enzima.

- Túneis/canais: espaço vazio que conecta o meio externo ao interior de uma cavidade.

- Poros: canal que transpassa a proteína de um ponto a outro, muito frequente em proteínas de membrana.

A seguir é apresentado um tutorial em três etapas para mostrar como dar os primeiros passos na busca por túneis e afins utilizando-se os softwares Caver Analyst (versão 1.0) e Mole (versão 2.13.9.6). Como exemplo de aplicação, a estrutura utilizada será a da enzima CYP1A P450 (código PDB: 4I8V), presente em mamíferos e pertencente à subfamília das enzimas citocromo P450. Essas enzimas constituem um grupo bastante grande e diverso de proteínas responsáveis por catalisar a reação de oxidação de compostos endógenos e exógenos (estranhos à vida), geralmente a partir de uma reação de monooxigenase. Estudo recente identificou um conjunto de canais de acesso ao sítio ativo (quatro ao todo), os quais, por possuírem distintas características entre si, acabam influenciando diretamente na maneira da enzima controlar suas especificidades com relação a diferentes substratos.

Tutorial - Mole

Por apresentar uma operacionalidade mais simples, opte por começar sua busca por espaços vazios a partir de Mole. Este possui uma interface de interação bastante autoexplicativa (posicionando o mouse sobre um parâmetro, surge na tela uma explicação sobre o significado dele) e que ajuda o usuário já ir se familiarizando com a linguagem/nomenclatura do assunto. 

Mole também está disponível na internet. Este tutorial focará sua explicação no programa, mas praticamente tudo o que será visto aqui pode ser feito da mesma forma no servidor online. No entanto, além de você ter o Java instalado em seu computador, as configurações de permissão dele devem ser modificadas. Do contrário, problemas com relação à segurança do computador impedirão que o servidor funcione corretamente. Em “Painel de Controle Java” e dentro da aba “Segurança”, abaixe o nível de segurança para um inferior (médio). Se esta opção não estiver disponível em seu computador, outra maneira de habilitar Mole consiste em introduzir o site de seu servidor na lista de sites permitidos. Ainda na aba “Segurança”, edite a seção “Lista de Exceções de Sites”, introduzindo o site em questão, como mostra a figura 2.

    

                                 Figura 2. Configurações do “Painel de Controle Java”.

Etapa 1: Carregando uma estrutura

Ao abrir o programa, o qual possui ícone de abertura dentro da pasta criada pelo processo de instalação, a tela inicial (figura 3) oferece três possibilidades de escolha. Geralmente a abertura de estruturas já processadas por outros programas ocorre com inúmeros erros de segmentação, deformando-as. Portanto, pelo menos por enquanto, opte por baixar a estrutura direto a partir do portal Protein Data Bank. No lugar da opção sugerida pelo programa (1TQN), digite “4I8V” e clique em “Download”.

                                             Figura 3. Interface de abertura de Mole.

Na tela de visualização, quatro cadeias aparecerão, conforme também indicado na primeira seção da aba “Refinement”, denominada “Chains”. Por enquanto, trabalhe apenas com a cadeia “A”, desmarcando as demais opções.

Etapa 2: Localizando espaços vazios

Automaticamente, ao carregar a estrutura, o programa já é capaz de detectar algumas cavidades. Ao todo, só a cadeia “A” possui doze cavidades (além de mais nove cavidades internas) a partir dos parâmetros iniciais sugeridos pelo programa na aba “Refinement”, à direita na tela. Para detectar outras cavidades ou ainda mais tipos de espaços vazios, um ajuste nos parâmetros é necessário. Experimente variar os parâmetros das seções “Cavity Parameters” e “Tunel Parameters” de forma aleatória até achar os espaços vazios de interesse. Resíduos na seção “Active Residues” também podem ser marcados ou desmarcados para filtrar a pesquisa. A mudança de parâmetros imediatamente inicia um novo processo de busca e os resultados aparecem em questão de segundos na aba “Results”, logo abaixo.

Porém, uma maneira mais eficiente de fazer isso é indicar um ponto de partida. Como já citado anteriormente, a estrutura 4I8V possui ao todo quatro canais que permitem que diferentes tipos de substratos adentrem a cavidade catalítica da enzima. De acordo com a literatura, um desses túneis (saída entre as hélices E, F e I e a folha beta-5) permite que moléculas de solvente consigam atingir a cavidade catalítica de CYP1A P450. Ao tomar inicialmente um ponto entre alguns resíduos convenientes (GLU 226, GLY 225, LEU 496, THR 497 e MET 498) em “Residues” da seção “Specific Point/Residue(s)”, o túnel (também conhecido como túnel S) é localizado, aparecendo na tela (destacado em laranja), como mostra a figura 4.

Figura 4. Resultado final do procedimento de procura de túneis a partir de um ponto inicial.

Por fim, uma forma ainda mais fácil de achar túneis seria a partir do próprio programa via um processo automático. Ao clicar em “Auto” do grupo “Tunnels”, o programa é capaz de achar automaticamente os túneis mais prováveis.

Etapa 3: Analisando as propriedades

Após localizar os túneis ou outros espaços vazios de interesse, chegou a hora de analisá‑los. A partir da aba “Results”, ao lado de cada espaço vazio identificado, aparecem informações que permitirão conhecê-los melhor. Clique em “Details” para ter acesso a essas informações. Tomando o túnel S como exemplo, tal janela traz o perfil do túnel, mostrando como o raio do mesmo varia ao longo de toda a sua extensão. Clicando em “Lining and Properties” na mesma janela, diversas propriedades físico-químicas do túnel são mostradas (hidrofobicidade, polaridade, mutabilidade, etc.), sendo toda a informação detalhada entre os diversos segmentos de resíduos que compõem o túnel.

A                                                B

Figura 5. Janelas abertas a partir de “Details”, trazendo o perfil de variação do raio do túnel S em função da distância (A) e as propriedades físico-químicas do mesmo (B).

Um pouco sobre Caver Analyst

Agora que você já tem alguma noção de como achar túneis e afins, realizar o mesmo tipo de busca a partir de Caver Analyst torna-se muito mais fácil. Ao contrário do que ocorre em Mole, uma estrutura já previamente processada em outro programa pode ser aberta sem problemas. Caver Analyst possui abas específicas para computar túneis e cavidades (poros ou pockets são exclusivos de Mole). A partir das abas “Tunnel” e “Cavities”, janelas podem ser abertas para a escolha dos parâmetros de busca. Alguns termos diferentes aparecem, mas são detalhados no manual do programa (disponível nas referências abaixo).

O grande diferencial do programa é que ele também pode efetuar busca por túneis ou cavidades a partir de simulações de dinâmica molecular ou um conjunto menor de frames de uma simulação. Dependendo da configuração de seu computador, o programa não é capaz de abrir trajetórias muito grandes. Por isso, na dúvida, selecione um trecho mais interessante de sua trajetória (contendo menos frames) para que o programa possa funcionar de maneira mais rápida.

Ao contrário de Mole, Caver Analyst não contabiliza propriedades físico-químicas, mas alguma análise pode ser feita a partir da geração de gráficos que são capazes de mostrar como diferentes aspectos variam ao longo da extensão do túnel. Uma opção particular do programa são os mapas de calor (figura 6), os quais mostram como o raio de um túnel varia ao longo de toda a extensão do mesmo durante uma trajetória.

Figura 6. Exemplo de mapa de calor gerado em Caver Analyst.

Referências

1)   URBAN, F.; TRUAN, G.; POMPON, D. Access channels to the buried active site control substrate specificityin CYP1A P450 enzymes. Biochimica et Biophysica Acta, 2014.

2)      CAVER ANALYST USER GUIDE – VERSION 1.0 BETA. 43 p.

3)      MOLE 2.0 – USER MANUAL. 39 p. 

TUTORIAL - GROMACS (VERSÃO 4.6.7)

SIMULAÇÃO DE UMA PEQUENA MOLÉCULA EM SOLVENTE NÃO AQUOSO

Passo 1: Preparar um arquivo em formato .pdb da estrutura desejada

A partir de qualquer programa que permita editar moléculas (Discovery Studio, Avogadro, etc.), desenhe a estrutura de interesse e salve-a em formato .pdb.

Passo 2: Obter arquivos .pdb de acordo com o campo de forças OPLS-AA

Servidor LigParGen: carregue o arquivo .pdb criado e obtenha os arquivos .itp e .pdb a partir dos parâmetros do campo de forças OPLSA-AA. 

Servidor Virtual Chemistry: obtenha a caixa de simulação do solvente escolhido (formato .pdb) e seu correspondente arquivo .itp, ambos de acordo com o campo de forças OPLS-AA.

Passo 3: Montar arquivo de topologia da molécula (3-aminopropanol):

Em um editor de textos apropriado (Notepad++), monte o arquivo de topologia como mostrado a seguir. Salve o arquivo no formato apropriado (.top).

; The force field files to be included

#include "oplsaa.ff/forcefield.itp"

; Include H9 topology

#include "H9_E30742.itp"

; Include trichloromethane topology

#include "trichloromethane.itp"

[ system ]

H9 in trichloromethane

[ molecules ]

;molecule name number

H9               1

Passo 4: Fazer a simulação

Preparação

a) Caixa de simulação

editconf -f H9_E30742-01.pdb -o H9_E30742-02.pdb -bt cubic -c -d 1.2

b) Solvatação

genbox -cp H9_E30742-02.pdb -o H9_E30742-03.pdb -p topologia.top -cs trichloromethane_T298.15.pdb

Após isso é necessário atualizar o arquivo de topologia adicionando uma linha com o número de moléculas de solvente que o programa adicionou na caixa.

; The force field files to be included

#include "oplsaa.ff/forcefield.itp"

; Include H9 topology

#include "H9_E30742.itp"

; Include trichloromethane topology

#include "trichloromethane.itp"

[ system ]

H9 in trichloromethane

[ molecules ]

;molecule name number

H9               1

trichloromethane 129

Minimização de energia

grompp -f minimizacao.mdp -o minimizacao.tpr -c H9_E30742-03.pdb -p topologia.top

mdrun -s minimizacao.tpr -v

Restrição NVT

grompp -f restricaoNVT.mdp -o restricaoNVT.tpr -c confout.gro -p topologia.top

mdrun -s restricaoNVT.tpr -v

Restrição NPT

grompp -f restricaoNPT.mdp -o restricaoNPT.tpr -c confout.gro -p topologia.top -t state.cpt

mdrun -s restricaoNPT.tpr -v

Simulação

grompp -f dinamica.mdp -o dinamica.tpr -c confout.gro -p topologia.top -t state.cpt

mdrun -s dinamica.tpr -v

Análise

Etapa 1: manter molécula inteira e centralizada

trjconv -f dinamica.trr -s dinamica.tpr -o md1.xtc -pbc whole -center

Grupo para centralização: 3-aminopropanol (H9)

Grupo para output: todo o sistema (System)

Etapa 2: remover translação e rotação em torno do centro de massa

trjconv -f md1.xtc -s dinamica.tpr -o md2.xtc -fit rot+trans

Quando solicitado pelo programa, escolha as mesmas opções da etapa anterior.

Visualização no VMD: carregue md2.xtc sobre confout.gro

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). 

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. Para os resíduos terminais, escolha a opção 0 (NH3+) para valina e 2 (COOH) para cisteína (não é recomendada a forma zwiteriônica dos resíduos terminais para um POLIpeptídeo). Dessa forma, a carga resultante adicionada à estrutura de 1TSK será de +6.

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:

  genion –s  minimizacao.tpr  –o  proteina04.gro  –p  topologia.top  –pname NA –nname CL –neutral –conc 0.15

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 10⁵ 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.

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: calcula como o Root Mean Square Deviation (RMSD) varia com o tempo. Através desse parâmetro, pode-se avaliar o quanto uma proteína deforma-se em relação a sua estrutura inicial durante a simulação. Para acioná-lo, digite:

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)

✓ 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.