Metadinámica

La metadinámica (MTD; también abreviada como METAD o MetaD) es un método de simulación por computadora en física computacional , química y biología . Se utiliza para estimar la energía libre y otras funciones de estado de un sistema , donde la ergodicidad se ve obstaculizada por la forma del paisaje energético del sistema . Fue sugerida por primera vez por Alessandro Laio y Michele Parrinello en 2002 ^[1] y generalmente se aplica dentro de las simulaciones de dinámica molecular . La MTD se parece mucho a una serie de métodos más nuevos, como la dinámica molecular sesgada adaptativamente, ^[2] las fuerzas de coordenadas de reacción adaptativas ^[3] y el muestreo de paraguas de elevación local. ^[4] Más recientemente, tanto la metadinámica original como la bien templada ^[5] se derivaron en el contexto del muestreo de importancia y se demostró que eran un caso especial del entorno de potencial de sesgo adaptativo. ^[6] La MTD está relacionada con el muestreo de Wang-Landau . ^[7]

La técnica se basa en una gran cantidad de métodos relacionados que incluyen (en orden cronológico) la deflación, ^[8] la tunelización, ^[9] la búsqueda tabú , ^[10] la elevación local , ^[11] la inundación conformacional, ^[12] los métodos de Engkvist-Karlström ^[13] y Adaptive Biasing Force . ^[14]

La metadinámica se ha descrito informalmente como "llenar los pozos de energía libre con arena computacional". ^[15] El algoritmo supone que el sistema puede describirse mediante unas pocas variables colectivas (CV). Durante la simulación, se calcula la ubicación del sistema en el espacio determinado por las variables colectivas y se añade un potencial gaussiano positivo al paisaje energético real del sistema. De esta manera, se desanima al sistema a volver al punto anterior. Durante la evolución de la simulación, se suman cada vez más gaussianas, lo que desanima cada vez más al sistema a volver a sus pasos anteriores, hasta que el sistema explora el paisaje energético completo; en este punto, la energía libre modificada se convierte en una constante en función de las variables colectivas, lo que es la razón por la que las variables colectivas comienzan a fluctuar fuertemente. En este punto, el paisaje energético se puede recuperar como el opuesto de la suma de todas las gaussianas.

El intervalo de tiempo entre la adición de dos funciones gaussianas, así como la altura y el ancho de la función gaussiana, se ajustan para optimizar la relación entre precisión y costo computacional. Con solo cambiar el tamaño de la función gaussiana, la metadinámica se puede ajustar para producir muy rápidamente un mapa aproximado del paisaje energético utilizando funciones gaussianas grandes, o se puede utilizar para una descripción más detallada utilizando funciones gaussianas más pequeñas. ^[1] Por lo general, la metadinámica bien temperada ^[5] se utiliza para cambiar el tamaño de la función gaussiana de manera adaptativa. Además, el ancho de la función gaussiana se puede adaptar con la metadinámica gaussiana adaptativa. ^[16]

La metadinámica tiene la ventaja, sobre métodos como el muestreo de paraguas adaptativo , de no requerir una estimación inicial del paisaje energético a explorar. ^[1] Sin embargo, no es trivial elegir las variables colectivas adecuadas para una simulación compleja. Por lo general, se requieren varios ensayos para encontrar un buen conjunto de variables colectivas, pero existen varios procedimientos automáticos propuestos: coordenadas esenciales , ^[17] Sketch-Map, ^[18] y variables colectivas no lineales basadas en datos. ^[19]

Las simulaciones metadinámicas independientes (réplicas) se pueden acoplar entre sí para mejorar la usabilidad y el rendimiento paralelo. Se han propuesto varios métodos de este tipo: el MTD de caminantes múltiples ^[20] , el MTD de templado paralelo ^[21], el MTD de intercambio de sesgo ^[22] y el MTD de templado de variable colectiva ^[23] . Los últimos tres son similares al método de templado paralelo y utilizan intercambios de réplicas para mejorar el muestreo. Normalmente, se utiliza el algoritmo Metropolis-Hastings para los intercambios de réplicas, pero los algoritmos de intercambio infinito ^[24] y Suwa-Todo ^[25] ofrecen mejores tasas de intercambio de réplicas ^{[26] .}

Las simulaciones MTD típicas (de réplica única) pueden incluir hasta 3 CV, incluso utilizando el enfoque de réplica múltiple, es difícil superar los 8 CV en la práctica. Esta limitación proviene del potencial de sesgo, construido agregando funciones gaussianas (núcleos). Es un caso especial del estimador de densidad de núcleos (KDE). El número de núcleos requeridos, para una precisión KDE constante, aumenta exponencialmente con el número de dimensiones. Por lo tanto, la longitud de la simulación MTD tiene que aumentar exponencialmente con el número de CV para mantener la misma precisión del potencial de sesgo. Además, el potencial de sesgo, para una evaluación rápida, generalmente se aproxima con una cuadrícula regular . ^{[27] La} ​​memoria requerida para almacenar la cuadrícula también aumenta exponencialmente con el número de dimensiones (CV).

Una generalización de alta dimensión de la metadinámica es NN2B. ^[28] Se basa en dos algoritmos de aprendizaje automático : el estimador de densidad del vecino más cercano (NNDE) y la red neuronal artificial (ANN). NNDE reemplaza a KDE para estimar las actualizaciones del potencial de sesgo a partir de simulaciones sesgadas cortas, mientras que ANN se utiliza para aproximar el potencial de sesgo resultante. ANN es una representación de memoria eficiente de funciones de alta dimensión, donde las derivadas (fuerzas de sesgo) se calculan de manera efectiva con el algoritmo de retropropagación . ^{[28] [29]}

Un método alternativo, que explota las ANN para el potencial de sesgo adaptativo, utiliza fuerzas potenciales medias para la estimación. ^[30] Este método también es una generalización de alta dimensión del método de fuerza de sesgo adaptativo (ABF). ^[31] Además, el entrenamiento de las ANN se mejora utilizando la regularización bayesiana, ^[32] y el error de aproximación se puede inferir mediante el entrenamiento de un conjunto de ANN. ^[30]

En 2015, White, Dama y Voth introdujeron la metadinámica dirigida por experimentos, un método que permite dar forma a las simulaciones de dinámica molecular para que coincidan con una superficie de energía libre deseada . Esta técnica guía la simulación hacia conformaciones que se alinean con los datos experimentales, lo que mejora nuestra comprensión de los sistemas moleculares complejos y su comportamiento. ^[33]

En 2020, se propuso una evolución de la metadinámica, el método de muestreo mejorado con probabilidad sobre la marcha (OPES), ^{[34] [35] [36]} que ahora es el método de elección del grupo de investigación de Michele Parrinello . ^[37] El método OPES tiene solo unos pocos parámetros robustos, converge más rápido que la metadinámica y tiene un esquema de reponderación sencillo. ^[38] En 2024, se desarrolló una variante de intercambio de réplicas de OPES, llamada OneOPES, ^[39] diseñada para explotar un gradiente térmico y múltiples CV para muestrear grandes sistemas bioquímicos con varios grados de libertad. Esta variante tiene como objetivo abordar el desafío de describir dichos sistemas, donde los numerosos grados de libertad a menudo son difíciles de capturar con solo unos pocos CV. OPES se ha implementado en la biblioteca PLUMED desde la versión 2.7. ^[40]

Supongamos que tenemos un sistema clásico de partículas con posiciones en las coordenadas cartesianas . La interacción de las partículas se describe con una función potencial . La forma de la función potencial (por ejemplo, dos mínimos locales separados por una barrera de alta energía) impide un muestreo ergódico con dinámica molecular o métodos de Monte Carlo . ${\textstyle N}$${\textstyle \{{\vec {r}}_{i}\}}$ ${\textstyle (i\en 1...N)}$ ${\textstyle ({\vec {r}}_{i}\in \mathbb {R} ^{3})}$${\textstyle V\equiv V(\{{\vec {r}}_{i}\})}$

Una idea general de MTD es mejorar el muestreo del sistema desalentando la revisión de los estados muestreados. Esto se logra aumentando el hamiltoniano del sistema con un potencial de sesgo : ${\textstyle H}$${\displaystyle V_{\text{sesgo}}}$

${\displaystyle H=T+V+V_{\text{sesgo}}}$.

El potencial de polarización es una función de variables colectivas . Una variable colectiva es una función de las posiciones de las partículas . El potencial de polarización se actualiza continuamente añadiendo polarización a una tasa , donde es un valor instantáneo de la variable colectiva en el tiempo : ${\textstyle (V_{\text{sesgo}}\equiv V_{\text{sesgo}}({\vec {s}}\,))}$${\displaystyle ({\vec {s}}\equiv {\vec {s}}(\{{\vec {r}}_{i}\}))}$${\estilo de visualización \omega}$${\displaystyle {\vec {s}}_{t}}$${\estilo de visualización t}$

${\displaystyle {\frac {\partial V_{\text{bias}}({\vec {s}}\,)}{\partial t}}=\omega \,\delta (|{\vec {s}}-{\vec {s}}_{t}|)}$.

En un tiempo de simulación infinitamente largo , el potencial de polarización acumulado converge a energía libre con signo opuesto (y constante irrelevante ): ${\displaystyle t_{\text{sim}}}$${\displaystyle C}$

${\displaystyle V_{\text{bias}}({\vec {s}}\,)=\!\!\int _{0}^{t_{\text{sim}}}\!\!\!\omega \,\delta (|{\vec {s}}-{\vec {s}}_{t}|)\;dt\quad \Rightarrow \quad F({\vec {s}}\,)=-\!\!\!\!\lim _{t_{\text{sim}}\to \infty }\!\!V_{\text{bias}}({\vec {s}}\,)+C}$

Para lograr una implementación computacionalmente eficiente, el proceso de actualización se discretiza en intervalos de tiempo ( denota la función base ) y la función se reemplaza con una función kernel positiva localizada . El potencial de sesgo se convierte en una suma de las funciones kernel centradas en los valores de las variables colectivas instantáneas en el tiempo : ${\displaystyle \tau }$${\displaystyle \lfloor \;\rfloor }$${\displaystyle \delta }$ ${\displaystyle K}$${\displaystyle {\vec {s}}_{j}}$${\displaystyle \tau j}$

${\displaystyle V_{\text{bias}}({\vec {s}}\,)\approx \tau \!\!\!\sum _{j=0}^{\left\lfloor {\frac {t_{\text{sim}}}{\tau }}\right\rfloor }\!\!\omega \,K(|{\vec {s}}-{\vec {s}}_{j}|)}$.

Normalmente, el núcleo es una función gaussiana multidimensional , cuya matriz de covarianza solo tiene elementos diagonales distintos de cero:

${\displaystyle V_{\text{bias}}({\vec {s}}\,)\approx \tau \!\!\!\sum _{j=0}^{\left\lfloor {\frac {t_{\text{sim}}}{\tau }}\right\rfloor }\!\!\omega \exp \!\!\left(\!-{\frac {1}{2}}\left|{\frac {{\vec {s}}-{\vec {s}}_{j}}{\vec {\sigma }}}\right|^{2}\right)}$.

Los parámetros , , y se determinan a priori y se mantienen constantes durante la simulación. ${\displaystyle \tau }$${\displaystyle \omega }$${\displaystyle {\vec {\sigma }}}$

A continuación se muestra un pseudocódigo de MTD basado en dinámica molecular (MD), donde y son las posiciones y velocidades del sistema de partículas, respectivamente. El sesgo se actualiza cada paso de MD y su contribución a las fuerzas del sistema es . ${\displaystyle \{{\vec {r}}\}}$${\displaystyle \{{\vec {v}}\}}$${\displaystyle N}$${\displaystyle V_{\text{bias}}}$${\displaystyle n=\tau /\Delta t}$${\displaystyle \{{\vec {F}}\,\}}$${\displaystyle \{{\vec {F}}_{\text{bias}}\}}$

establecer inicial y establecer${\displaystyle \{{\vec {r}}\}}$${\displaystyle \{{\vec {v}}\}}$ 
 ${\displaystyle V_{\text{bias}}({\vec {s}}\,):=0}$

Cada paso de MD:
     calcula los valores de CV:
        ${\displaystyle {\vec {s}}_{t}:={\vec {s}}(\{{\vec {r}}\})}$
    
    Cada  paso del MD:
         actualizar el potencial de sesgo:
            ${\displaystyle n}$${\displaystyle V_{\text{bias}}({\vec {s}}\,):=V_{\text{bias}}({\vec {s}}\,)+\tau \omega \exp \!\!\left(\!-{\frac {1}{2}}\left|{\frac {{\vec {s}}-{\vec {s}}_{t}}{\vec {\sigma }}}\right|^{2}\right)}$
    
    Calcular fuerzas atómicas:
        ${\displaystyle {\vec {F}}_{i}:=-{\frac {\partial V(\{{\vec {r}}\,\})}{\partial {\vec {r}}_{i}}}\overbrace {\left.-{\frac {\partial V_{\text{bias}}({\vec {s}}\,)}{\partial {\vec {s}}}}\right|_{{\vec {s}}_{t}}\!\!\!{\frac {\partial {\vec {s}}(\{{\vec {r}}\,\})}{\partial {\vec {r}}_{i}}}} ^{{\vec {F}}_{{\text{bias}},i}}}$
    
    propagar  y por${\displaystyle \{{\vec {r}}\}}$${\displaystyle \{{\vec {v}}\}}$${\displaystyle \Delta t}$

El tamaño finito del núcleo hace que el potencial de polarización fluctúe alrededor de un valor medio. Se puede obtener una energía libre convergente promediando el potencial de polarización. El promedio comienza desde , cuando el movimiento a lo largo de la variable colectiva se vuelve difusivo: ${\displaystyle t_{\text{diff}}}$

${\displaystyle {\bar {F}}({\vec {s}})=-{\frac {1}{t_{\text{sim}}-t_{\text{diff}}}}\int _{t_{\text{diff}}}^{t_{\text{sim}}}\!\!\!\!\!V_{\text{bias}}({\vec {s}},t)\,dt+C}$

La metadinámica se ha utilizado para estudiar:

• plegamiento de proteínas ^[22]
• reacciones químicas ^[41]
• acoplamiento molecular ^{[42] [43]}
• transiciones de fase . ^[44]
• encapsulación de ADN en nanotubos de carbono de pared simple hidrófobos ^[45] e hidrófilos ^{[46] .}

PLUMED ^[47] es una biblioteca de código abierto que implementa muchos algoritmos MTD y variables colectivas . Tiene un diseño orientado a objetos flexible ^{[48] [49]} y se puede interconectar con varios programas MD ( AMBER , GROMACS , LAMMPS , NAMD , Quantum ESPRESSO , DL_POLY_4, CP2K y OpenMM). ^{[50] [51]}

Existen otras implementaciones de MTD en el Módulo de Variables Colectivas ^[52] (para LAMMPS , NAMD y GROMACS ), ORAC , CP2K , ^[53] EDM, ^[54] y Desmond .

• Introducción a la metadinámica
• CON PENACHO
• Sitio web del módulo Colvars (NAMD, LAMMPS, Gromacs)
• Película visual de metadinámica
• Muestreo de probabilidad mejorada sobre la marcha (OPES)

• Elevación local
• Templado paralelo
• Muestreo de paraguas