Dinámica Molecular

Preliminares

Para llevar a cabo la dinámica molecular es necesario preparar el sistema a simular. El sistema estará representado por los siguientes dos archivos:

  1. 📃 .rst7 que contiene:
    • Las coordenadas de los átomos del sistema.
  2. 📃 .prmtop:
    • Los parámetros del campo de fuerza.
    • Información sobre la topología del sistema.
    • Propiedades de los átomos
      • Carga, masa, tipo de átomo

Información contenida en cada tipo de archivo

Los archivos .rst7 y .prmtop se construyen a partir de la información contenida en los archivos de entrada.

Estructura (Input) Parámetros (Input) Output Files
Archivo .pdb .mol2 .lib
.mol2 		.dat
.frcmod 		.prmtop
.parm7 		.inpcrd
.rst7
Moléculas ✅ ✅ ✅ — ✅ —
Nombre Átomos ✅ ✅ ✅ — ✅ —
Tipo de Átomos — ✅ ✅ ✅ ✅ —
Cargas — ✅ ✅ — ✅ —
Conexiones — ✅ ✅ — ✅ —
Coordenadas ✅ ✅ ✅ — — ✅
Masa — — — ✅ ✅ —
Parm. Enlazantes — — — ✅ ✅ —
Parm. No enlazantes — — — ✅ ✅ —

Creación de los archivos de topología .prmtop y coordenadas rst7

El programa LEaP permite utilizar la información de los archivos de entrada: estructura y parámetros de campo de fuerza, para generar los archivos .rst7 y .prmtop.

Campo de Fuerza de Amber

El objetivo de la preparación de los archivos es contar con la información necesaria respecto a las moléculas del sistema para calcular su energía potencial (\(V\)) de acuerdo con el campo de fuerza de amber durante la simulación de la dinámica molecular:

\[ \begin{aligned} V(r^{N}) & =\sum _{i\in {\text{bonds}}}{\color{red}k_{b}}_{i}(l_{i}-\color{red}l_{i}^{0})^{2} \\ & +\sum _{i\in {\text{angles}}}{\color{red}k_{a}}_{i}(\theta _{i}-\color{red}\theta _{i}^{0})^{2} \\ & +\sum _{i\in {\text{torsions}}}\sum _{n}{\frac {1}{2}}\color{red}V_{i,n}[1+\cos(\color{red}n\omega _{i}-\color{red}\gamma _{i})] \\ & +\sum _{j=1}^{N-1}\sum _{i=j+1}^{N}{\biggl \{}\color{red}\epsilon _{ij}{\biggl [}\left({\frac {\color{red}r_{ij}^{0}}{r_{ij}}} \right)^{12}-2\left({\frac {\color{red}r_{ij}^{0}}{r_{ij}}}\right)^{6}{\biggr ]}+{\frac {q_{i}q_{j}}{4\pi \epsilon _{0}r_{ij}}}{\biggr \}} \end{aligned} \]

Las constantes de fuerza (\(k_{a,b}, V_n, \epsilon\)), constantes de equilibrio (\(\theta^0, l^0, \gamma, r^0\)), y las cargas parciales (\(q\)), así como la masa de cada tipo de átomo son recuperadas de las librerías de campo de fuerza de AMBER y almacenadas en el archivo .prmtop.

Biblioteca de parámetros

  1. Proteínas y péptidos -> leaprc.protein.ff14SB:
    1. parm10.dat: Archivo principal con la mayoría de parámetros del campo de fuerza.
      1. Tipos de átomos
      2. Masas
      3. Cargas \(q_i\)
      4. Parámetros de Lennard-Jones (\(\epsilon_{ij}, r_{ij}\))
      5. Enlaces (\(k_a, l^0\))
      6. Ángulos (\(k_b, \theta^0\))
      7. Torsiones (\(V_n, n, \gamma\))
    2. frcmod.ff14SB: Suplementa a parm10.dat
    3. amino12.lib
  2. Moléculas orgánicas (general AMBER force field ) -> leaprc.gaff:
    1. gaff.dat
  3. Solvente -> leaprc.water.tip3p
    1. atomic_ions.lib
    2. solvents.lib
    3. frcmod.ionsjc_tip3p
    4. frcmod.ions234lm_126_tip3p

Estos archivos se encuentran el siguiente directorio de la máquina virtual:

/home/ssb/miniconda3/envs/amber/dat/leap/

Tutorial

Obtención y preparación de la proteína y el ligando:

Descarga de las moléculas:

  1. Definir una carpeta de trabajo para guardar los archivos que se irán generando.
mkdir wd_md
  1. Descarga los archivos de la proteína y el ligando:
  2. Ambos archivos pertenecen al modelo 4fkw.
  3. Los residuos faltantes de la proteína (loops o regiones sin estructura secundaria definida) fueron añadidos utilizando Modeller, en un proceso conocido como Model/Refine Loops. Puedes revisar más sobre el tema de modelado de proteínas en los siguientes links:
    1. Chimera Interface to Modeller
      1. Chimera Modeller tutorial
    2. Modeller Tutorial
    3. Swiss-Model -> Servidor para modelado por homología.
    4. ITASSER -> Servdor para modelado por enhebrado.
    5. Más servidores para modelado de estructura 3D de proteínas.
📂 wd_md
│  🗒 prot_unprep.pdb
│  🗒 lig_unprep.pdb

Preparación de la proteína con PDB2PQR:

  1. Tener el ambiente amber activado y localizarse en la carpeta de trabajo (wd_md).
  2. Ejecutar PDB2PQR con alguno de los siguientes comandos (dependiendo de la versión):
# Para la versión más reciente
pdb2pqr30 --ff='AMBER' --ffout='AMBER' \
    --with-ph=7.0 --drop-water --keep-chain \
    --pdb-output prot.pdb \
    prot_unprep.pdb pqr_file.pqr
# Para versiones anteriores
pdb2pqr --ff=amber --with-ph=7.0 --ffout=amber \
    --ph-calc-method=propka \
    --chain prot_unprep.pdb prot.pdb

Esto creará un nuevo archivo prot.pdb con la nomenclatura correcta en Amber y con los estados de protonación a pH 7 de los residuos ionizables.

Convertir el ligando a formato mol2

  1. Ejecutar el siguiente comando usando Open Babel:
obabel -ipdb lig_unprep.pdb -omol2 -O lig.mol2 -p 7.0

Esto creará un nuevo archivo lig.mol2 con la molécula protonizada a un pH de 7.0.

  1. En este punto, vamos a necesitar únicamente el archivo prot.pdb y el archivo lig.mol2 que ya acaban de ser creados.
📂 wd_md
│  📜 prot_unprep.pdb
│  📜 lig_unprep.pdb
|  📜 prot.pdb.propka
│  🗒 prot.pdb
│  🗒 lig.mol2

Preparación de parámetros y de coordenadas del ligando

Un tutorial más a detalle se puede consultar aquí.

Uso de Antechamber para definir la estructura y nomenclatura del ligando

  1. Desde la terminal con el ambiente de conda amber activado y estando en el directorio de trabajo (wd_md), ejecutar:
antechamber -i lig_unprep.mol2 -fi mol2 /
-o LIG_def.mol2 -fo mol2 -c bcc -s 2 -rn LIG

Esto permite usar la herramienta Antechamber para generar los archivos de topología del ligando, necesarios para la dinámica. El proceso suele tardar un poco. Una vez finalizado, se crearán varios archivos, pero de todos ellos sólo nos interesa el archivo LIG_def.mol2, que define a la molécula como un residuo y define los tipos de átomos, la nomenclatura correspondiente y las cargas parciales adecuadas para cada átomo.

El parámetro -c bcc indica a antechamber usar el modelo de cargas AM1-BCC (semi-empirical (AM1) with bond charge correction (BCC)) para calcular las cargas parciales de la molécula. La opción -rn indica renombrar el nombre de la molécula a LIG. Finalmente, -s 2define una opción de verbose que nos permite recibir información en la terminal del proceso de ejecución.

📂 wd_md
│  📜 ANTECHAMBER_*
│  📜 ATOMTYPE.INF
│  📜 sqm.*
│  📜 ...
│  🗒 prot.pdb
│  🗒 LIG_def.mol2

Uso de parmcheck para verificar los parámetros del campo de fuerza

  1. Ahora validamos que los parámetros del campo de fuerza necesarios para la molécula estén disponibles en las librerías de AMBER. Para ello ejecutar lo siguiente:
parmchk2 -i LIG_def.mol2 -f mol2 -o LIG_def.frcmod

Esto crea un archivo LIG_def.frcmod donde reporta los tipos de enlaces en la molécula que no están directamente en el campo de fuerza, pero que pueden ser representados por otros que ya se encuentran en el campo de fuerza GAFF (General AMBER Force Field). Si abres el archivo LIG_def.frcmod podrás ver ahí qué enlaces son los que requieren dicha reasignación.

General AMBER force field (GAFF) for rational drug design

  • Más información sobre los parámetros del campo de fuerza GAFF Aquí 
📂 wd_md
│  📜 ANTECHAMBER_*
│  📜 ATOMTYPE.INF
│  📜 sqm.*
│  📜 ...
│  🗒 prot.pdb
│  🗒 LIG_def.mol2
│  🗒 LIG_def.frcmod

Creación de la librería de parámetros con tleap

Los siguientes pasos son para crear una librería de parámetros del ligando que permitan construir el sistema.

  1. En la terminal, ejecutar:
tleap -f oldff/leaprc.ff14SB

Esto abrirá la consola de LEaP precargando el campo de fuerza ff14SB (más sobre los campos de fuerza de AMBER: AmberModels)

  1. Estando en la terminal de tleap, ejecutar uno a uno los siguientes comandos:
  • source leaprc.gaff
  • LIG = loadmol2 LIG_def.mol2
  • loadamberparams LIG_def.frcmod
  • saveoff LIG LIG_def.lib
  • quit

  1. Lo anterior creará la librería LIG_def.lib

  2. Revisa el archivo.

📂 wd_md
│  📜 ANTECHAMBER_*
│  📜 ATOMTYPE.INF
│  📜 sqm.*
│  📜 ...
│  🗒 prot.pdb
│  🗒 LIG_def.mol2
│  🗒 LIG_def.frcmod
│  🗒 LIG_def.lib

Construcción y solvatación del Sistema

Entra nuevamente a tleap desde la terminal, esta vez tecleando unicamente (o hay que precargar ningún campo de fuerza esta vez):

tleap

Dentro de tleap, ejecutar uno a uno los siguientes comandos:

  1. Carga de las librerías del campo de fuerza - Carga de las librerías necesarias tanto como para la proteína como para el solvente del sistema:
source leaprc.gaff
source leaprc.protein.ff14SB
source leaprc.water.tip3p
  1. Carga de las librerías del ligando - Carga de las librerías del ligando, las cuales se crearon en el paso anterior:
loadamberparams LIG_def.frcmod
loadoff LIG_def.lib
  1. Creación del complejo - Se leen los archivos pdb y mol2 de la proteína y ligando, respectivamente, y se combinan en un sólo objeto:
protein = loadpdb prot.pdb
LIG = loadmol2 LIG_def.mol2
system = combine {protein LIG}
  1. Solvatación del sistema - Se solvata el complejo proteína-ligando utilizando agua TIP3P.
solvateOct system TIP3PBOX 12
> Si se quiere solvatar en caja de agua sustituir `solvateOct` por `solvatebox`.
> El número `12` indica la distancia en A del *padding* de la caja de agua.
  1. Neutralización del sistema - Se neutraliza todo el sistema usando los iones correspondientes:
addions system Cl- 0
> Dado que en este caso el sistema tiene carga positiva, lo que hacemos en neutralizar con iones `Cl-`. El `0` indica que se añadan tantos iones como sean necesarios para alcanzar una carga total igual a 0 (neutra).
  1. Finalmente, guardar los archivos de topología (prmtop) y de coordenadas (rst7):
saveamberparm system prot_LIG.prmtop prot_LIG.rst7
quit
  1. Los archivos prot_LIG.prmtop y prot_LIG.rst7 son los que servirán como entrada para llevar a cabo la dinámica.
📂 wd_md
│  📜 ...
│  🗒 prot_LIG.prmtop
│  🗒 prot_LIG.rst7