molecular-dynamics

Dinámica Molecular

Descripción General

La simulación de dinámica molecular (MD) modela computacionalmente la evolución temporal de sistemas moleculares integrando las ecuaciones de movimiento de Newton. Esta habilidad cubre dos herramientas complementarias:

• OpenMM (https://openmm.org/): Motor de simulación MD de alto rendimiento con soporte de GPU, API en Python y soporte flexible de campos de fuerza.
• MDAnalysis (https://mdanalysis.org/): Librería de Python para leer, escribir y analizar trayectorias MD de todos los principales paquetes de simulación.

Instalación:

conda install -c conda-forge openmm mdanalysis nglview
# o
uv pip install openmm mdanalysis

Cuándo Usar Esta Habilidad

Usa dinámica molecular cuando:

• Análisis de estabilidad de proteínas: ¿Cómo afecta una mutación a la dinámica de la proteína?
• Simulaciones de unión de fármacos: Caracteriza el modo de unión y el tiempo de residencia de un ligando.
• Muestreo conformacional: Explora la flexibilidad de proteínas y cambios conformacionales.
• Interacción proteína-proteína: Modela la dinámica de la interfaz y la energética de unión.
• Análisis RMSD/RMSF: Cuantifica fluctuaciones estructurales a partir de una estructura de referencia.
• Estimación de energía libre: Calcula la energía libre de unión o la energía libre conformacional.
• Simulaciones de membrana: Modela proteínas en bicapas lipídicas.
• Proteínas intrínsecamente desordenadas: Estudia conjuntos conformacionales de IDR.

Flujo de Trabajo Principal: Simulación con OpenMM

1. Preparación del Sistema

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

def preparar_sistema_desde_pdb(archivo_pdb, nombre_campo_fuerza="amber14-all.xml",
                              modelo_agua="amber14/tip3pfb.xml"):
    """
    Prepara un sistema OpenMM a partir de un archivo PDB.

    Args:
        archivo_pdb: Ruta al archivo PDB limpio (usa PDBFixer para archivos PDB crudos).
        nombre_campo_fuerza: Archivo XML del campo de fuerza.
        modelo_agua: Archivo XML del modelo de agua.

    Returns:
        pdb, campo_de_fuerza, sistema, topología
    """
    # Cargar PDB
    pdb = PDBFile(archivo_pdb)

    # Cargar campo de fuerza
    forcefield = ForceField(nombre_campo_fuerza, modelo_agua)

    # Agregar hidrógenos y solvatar
    modeller = Modeller(pdb.topology, pdb.positions)
    modeller.addHydrogens(forcefield)

    # Agregar caja de solvente (relleno de 10 Å, NaCl 150 mM)
    modeller.addSolvent(
        forcefield,
        model='tip3p',
        padding=10*angstroms,
        ionicStrength=0.15*molar
    )

    print(f"Sistema: {modeller.topology.getNumAtoms()} átomos, "
          f"{modeller.topology.getNumResidues()} residuos")

    # Crear sistema
    system = forcefield.createSystem(
        modeller.topology,
        nonbondedMethod=PME,         # Malla de partículas Ewald para electrostática de largo alcance
        nonbondedCutoff=1.0*nanometer,
        constraints=HBonds,           # Restringir enlaces de hidrógeno (permite paso de tiempo de 2 fs)
        rigidWater=True,
        ewaldErrorTolerance=0.0005
    )

    return modeller, system

2. Minimización de Energía

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

def minimizar_energia(modeller, system, salida_pdb="minimized.pdb",
                     max_iteraciones=1000, tolerancia=10.0):
    """
    Minimiza energéticamente el sistema para eliminar choques estéricos.

    Args:
        modeller: Objeto Modeller con topología y posiciones.
        system: Sistema OpenMM.
        salida_pdb: Ruta para guardar la estructura minimizada.
        max_iteraciones: Pasos máximos de minimización.
        tolerancia: Criterio de convergencia en kJ/mol/nm.

    Returns:
        objeto simulation con posiciones minimizadas.
    """
    # Configurar integrador (no importa para la minimización)
    integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)

    # Crear simulación
    # Usar GPU si está disponible (CUDA o OpenCL), de lo contrario usar CPU
    try:
        platform = Platform.getPlatformByName('CUDA')
        properties = {'DeviceIndex': '0', 'Precision': 'mixed'}
    except Exception:
        try:
            platform = Platform.getPlatformByName('OpenCL')
            properties = {}
        except Exception:
            platform = Platform.getPlatformByName('CPU')
            properties = {}

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

    # Verificar energía inicial
    state = simulation.context.getState(getEnergy=True)
    print(f"Energía inicial: {state.getPotentialEnergy()}")

    # Minimizar
    simulation.minimizeEnergy(
        tolerance=tolerancia*kilojoules_per_mole/nanometer,
        maxIterations=max_iteraciones
    )

    state = simulation.context.getState(getEnergy=True, getPositions=True)
    print(f"Energía minimizada: {state.getPotentialEnergy()}")

    # Guardar estructura minimizada
    with open(salida_pdb, 'w') as f:
        PDBFile.writeFile(simulation.topology, state.getPositions(), f)

    return simulation

3. Equilibración NVT

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

def ejecutar_equilibracion_nvt(simulation, n_pasos=50000, temperatura=300,
                            intervalo_reporte=1000, prefijo_salida="nvt"):
    """
    Equilibración NVT: N, V, T constantes.
    Equilibra las velocidades a la temperatura objetivo.

    Args:
        simulation: Simulación OpenMM (después de minimización).
        n_pasos: Número de pasos MD (50000 × 2 fs = 100 ps).
        temperatura: Temperatura en Kelvin.
        intervalo_reporte: Pasos entre reportes de datos.
        prefijo_salida: Prefijo de archivos para trayectoria y registro.
    """
    # Agregar restricciones de posición para el esqueleto durante NVT
    # (Opcional: restringir átomos pesados)

    # Fijar temperatura
    simulation.context.setVelocitiesToTemperature(temperatura*kelvin)

    # Agregar reporters
    simulation.reporters = []

    # Archivo de registro
    simulation.reporters.append(
        StateDataReporter(
            f"{prefijo_salida}_log.txt",
            intervalo_reporte,
            step=True,
            potentialEnergy=True,
            kineticEnergy=True,
            temperature=True,
            volume=True,
            speed=True
        )
    )

    # Trayectoria DCD (formato binario compacto)
    simulation.reporters.append(
        DCDReporter(f"{prefijo_salida}_traj.dcd", intervalo_reporte)
    )

    print(f"Ejecutando equilibración NVT: {n_pasos} pasos ({n_pasos*2/1000:.1f} ps)")
    simulation.step(n_pasos)
    print("Equilibración NVT completada")

    return simulation

4. Equilibración NPT y Producción

def ejecutar_produccion_npt(simulation, n_pasos=500000, temperatura=300, presion=1.0,
                        intervalo_reporte=5000, prefijo_salida="npt"):
    """
    Ejecución de producción NPT: N, P, T constantes.

    Args:
        n_pasos: Pasos de producción (500000 × 2 fs = 1 ns).
        temperatura: Temperatura en Kelvin.
        presion: Presión en bar.
        intervalo_reporte: Pasos entre reportes.
    """
    # Agregar barostato Monte Carlo para control de presión
    system = simulation.context.getSystem()
    system.addForce(MonteCarloBarostat(presion*bar, temperatura*kelvin, 25))
    simulation.context.reinitialize(preserveState=True)

    # Actualizar reporters
    simulation.reporters = []
    simulation.reporters.append(
        StateDataReporter(
            f"{prefijo_salida}_log.txt",
            intervalo_reporte,
            step=True,
            potentialEnergy=True,
            temperature=True,
            density=True,
            speed=True
        )
    )
    simulation.reporters.append(
        DCDReporter(f"{prefijo_salida}_traj.dcd", intervalo_reporte)
    )

    # Guardar puntos de control
    simulation.reporters.append(
        CheckpointReporter(f"{prefijo_salida}_checkpoint.chk", 50000)
    )

    print(f"Ejecutando producción NPT: {n_pasos} pasos ({n_pasos*2/1000000:.2f} ns)")
    simulation.step(n_pasos)
    print("MD de producción completada")
    return simulation

Análisis de Trayectorias con MDAnalysis

1. Cargar Trayectoria

import MDAnalysis as mda
from MDAnalysis.analysis import rms, align, contacts
import numpy as np
import matplotlib.pyplot as plt

def cargar_trayectoria(archivo_topologia, archivo_trayectoria):
    """
    Carga una trayectoria MD con MDAnalysis.

    Args:
        archivo_topologia: Archivo PDB, PSF u otro archivo de topología.
        archivo_trayectoria: Trayectoria DCD, XTC, TRR u otro formato.
    """
    u = mda.Universe(archivo_topologia, archivo_trayectoria)
    print(f"Universo: {u.atoms.n_atoms} átomos, {u.trajectory.n_frames} cuadros")
    print(f"Rango de tiempo: 0 a {u.trajectory.totaltime:.0f} ps")
    return u

2. Análisis RMSD

def calcular_rmsd(u, seleccion="backbone", cuadro_referencia=0):
    """
    Calcula el RMSD de átomos seleccionados con respecto al cuadro de referencia.

    Args:
        u: Universo MDAnalysis.
        seleccion: Cadena de selección de átomos (sintaxis MDAnalysis).
        cuadro_referencia: Índice del cuadro para la estructura de referencia.

    Returns:
        Arreglo numpy de valores (tiempo, rmsd).
    """
    # Alinear trayectoria para minimizar RMSD
    alineador = align.AlignTraj(u, u, select=seleccion, in_memory=True)
    alineador.run()

    # Calcular RMSD
    R = rms.RMSD(u, select=seleccion, ref_frame=cuadro_referencia)
    R.run()

    datos_rmsd = R.results.rmsd  # columnas: cuadro, tiempo, RMSD
    return datos_rmsd

def graficar_rmsd(datos_rmsd, titulo="RMSD en el tiempo", archivo_salida="rmsd.png"):
    """Grafica el RMSD a lo largo del tiempo de simulación."""
    fig, ax = plt.subplots(figsize=(10, 4))
    ax.plot(datos_rmsd[:, 1] / 1000, datos_rmsd[:, 2], 'b-', linewidth=0.5)
    ax.set_xlabel("Tiempo (ns)")
    ax.set_ylabel("RMSD (Å)")
    ax.set_title(titulo)
    ax.axhline(datos_rmsd[:, 2].mean(), color='r', linestyle='--',
               label=f'Media: {datos_rmsd[:, 2].mean():.2f} Å')
    ax.legend()
    plt.tight_layout()
    plt.savefig(archivo_salida, dpi=150)
    return fig

3. Análisis RMSF (Flexibilidad por Residuo)

def calcular_rmsf(u, seleccion="backbone", cuadro_inicio=0):
    """
    Calcula el RMSF por residuo (flexibilidad).

    Returns:
        Arreglos de residuos y valores RMSF.
    """
    # Seleccionar átomos
    atomos = u.select_atoms(seleccion)

    # Calcular RMSF
    R = rms.RMSF(atomos)
    R.run(start=cuadro_inicio)

    # Promediar por residuo
    residuos = []
    rmsf_por_res = []
    for res in u.select_atoms(seleccion).residues:
        atomos_res = res.atoms.intersection(atomos)
        if len(atomos_res) > 0:
            residuos.append(res.resid)
            rmsf_por_res.append(R.results.rmsf[atomos_res.indices].mean())

    return np.array(residuos), np.array(rmsf_por_res)

4. Contactos Proteína-Ligando

def analizar_contactos(u, seleccion_proteina="protein", seleccion_ligando="resname LIG",
                      radio=4.5, cuadro_inicio=0):
    """
    Rastrea contactos proteína-ligando a lo largo de la trayectoria.

    Args:
        radio: Distancia de corte de contacto en Ángstroms.
    """
    proteina = u.select_atoms(seleccion_proteina)
    ligando = u.select_atoms(seleccion_ligando)

    cuadros_contacto = []
    for ts in u.trajectory[cuadro_inicio:]:
        # Encontrar átomos de proteína dentro del radio del ligando
        distancias = contacts.contact_matrix(
            proteina.positions, ligando.positions, radio
        )
        residuos_contacto = set()
        for i in range(distancias.shape[0]):
            if distancias[i].any():
                residuos_contacto.add(proteina.atoms[i].resid)
        cuadros_contacto.append(residuos_contacto)

    return cuadros_contacto

Guía de Selección de Campos de Fuerza

Sistema	Campo de Fuerza Recomendado	Modelo de Agua
Proteínas estándar	AMBER14 (amber14-all.xml)	TIP3P-FB
Proteínas + pequeñas moléculas	AMBER14 + GAFF2	TIP3P-FB
Proteínas de membrana	CHARMM36m	TIP3P
Ácidos nucleicos	AMBER99-bsc1 o AMBER14	TIP3P
Proteínas desordenadas	ff19SB o CHARMM36m	TIP3P

Herramientas de Preparación del Sistema

PDBFixer (para archivos PDB crudos)

from pdbfixer import PDBFixer
from openmm.app import PDBFile

def arreglar_pdb(pdb_entrada, pdb_salida, ph=7.0):
    """Arregla problemas comunes de PDB: residuos faltantes, átomos, agregar H, estandarizar."""
    fixer = PDBFixer(filename=pdb_entrada)
    fixer.findMissingResidues()
    fixer.findNonstandardResidues()
    fixer.replaceNonstandardResidues()
    fixer.removeHeterogens(True)    # Eliminar agua/ligandos
    fixer.findMissingAtoms()
    fixer.addMissingAtoms()
    fixer.addMissingHydrogens(ph)

    with open(pdb_salida, 'w') as f:
        PDBFile.writeFile(fixer.topology, fixer.positions, f)

    return pdb_salida

GAFF2 para Pequeñas Moléculas (a través de OpenFF Toolkit)

# Para parametrización de ligandos, usar OpenFF toolkit o ACPYPE
# uv pip install openff-toolkit
from openff.toolkit import Molecule, ForceField as OFFForceField
from openff.interchange import Interchange

def parametrizar_ligando(smiles, nombre_ff="openff-2.0.0.offxml"):
    """Genera parámetros GAFF2/OpenFF para una pequeña molécula."""
    mol = Molecule.from_smiles(smiles)
    mol.generate_conformers(n_conformers=1)

    off_ff = OFFForceField(nombre_ff)
    interchange = off_ff.create_interchange(mol.to_topology())
    return interchange

Mejores Prácticas

• Siempre minimizar antes de MD: Las estructuras PDB crudas tienen choques estéricos.
• Equilibrar antes de la producción: NVT (50–100 ps) → NPT (100–500 ps) → Producción.
• Usar GPU: Las simulaciones son 10–100× más rápidas en GPU (CUDA/OpenCL).
• Paso de tiempo de 2 fs con restricciones HBonds: Estándar; usar 4 fs con HMR (repartición de masa de hidrógeno).
• Analizar solo la trayectoria equilibrada: Descartar el primer 20–50% como equilibración.
• Guardar puntos de control: Las ejecuciones MD pueden fallar; los puntos de control permiten reiniciar.
• Condiciones de frontera periódicas: Requeridas para sistemas solvatados.
• PME para electrostática: Más preciso que métodos de corte para sistemas cargados.

Recursos Adicionales

• Documentación de OpenMM: https://openmm.org/documentation.html
• Guía de usuario de MDAnalysis: https://docs.mdanalysis.org/
• GROMACS (motor MD alternativo): https://manual.gromacs.org/
• NAMD (alternativa): https://www.ks.uiuc.edu/Research/namd/
• CHARMM-GUI (preparación de sistemas): https://www.charmm-gui.org/