Dynamique Moléculaire

Que fait Dynamique Moléculaire ?

Aperçu

La simulation de dynamique moléculaire (DM) modélise informatiquement l'évolution temporelle de systèmes moléculaires en intégrant les équations du mouvement de Newton. Cette compétence couvre deux outils complémentaires :

• OpenMM (https://openmm.org/) : Moteur de simulation DM haute performance avec support GPU, API Python et prise en charge flexible des champs de force
• MDAnalysis (https://mdanalysis.org/) : Bibliothèque Python pour lire, écrire et analyser les trajectoires DM de tous les principaux packages de simulation

Installation :

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

Quand utiliser cette compétence

Utilisez la dynamique moléculaire lorsque :

• Analyse de stabilité protéique : Comment une mutation affecte-t-elle la dynamique d'une protéine ?
• Simulations de liaison de médicaments : Caractériser le mode de liaison et le temps de résidence d'un ligand
• Échantillonnage conformational : Explorer la flexibilité protéique et les changements conformationnels
• Interaction protéine-protéine : Modéliser la dynamique d'interface et l'énergétique de liaison
• Analyse RMSD/RMSF : Quantifier les fluctuations structurales par rapport à une structure de référence
• Estimation d'énergie libre : Calculer l'énergie libre de liaison ou l'énergie libre conformationnelle
• Simulations membranaires : Modéliser des protéines dans des bicouches lipidiques
• Protéines intrinsèquement désordonnées : Étudier les ensembles conformationnels des PID

Flux de travail de base : Simulation OpenMM

1. Préparation du système

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

def prepare_system_from_pdb(pdb_file, forcefield_name="amber14-all.xml",
                              water_model="amber14/tip3pfb.xml"):
    """
    Préparer un système OpenMM à partir d'un fichier PDB.

    Args :
        pdb_file : Chemin vers le fichier PDB nettoyé (utilisez PDBFixer pour les fichiers PDB bruts)
        forcefield_name : Fichier XML du champ de force
        water_model : Fichier XML du modèle d'eau

    Retourne :
        pdb, forcefield, system, topology
    """
    # Charger le PDB
    pdb = PDBFile(pdb_file)

    # Charger le champ de force
    forcefield = ForceField(forcefield_name, water_model)

    # Ajouter les hydrogènes et solvater
    modeller = Modeller(pdb.topology, pdb.positions)
    modeller.addHydrogens(forcefield)

    # Ajouter une boîte de solvant (rembourrage de 10 Å, 150 mM NaCl)
    modeller.addSolvent(
        forcefield,
        model='tip3p',
        padding=10*angstroms,
        ionicStrength=0.15*molar
    )

    print(f"Système : {modeller.topology.getNumAtoms()} atomes, "
          f"{modeller.topology.getNumResidues()} résidus")

    # Créer le système
    system = forcefield.createSystem(
        modeller.topology,
        nonbondedMethod=PME,         # Particle Mesh Ewald pour l'électrostatique à longue portée
        nonbondedCutoff=1.0*nanometer,
        constraints=HBonds,           # Contraindre les liaisons hydrogène (permet un pas de 2 fs)
        rigidWater=True,
        ewaldErrorTolerance=0.0005
    )

    return modeller, system

2. Minimisation d'énergie

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

def minimize_energy(modeller, system, output_pdb="minimized.pdb",
                     max_iterations=1000, tolerance=10.0):
    """
    Minimiser l'énergie du système pour éliminer les conflits stériques.

    Args :
        modeller : Objet Modeller avec topologie et positions
        system : Système OpenMM
        output_pdb : Chemin pour enregistrer la structure minimisée
        max_iterations : Nombre maximal d'étapes de minimisation
        tolerance : Critère de convergence en kJ/mol/nm

    Retourne :
        objet simulation avec positions minimisées
    """
    # Configurer l'intégrateur (n'a pas d'importance pour la minimisation)
    integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)

    # Créer la simulation
    # Utiliser le GPU si disponible (CUDA ou OpenCL), sinon 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)

    # Vérifier l'énergie initiale
    state = simulation.context.getState(getEnergy=True)
    print(f"Énergie initiale : {state.getPotentialEnergy()}")

    # Minimiser
    simulation.minimizeEnergy(
        tolerance=tolerance*kilojoules_per_mole/nanometer,
        maxIterations=max_iterations
    )

    state = simulation.context.getState(getEnergy=True, getPositions=True)
    print(f"Énergie minimisée : {state.getPotentialEnergy()}")

    # Enregistrer la structure minimisée
    with open(output_pdb, 'w') as f:
        PDBFile.writeFile(simulation.topology, state.getPositions(), f)

    return simulation

3. Équilibration NVT

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

def run_nvt_equilibration(simulation, n_steps=50000, temperature=300,
                            report_interval=1000, output_prefix="nvt"):
    """
    Équilibration NVT : N, V, T constants.
    Équilibrer les vitesses à la température cible.

    Args :
        simulation : Simulation OpenMM (après minimisation)
        n_steps : Nombre d'étapes MD (50000 × 2 fs = 100 ps)
        temperature : Température en Kelvin
        report_interval : Pas entre les rapports de données
        output_prefix : Préfixe de fichier pour la trajectoire et le journal
    """
    # Ajouter des contraintes de position pour le squelette pendant NVT
    # (Optionnel : contraindre les atomes lourds)

    # Définir la température
    simulation.context.setVelocitiesToTemperature(temperature*kelvin)

    # Ajouter des rapporteurs
    simulation.reporters = []

    # Fichier journal
    simulation.reporters.append(
        StateDataReporter(
            f"{output_prefix}_log.txt",
            report_interval,
            step=True,
            potentialEnergy=True,
            kineticEnergy=True,
            temperature=True,
            volume=True,
            speed=True
        )
    )

    # Trajectoire DCD (format binaire compact)
    simulation.reporters.append(
        DCDReporter(f"{output_prefix}_traj.dcd", report_interval)
    )

    print(f"Exécution de l'équilibration NVT : {n_steps} étapes ({n_steps*2/1000:.1f} ps)")
    simulation.step(n_steps)
    print("Équilibration NVT terminée")

    return simulation

4. Équilibration et Production NPT

def run_npt_production(simulation, n_steps=500000, temperature=300, pressure=1.0,
                        report_interval=5000, output_prefix="npt"):
    """
    Simulation de production NPT : N, P, T constants.

    Args :
        n_steps : Étapes de production (500000 × 2 fs = 1 ns)
        temperature : Température en Kelvin
        pressure : Pression en bar
        report_interval : Pas entre les rapports
    """
    # Ajouter un barostat Monte-Carlo pour le contrôle de la pression
    system = simulation.context.getSystem()
    system.addForce(MonteCarloBarostat(pressure*bar, temperature*kelvin, 25))
    simulation.context.reinitialize(preserveState=True)

    # Mettre à jour les rapporteurs
    simulation.reporters = []
    simulation.reporters.append(
        StateDataReporter(
            f"{output_prefix}_log.txt",
            report_interval,
            step=True,
            potentialEnergy=True,
            temperature=True,
            density=True,
            speed=True
        )
    )
    simulation.reporters.append(
        DCDReporter(f"{output_prefix}_traj.dcd", report_interval)
    )

    # Enregistrer des points de contrôle
    simulation.reporters.append(
        CheckpointReporter(f"{output_prefix}_checkpoint.chk", 50000)
    )

    print(f"Exécution de la production NPT : {n_steps} étapes ({n_steps*2/1000000:.2f} ns)")
    simulation.step(n_steps)
    print("Production DM terminée")
    return simulation

Analyse de trajectoire avec MDAnalysis

1. Charger la trajectoire

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

def load_trajectory(topology_file, trajectory_file):
    """
    Charger une trajectoire DM avec MDAnalysis.

    Args :
        topology_file : Fichier PDB, PSF ou autre fichier de topologie
        trajectory_file : Trajectoire DCD, XTC, TRR ou autre
    """
    u = mda.Universe(topology_file, trajectory_file)
    print(f"Univers : {u.atoms.n_atoms} atomes, {u.trajectory.n_frames} images")
    print(f"Plage de temps : 0 à {u.trajectory.totaltime:.0f} ps")
    return u

2. Analyse RMSD

def compute_rmsd(u, selection="backbone", reference_frame=0):
    """
    Calculer le RMSD des atomes sélectionnés par rapport à l'image de référence.

    Args :
        u : Univers MDAnalysis
        selection : Chaîne de sélection d'atomes (syntaxe MDAnalysis)
        reference_frame : Index de l'image pour la structure de référence

    Retourne :
        Tableau numpy de valeurs (temps, rmsd)
    """
    # Aligner la trajectoire pour minimiser le RMSD
    aligner = align.AlignTraj(u, u, select=selection, in_memory=True)
    aligner.run()

    # Calculer le RMSD
    R = rms.RMSD(u, select=selection, ref_frame=reference_frame)
    R.run()

    rmsd_data = R.results.rmsd  # colonnes : image, temps, RMSD
    return rmsd_data

def plot_rmsd(rmsd_data, title="RMSD au fil du temps", output_file="rmsd.png"):
    """Tracer le RMSD au fil du temps de simulation."""
    fig, ax = plt.subplots(figsize=(10, 4))
    ax.plot(rmsd_data[:, 1] / 1000, rmsd_data[:, 2], 'b-', linewidth=0.5)
    ax.set_xlabel("Temps (ns)")
    ax.set_ylabel("RMSD (Å)")
    ax.set_title(title)
    ax.axhline(rmsd_data[:, 2].mean(), color='r', linestyle='--',
               label=f'Moyenne : {rmsd_data[:, 2].mean():.2f} Å')
    ax.legend()
    plt.tight_layout()
    plt.savefig(output_file, dpi=150)
    return fig

3. Analyse RMSF (Flexibilité par résidu)

def compute_rmsf(u, selection="backbone", start_frame=0):
    """
    Calculer le RMSF par résidu (flexibilité).

    Retourne :
        Tableaux resids, rmsf_values
    """
    # Sélectionner les atomes
    atoms = u.select_atoms(selection)

    # Calculer le RMSF
    R = rms.RMSF(atoms)
    R.run(start=start_frame)

    # Moyenne par résidu
    resids = []
    rmsf_per_res = []
    for res in u.select_atoms(selection).residues:
        res_atoms = res.atoms.intersection(atoms)
        if len(res_atoms) > 0:
            resids.append(res.resid)
            rmsf_per_res.append(R.results.rmsf[res_atoms.indices].mean())

    return np.array(resids), np.array(rmsf_per_res)

4. Contacts Protéine-Ligand

def analyze_contacts(u, protein_sel="protein", ligand_sel="resname LIG",
                      radius=4.5, start_frame=0):
    """
    Suivre les contacts protéine-ligand au fil de la trajectoire.

    Args :
        radius : Distance de contact seuil en Angströms
    """
    protein = u.select_atoms(protein_sel)
    ligand = u.select_atoms(ligand_sel)

    contact_frames = []
    for ts in u.trajectory[start_frame:]:
        # Trouver les atomes de la protéine dans le rayon du ligand
        distances = contacts.contact_matrix(
            protein.positions, ligand.positions, radius
        )
        contact_residues = set()
        for i in range(distances.shape[0]):
            if distances[i].any():
                contact_residues.add(protein.atoms[i].resid)
        contact_frames.append(contact_residues)

    return contact_frames

Guide de sélection du champ de force

Système	Champ de force recommandé	Modèle d'eau
Protéines standard	AMBER14 (amber14-all.xml)	TIP3P-FB
Protéines + petites molécules	AMBER14 + GAFF2	TIP3P-FB
Protéines membranaires	CHARMM36m	TIP3P
Acides nucléiques	AMBER99-bsc1 ou AMBER14	TIP3P
Protéines désordonnées	ff19SB ou CHARMM36m	TIP3P

Outils de préparation du système

PDBFixer (pour fichiers PDB bruts)

from pdbfixer import PDBFixer
from openmm.app import PDBFile

def fix_pdb(input_pdb, output_pdb, ph=7.0):
    """Corriger les problèmes courants de PDB : résidus manquants, atomes, ajout H, normalisation."""
    fixer = PDBFixer(filename=input_pdb)
    fixer.findMissingResidues()
    fixer.findNonstandardResidues()
    fixer.replaceNonstandardResidues()
    fixer.removeHeterogens(True)    # Supprimer eau/ligands
    fixer.findMissingAtoms()
    fixer.addMissingAtoms()
    fixer.addMissingHydrogens(ph)

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

    return output_pdb

GAFF2 pour petites molécules (via OpenFF Toolkit)

# Pour la paramétrisation des ligands, utilisez OpenFF toolkit ou ACPYPE
# uv pip install openff-toolkit
from openff.toolkit import Molecule, ForceField as OFFForceField
from openff.interchange import Interchange

def parameterize_ligand(smiles, ff_name="openff-2.0.0.offxml"):
    """Générer des paramètres GAFF2/OpenFF pour une petite molécule."""
    mol = Molecule.from_smiles(smiles)
    mol.generate_conformers(n_conformers=1)

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

Bonnes pratiques

• Toujours minimiser avant MD : Les structures PDB brutes présentent des conflits stériques
• Équilibrer avant la production : NVT (50–100 ps) → NPT (100–500 ps) → Production
• Utiliser le GPU : Les simulations sont 10–100× plus rapides sur GPU (CUDA/OpenCL)
• Pas de 2 fs avec contraintes HBonds : Standard ; utilisez 4 fs avec HMR (redistribution de masse des hydrogènes)
• Analyser uniquement la trajectoire équilibrée : Éliminez les premiers 20–50 % comme équilibration
• Enregistrer des points de contrôle : Les runs MD peuvent échouer ; les points de contrôle permettent le redémarrage
• Conditions aux limites périodiques : Requises pour les systèmes solvatés
• PME pour l'électrostatique : Plus précis que les méthodes de coupure pour les systèmes chargés

Ressources supplémentaires

• Documentation OpenMM : https://openmm.org/documentation.html
• Guide utilisateur MDAnalysis : https://docs.mdanalysis.org/
• GROMACS (moteur MD alternatif) : https://manual.gromacs.org/
• NAMD (alternatif) : https://www.ks.uiuc.edu/Research/namd/