F09 La siguiente pregunta: el acoplamiento es una foto fija

En F09, se calculó con AutoDock Vina cómo se une el ligando al sitio de unión en cuanto a su orientación. Sin embargo, ese resultado es, estrictamente hablando, una sola foto fija. En la realidad, las proteínas y los ligandos están constantemente vibrando a temperatura corporal, no quedan completamente fijos tras la unión y algunas interacciones se forman y rompen repetidamente en intervalos breves. Aquí surge una pregunta que los modelos de predicción estructural y acoplamiento tratados en F02～F08 no logran responder por igual: cómo se mueve esta estructura a lo largo del tiempo.

La Dinámica Molecular (Molecular Dynamics, MD) responde directamente a esta pregunta utilizando la herramienta más fundamental de la física: las ecuaciones de movimiento de Newton. GROMACS y OpenMM son los dos motores representativos que se utilizan en la práctica para realizar estos cálculos.

Principio — Aplicar las ecuaciones de Newton átomo por átomo

Campo de fuerzas (Force Field): las reglas que definen la fuerza entre átomos

El punto de partida de la MD es calcular la fuerza que actúa entre todos los pares de átomos. El conjunto de reglas que definen esta fuerza se denomina campo de fuerzas. Aunque tiene una composición similar a los términos de van der Waals y electrostáticos del acoplamiento tratado en F09, el campo de fuerzas de MD también incluye términos internos al enlace, como la longitud del enlace, el ángulo de enlace y el ángulo diedro.

$U_{\text{total}} = \underbrace{\sum_{\text{bonds}} k_b(r-r_0)^2}_{\text{estiramiento del enlace}} + \underbrace{\sum_{\text{angles}} k_\theta(\theta-\theta_0)^2}_{\text{flexión del ángulo de enlace}} + \underbrace{\sum_{\text{dihedrals}} k_\phi[1+\cos(n\phi-\delta)]}_{\text{rotación del ángulo diedro}} + \underbrace{\sum_{i<j}\left[\frac{A_{ij}}{r_{ij}^{12}} - \frac{B_{ij}}{r_{ij}^{6}} + \frac{q_iq_j}{4\pi\epsilon_0 r_{ij}}\right]}_{\text{interacciones no enlazadas}}$

Los tres primeros términos aproximan la vibración y rotación de vecinos conectados por enlaces químicos mediante un modelo de resorte, mientras que el último término corresponde a las interacciones no enlazadas, análogas al término de Lennard-Jones + Coulomb de F09. Si has oído hablar de nombres como AMBER o CHARMM, esos son conjuntos específicos de campos de fuerzas donde los parámetros $k_b, k_\theta, A_{ij}, B_{ij}$ de estos términos se ajustan con precisión mediante experimentos y cálculos de química cuántica.

Ecuaciones de Newton e integración numérica

La fuerza que actúa sobre cada átomo derivada de esta energía potencial $U$ es el negativo del gradiente de la energía.

$F_i = -\nabla_i U, \qquad m_i \frac{d^2 r_i}{dt^2} = F_i$

Esta es la segunda ley de Newton pura. El problema radica en que este sistema de ecuaciones diferenciales, con miles a decenas de miles de átomos entrelazados, no puede resolverse analíticamente. Por ello, se divide el tiempo en intervalos muy pequeños ($\Delta t$, típicamente 1～2 femtosegundos) y se utiliza integración numérica para actualizar la posición y la velocidad en cada paso. El método más ampliamente utilizado es la integración leapfrog o velocity Verlet.

$v\left(t+\frac{\Delta t}{2}\right) = v\left(t-\frac{\Delta t}{2}\right) + \frac{F(t)}{m}\Delta t$

$r(t+\Delta t) = r(t) + v\left(t+\frac{\Delta t}{2}\right)\Delta t$

Este método, que actualiza la velocidad con un desfase de medio paso (leapfrog, "salto de rana"), posee buenas propiedades de conservación de energía, lo que evita la acumulación de errores en simulaciones a largo plazo.

Ejemplo de cálculo manual: ¿Por qué el paso de tiempo debe estar en unidades de femtosegundos?

El período de vibración de los enlaces que involucran átomos de hidrógeno es aproximadamente de 10 femtosegundos ($10^{-14}$ segundos). Para que la integración numérica funcione de manera estable, el paso de tiempo $\Delta t$ debe ser significativamente menor que este período de vibración más rápido. Según un criterio de estabilidad aproximado, si se establece el paso de tiempo como aproximadamente 1/10 del período de vibración,

$\Delta t \lesssim \frac{10\,\text{fs}}{10} = 1\,\text{fs}$

Es decir, la razón para elegir un paso de tiempo de 1～2 femtosegundos no es arbitraria, sino una restricción física para no omitir el enlace que vibra más rápidamente en el sistema (principalmente los enlaces que involucran hidrógeno). Si se desea obtener una trayectoria de 1 nanosegundo ($10^{-9}$ segundos), el número de pasos necesarios es

$\frac{1\,\text{ns}}{1\,\text{fs}} = \frac{10^{-9}}{10^{-15}} = 10^{6}\,\text{steps}$

Es decir, 1 millón de pasos. Dado que el título de esta sección menciona "trayectorias de 100 ns", para alcanzar escalas de tiempo significativas en la práctica (decenas a cientos de nanosegundos), el número de pasos aumenta hasta el orden de los miles de millones, por lo que en este cálculo manual se puede experimentar directamente por qué la dinámica molecular es tan intensiva en términos computacionales.

Práctica: Simulación de un péptido corto con OpenMM

# Ejecutar en Colab T4
!pip install -q openmm

from openmm.app import *
from openmm import *
from openmm.unit import *

# Ejemplo: estructura cristalina pequeña crambin (1CRN) descargada realmente en la sesión de Colab
from urllib.request import urlretrieve
urlretrieve("https://files.rcsb.org/download/1CRN.pdb", "1CRN.pdb")
pdb = PDBFile("1CRN.pdb")
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')

# Añade los hidrógenos faltantes a la estructura cristalina según la plantilla del campo de fuerza.
modeller = Modeller(pdb.topology, pdb.positions)
modeller.addHydrogens(forcefield, pH=7.0)

# Simulación en vacío con fines educativos. El análisis real requiere configuración de solvente, iones y PME.
system = forcefield.createSystem(modeller.topology, nonbondedMethod=NoCutoff, constraints=HBonds)

integrator = LangevinMiddleIntegrator(
    300 * kelvin,        # Temperatura objetivo
    1 / picosecond,       # Coeficiente de fricción (intensidad del acoplamiento con el baño térmico)
    2 * femtoseconds       # Paso de tiempo — exactamente la misma escala que el cálculo manual anterior
)

simulation = Simulation(modeller.topology, system, integrator)
simulation.context.setPositions(modeller.positions)

# 1) Minimización de energía: eliminar primero las superposiciones atómicas anómalas (choque estérico) de la estructura inicial
simulation.minimizeEnergy()

# 2) Ejecución de producción breve: 10.000 pasos = 20 picosegundos (escala reducida para fines educativos)
simulation.reporters.append(
    StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True)
)
simulation.step(10000)

constraints=HBonds es un truco que fija por completo (mediante restricciones) las vibraciones de los enlaces que involucran hidrógeno. Como se vio en el cálculo manual anterior, las vibraciones más rápidas limitan el paso de tiempo; al bloquear estas vibraciones mediante restricciones, es posible aumentar el paso de tiempo hasta 2 femtosegundos, lo cual es ampliamente utilizado en la práctica.

Mapeo a CS

• Integración numérica: Los métodos de integración leapfrog/Verlet, que aproximan y resuelven ecuaciones diferenciales continuas mediante pasos de tiempo discretos, son conceptualmente idénticos a los integradores numéricos utilizados en simulaciones de cuerpos rígidos en motores físicos y gráficos por computadora.
• Exploración del espacio de estados: El proceso mediante el cual la simulación actualiza la posición y velocidad (estado) del sistema en cada paso, siguiendo una trayectoria, es estructuralmente similar a cómo un agente explora secuencialmente el espacio de estados en el aprendizaje por refuerzo.
• Aceleración basada en restricciones: El truco de fijar las vibraciones rápidas mediante restricciones para aumentar el paso de tiempo es un patrón de optimización general que consiste en aproximar y fijar los detalles finos en un grafo computacional para ahorrar operaciones totales.

Defectos comunes

• Iniciar la simulación sin minimización de energía: Si quedan superposiciones atómicas en la estructura inicial, las fuerzas se vuelven anormalmente grandes desde el primer paso, provocando que el sistema diverja (explotión). minimizeEnergy() es un paso esencial que no debe omitirse.
• Analizar directamente los resultados de la producción sin equilibración: La estructura inicial aún no se ha adaptado a la temperatura y presión objetivo. Incluir las partes iniciales de la trayectoria en el análisis, sin pasar por una fase suficiente de equilibración, introduce sesgos artificiales.
• Sacar conclusiones estadísticamente significativas a partir de simulaciones cortas: Los 20 picosegundos de este ejercicio son puramente una escala reducida con fines educativos. Para discutir la estabilidad de los enlaces o cambios estructurales reales, se requieren decenas a cientos de nanosegundos, y en algunos casos múltiples réplicas.
• Confundir el ejemplo en vacío con un entorno biológico real: El código anterior es un ejemplo mínimo que muestra la API y el flujo de integración. Los cálculos a nivel de publicación científica requieren agua e iones explícitos, condiciones de frontera periódicas, electrostática de largo alcance (PME), equilibración NVT/NPT, validación del campo de fuerzas/parámetros del ligando y réplicas independientes.

Para profundizar más

El texto ha sido reconstruido directamente por el equipo de investigación de BPD. Profundice consultando el artículo original y los recursos oficiales.

• Documentación y tutoriales oficiales de GROMACS: manual.gromacs.org — Flujo de trabajo completo para la selección del campo de fuerzas y preparación del sistema.
• Documentación oficial de OpenMM: openmm.org — API de Python y configuración de aceleración por GPU.
• Artículo sobre la fuerza de campo AMBER: Maier et al. (2015), ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters, J. Chem. Theory Comput.

Se ha completado el episodio 5 de F1.2 (predicción y acoplamiento de estructuras abiertas, MD). En los próximos cinco episodios (F11～F15), volveremos a la secuencia en sí misma, en lugar de las estructuras, para abordar cómo se utilizan los modelos de lenguaje de proteínas en el embebido, la generación y el diseño.