Blog

Dec 11, 2023

Máquina

Nature Computational Science volumen 3, páginas 334–345 (2023)Cite este artículo 8339 Accesos 8 Citas 44 Detalles de métricas altmétricas Autoorganización molecular impulsada por muchos cuerpos concertados

Nature Computational Science volumen 3, páginas 334–345 (2023)Cite este artículo

8339 Accesos

8 citas

44 altmétrica

Detalles de métricas

La autoorganización molecular impulsada por interacciones concertadas de muchos cuerpos produce las estructuras ordenadas que definen tanto la materia viva como la inanimada. Aquí presentamos un algoritmo de muestreo de ruta autónomo que integra el aprendizaje profundo y la teoría de la ruta de transición para descubrir el mecanismo de los fenómenos de autoorganización molecular. El algoritmo utiliza el resultado de trayectorias recién iniciadas para construir, validar y, si es necesario, actualizar modelos mecanicistas cuantitativos. Al cerrar el ciclo de aprendizaje, los modelos guían el muestreo para mejorar el muestreo de eventos de ensamblaje poco comunes. La regresión simbólica condensa el mecanismo aprendido en una forma interpretable por humanos en términos de observables físicos relevantes. Aplicado a la asociación de iones en solución, la formación de cristales de hidrato de gas, el plegamiento de polímeros y el ensamblaje de membrana-proteína, capturamos los movimientos de solventes de muchos cuerpos que gobiernan el proceso de ensamblaje, identificamos las variables de la teoría de la nucleación clásica, descubrimos el mecanismo de plegamiento en diferentes niveles de resolución y revelar vías de montaje en competencia. Las descripciones mecanicistas son transferibles a través de estados termodinámicos y espacio químico.

Comprender cómo cooperan las interacciones genéricas pero sutilmente orquestadas en la formación de estructuras complejas es la clave para dirigir el autoensamblaje molecular1,2. Como experimentos informáticos, las simulaciones de dinámica molecular (MD) nos prometen vistas atómicamente detalladas e imparciales de los procesos de autoorganización3. Sin embargo, la mayoría de los procesos de autoorganización colectiva son eventos raros que ocurren en escalas de tiempo de muchos órdenes de magnitud más largas que los rápidos movimientos moleculares que limitan el paso de integración de MD. El sistema pasa la mayor parte del tiempo en estados metaestables, y las transiciones estocásticas rápidas e infrecuentes entre estados rara vez se resuelven en simulaciones MD imparciales, si es que se resuelven. Estas rutas de transición (TP) son segmentos de trayectoria muy especiales que capturan el proceso de reorganización. Aprender mecanismos moleculares a partir de simulaciones requiere que la potencia computacional se centre en tomar muestras de TP4 y destilar modelos cuantitativos a partir de ellos5. Debido a la alta dimensionalidad del espacio de configuración, tanto el muestreo como la extracción de información son extremadamente desafiantes en la práctica. Nuestro algoritmo aborda ambos desafíos a la vez. Construye de forma autónoma y simultánea modelos mecanicistas cuantitativos de eventos moleculares complejos, valida los modelos sobre la marcha y los utiliza para acelerar el muestreo en órdenes de magnitud en comparación con la MD normal.

La mecánica estadística proporciona un marco general para obtener modelos mecanicistas de baja dimensión de eventos de autoorganización. En este artículo, nos centramos en sistemas que se reorganizan entre dos estados A y B (ensamblados o desensamblados, respectivamente), pero generalizar a un número arbitrario de estados es sencillo. Cada TP que conecta los dos estados contiene una secuencia de instantáneas que capturan el sistema durante su reorganización. En consecuencia, el conjunto de trayectorias de transición (TPE) es el mecanismo de mayor resolución. Como la transición es efectivamente estocástica, cuantificar su mecanismo requiere un marco probabilístico. Definimos el comprometidor pS(x) como la probabilidad de que una trayectoria entre en el estado S primero, con S = A o B, respectivamente, donde x es un vector de características que representan el punto inicial X en el espacio de configuración, y pA(x) + pB(x) = 1 para dinámica ergódica. El comprometidor pB informa sobre el progreso de la reacción A → B y predice el destino de la trayectoria en forma markoviana6,7, lo que la convierte en la coordenada de reacción ideal8,9. En el juego de ajedrez, uno puede pensar en el comprometidor como la probabilidad de, digamos, que las negras ganen para determinadas posiciones iniciales en el tablero en partidas repetidas10. Los requisitos mínimos para aplicaciones más allá de las simulaciones moleculares son (1) que exista una cantidad similar a un comprometidor y (2) que la dinámica del sistema pueda muestrearse repetidamente, al menos en la dirección directa. Por lo tanto, la probabilidad de diferentes eventos posibles (A, B,…) debería codificarse al menos en parte (y por lo tanto aprenderse en términos de) el estado instantáneo X del sistema y la dinámica del sistema debería ser susceptible de muestreo repetido, preferiblemente mediante una eficiente simulación por ordenador. Sin embargo, si uno puede preparar repetidamente un sistema real con un control satisfactorio sobre las condiciones iniciales, puede aprender a predecir el destino probable de este sistema dadas las condiciones iniciales observadas y controladas utilizando nuestro marco.

El muestreo de TP para eventos raros y el aprendizaje de la función de compromiso asociada pB (x) son dos desafíos sobresalientes e intrínsecamente conectados. Dado que los TP son extremadamente raros en un espacio de muy altas dimensiones, una búsqueda desinformada es inútil. Sin embargo, los TP convergen cerca de los estados de transición7, donde la trayectoria está menos comprometida con pA(x) = pB(x) = 1/2. Para la dinámica markoviana y reversible en el tiempo P(TP|x), la probabilidad de que una trayectoria que pasa por x sea un TP satisface P(TP|x) = 2pB(x)(1 − pB(x)), es decir , el comprometidor determina la probabilidad de muestrear un TP11. Por lo tanto, los desafíos de la extracción de información y el muestreo están entrelazados.

Para abordar estos desafíos duales, diseñamos un algoritmo iterativo que se basa en la teoría de la ruta de transición9 y el muestreo de la ruta de transición (TPS)12 con el espíritu del aprendizaje por refuerzo. El algoritmo aprende al autor de eventos raros en sistemas complejos de muchos cuerpos mediante la ejecución repetida de experimentos virtuales y utiliza el conocimiento adquirido para mejorar el muestreo de TP (Fig. 1a). En cada experimento, el algoritmo selecciona un punto X desde el cual disparar trayectorias (propagadas de acuerdo con la dinámica imparcial del modelo físico) para generar TP. Para garantizar un equilibrio detallado de TPS, el algoritmo selecciona estructuras de la ruta de transición actual para iniciar movimientos de disparo con velocidades iniciales rediseñadas de Maxwell-Boltzmann. Después de disparos repetidos desde diferentes puntos X, el algoritmo compara el número de TP generados con el número esperado en función de su conocimiento del autor del compromiso en ese punto. Sólo si la predicción es deficiente, el algoritmo vuelve a entrenar el modelo según el resultado de todos los experimentos virtuales, lo que evita el sobreajuste. A medida que aumenta el poder predictivo del modelo mecanicista, el algoritmo se vuelve más eficiente en el muestreo de TP al elegir puntos iniciales cerca de los estados de transición, es decir, de acuerdo con P (TP | x).

a, Mecanismo de aprendizaje mediante muestreo de ruta. El método itera entre el muestreo de rutas de transición de una configuración x entre los estados metaestables A y B (izquierda) y el aprendizaje del confirmador pB(x) (derecha). Una función de red neuronal de características moleculares (x1 a x4) modela al comprometidor. No se muestra el predictor de registro que forma la última capa. En la convergencia, la regresión simbólica destila una expresión interpretable que cuantifica el mecanismo molecular en términos de características seleccionadas (x1, x2) y constantes numéricas (a, c) conectadas por operaciones matemáticas (aquí: +, −, ×, exp). b, Instantáneas a lo largo de un TP que muestran la formación de un par de iones LiCl (de derecha a izquierda) en una simulación MD atomística. El agua se muestra como palos, el Li+ como una esfera pequeña y el Cl− como una esfera grande. Los átomos se colorean según su contribución al progreso de la reacción desde bajo (azul) hasta alto (rojo), cuantificado por su contribución al gradiente de la coordenada de reacción q(x|w). c, Autoconsistencia. Recuentos del número de eventos de transición generados (línea azul) y esperados (línea discontinua naranja). La línea verde muestra la diferencia acumulada entre los recuentos observados y esperados. El recuadro muestra un acercamiento a las primeras 1000 iteraciones. d, Validación del comprometidor aprendido. Correlación cruzada entre el confirmador predicho por la red entrenada y el confirmador obtenido mediante muestreo repetido de configuraciones moleculares en las que no se entrenó el modelo de comprometidor. El promedio de los confirmadores muestreados (línea azul) y su sd (sombreado en naranja) se calcularon en los contenedores del confirmador aprendido indicados por los pasos verticales. Como referencia, la línea roja indica la identidad. e, Transferibilidad del comprometidor aprendido. Representación del aprendizaje por transferencia y correlaciones cruzadas entre los comprometidores muestreados para el emparejamiento de iones NaCl y NaI y predicciones del comprometidor a partir de un modelo entrenado con datos para LiCl y ajustado mediante aprendizaje por transferencia utilizando solo 1000 resultados de disparo adicionales cada uno. Los colores y sd (indicados con un sombreado naranja) son como en d.

Datos fuente

El algoritmo aprende de sus repetidos intentos utilizando el aprendizaje profundo5,13,14 de forma autoconsistente. Modelamos el comprometidor pB(x) = 1/(1 + e−q(x|w)) con una red neuronal15 q(x|w) de pesos w. Con esta elección, q es una función invertible de pB y podemos ver ambas funciones como la coordenada de reacción. En cada intento de generar un TP, el algoritmo propaga dos trayectorias, una hacia adelante y otra hacia atrás en el tiempo, ejecutando simulaciones MD4. En caso de éxito, una trayectoria entra primero en el estado A y la otra en B, formando conjuntamente un nuevo TP. En este proceso de lanzamiento de moneda de Bernoulli, el algoritmo aprende tanto de los éxitos como de los fracasos. La probabilidad logarítmica negativa5 para k intentos define la función de pérdida \(l({{{\bf{w}}}}| {{{\bf{\uptheta }}}})=\mathop{\sum }\nolimits_ {i = 1}^{k}\log (1+{{\mathrm{e}}}^{{s}_{i}q({{{{\bf{x}}}}}_{i) }| {{{\bf{w}}}})})\), donde si = 1 si la trayectoria i iniciada desde Xi entra primero en A, y si = −1 si entra primero en B. El conjunto de entrenamiento θ contiene los k vectores de características xi asociados con los puntos de disparo Xi y los resultados si. Al entrenar la red q(x|w) para minimizar la pérdida l(w|θ), obtenemos una estimación de máxima verosimilitud del comprometidor5 que es diferenciable y permite métodos de análisis sofisticados16.

Luego utilizamos la regresión simbólica17 para condensar el mecanismo molecular en una forma interpretable por humanos (Fig. 1a) y obtener información física. Primero, un análisis de atribución de la red entrenada identifica un pequeño subconjunto z de las coordenadas de entrada x que dominan el resultado de la predicción de la red. Luego, la regresión simbólica destila expresiones matemáticas explícitas qsr(z|wsr) mediante el uso de un algoritmo genético que busca espacios de funciones y parámetros para minimizar la pérdida l(wsr|θ) en el conjunto de entrenamiento θ independiente del entrenamiento de la red neuronal anterior, donde el subíndice 'sr' indica regresión simbólica. Las expresiones analíticas resultantes nos proporcionan una lista de hipótesis para modelos cuantitativos de la física que gobierna el proceso de ensamblaje molecular. Para un examen más detallado, estas hipótesis se clasifican según una combinación de probabilidad estadística (es decir, qué tan bien explican todos los datos disponibles) y complejidad matemática.

La formación de pares de iones en agua es un proceso de ensamblaje paradigmático controlado por interacciones de muchos cuerpos en el medio disolvente circundante. Aunque la MD puede simular eficientemente el proceso, la reorganización colectiva de las moléculas de agua desafía la formulación de modelos mecanicistas cuantitativos hasta el día de hoy (Fig. 1b).

El algoritmo aprendió rápidamente cómo muestrear de manera óptima la formación de pares de iones. Para los pares de iones de litio (Li+) y cloruro (Cl−) en agua (Fig. 1b,c), la red utilizó la distancia interiónica rLiCl y 220 características moleculares x1,…, x220 que describen la disposición angular de los átomos de oxígeno e hidrógeno del agua. a una distancia específica de cada ion19. Estas coordenadas proporcionan una representación general del sistema que es invariante con respecto a las simetrías físicas y el intercambio de etiquetas atómicas. Después de las primeras 500 iteraciones, los números de TP previstos y observados coinciden (Fig. 1c). El muestreo es aproximadamente diez veces más rápido que el TPS convencional (Datos ampliados, figura 1). Observamos que esta aceleración se logra enteramente mejorando la eficiencia del muestreo de nuevas rutas de transición y sin sesgos en la dinámica subyacente. Además, validamos la función de compromiso aprendida comparando sus predicciones con simulaciones independientes. A partir de 1000 configuraciones no utilizadas en el entrenamiento, iniciamos 500 simulaciones independientes cada una y estimamos el pB del comprometidor muestreado como la fracción de trayectorias que ingresan por primera vez al estado libre. Los comprometidos previstos y muestreados están en acuerdo cuantitativo (Fig. 1d).

Contrariamente a una preocupación común por los modelos de aprendizaje automático, el mecanismo aprendido es general y, con ajustes mínimos, describe el ensamblaje de especies iónicas químicamente distintas. Realizamos aprendizaje por transferencia en cinco sistemas adicionales al permitir modificaciones solo en la última capa lineal de la red entrenada que contiene una sola neurona (Fig. 1e y Datos extendidos Fig. 2). Una cantidad muy limitada de nuevas transiciones simuladas es suficiente para ajustar la red que contiene el comprometidor LiCl para predecir correctamente el comprometidor para LiI, NaCl, NaI, CsCl y CsI.

También construimos modelos de múltiples iones que se extienden por el espacio químico. Como informadores sobre el tamaño y la energía de los iones, incluimos los parámetros tamaño de partícula σ y energía de dispersión ϵ del potencial de Lennard-Jones en los vectores de características x. Descubrimos que los modelos entrenados con las estadísticas de TP combinadas para diferentes combinaciones de pares iónicos pueden interpolar y extrapolar en el espacio químico, haciendo predicciones razonables para el mecanismo de asociación de especies iónicas en las que no se han entrenado (Datos ampliados, figura 3).

Los reordenamientos de disolventes desempeñan un papel fundamental en la determinación del ensamblaje de iones. El análisis de atribución para un modelo entrenado en el ensamblaje de LiCl, LiI, NaCl, NaI, CsCl y CsI identificó simultáneamente el rión de distancia interiónico y los parámetros de Lennard-Jones como cruciales (Fig. 2a). Además, las funciones de simetría que describen la geometría de las moléculas de agua alrededor del catión controlan el mecanismo de ensamblaje. Como la más importante de las 176 características moleculares utilizadas para describir el disolvente, x7 cuantifica la anisotropía del oxígeno a una distancia radial de 0,1 nm de los cationes (Fig. 2b). Para un ensamblaje exitoso de pares iónicos, estas moléculas de agua de la capa interna necesitan abrir espacio para el anión entrante. La importancia del reordenamiento del agua de la capa interna es consistente con un análisis visual que resalta los átomos en un TP de acuerdo con su contribución al gradiente del comprometidor (Fig. 1b).

a, Relevancia de entrada para las 179 coordenadas de entrada utilizadas para el aprendizaje profundo. Los primeros 176 describen la geometría de las moléculas de agua alrededor de cationes y aniones. Los restantes son el rión de distancia interiónico y los parámetros de Lennard-Jones, siendo σ el tamaño del ión. b, Representación esquemática de la reorientación de disolventes más importante. La función de simetría x7 informa la geometría de los átomos de oxígeno del agua (O, en azul) a 0,1 nm alrededor del catión (en rosa) (consulte el cuadro para ver la definición de x7, donde rij y rik son las distancias entre el catión central i y el oxígeno átomos j y k, y ϑijk es el ángulo formado por el catión central y dos átomos de oxígeno). c, Gráfico de Pareto de todos los modelos destilados por regresión simbólica. Cada punto corresponde a un modelo alternativo qsr(z|wsr), coloreado según el número de coordenadas de entrada (Nin) que utiliza. La cruz roja identifica el modelo óptimo en la rodilla del frente de Pareto. d, Modelo de iones múltiples a partir de regresión simbólica que describe el mecanismo de ensamblaje de LiCl, LiI, NaCl, NaI, CsCl y CsI en agua. El modelo, \(q\left({r}_{{\mathrm{ion}}},\sigma ,{x}_{7},{x}_{15};{\sigma }_{w} \right)\), es una función de la distancia interiónica (rión y tamaño de ion σ en unidades del parámetro de tamaño del agua σw = 0,315 nm) y la geometría del agua alrededor de los cationes (x7 y x15). e, Validación del modelo de ensamblaje de iones múltiples mediante correlación cruzada entre participantes muestreados no entrenados y la predicción para cada especie iónica por separado, aquí se muestra para LiCl y CsCl (consulte Datos ampliados en la Fig. 4a para todas las especies restantes). El promedio de los confirmadores muestreados (línea azul) y su sd (sombreado en naranja) se calcularon en los contenedores del confirmador aprendido indicados por los pasos verticales. Como referencia, la línea roja indica la identidad.

Datos fuente

La regresión simbólica proporciona un modelo multiion cuantitativo e interpretable del mecanismo de ensamblaje. En regresiones simbólicas independientes, variamos el número de entradas y la penalización por complejidad. Luego seleccionamos modelos en un diagrama de Pareto (Fig. 2c). Los modelos en el punto culminante del frente de Pareto ofrecen buenas compensaciones entre la calidad del modelo, medida por la pérdida, y la complejidad del modelo, medida por el número de operaciones matemáticas.

El modelo destilado de iones múltiples es interpretable y proporciona información física sobre el ensamblaje de iones monovalentes en agua (Fig. 2d). En el término principal de q, un parámetro de tamaño de ion escalado σ se resta de la distancia interiónica, de acuerdo con la intuición física. En el segundo término, el tamaño del ion modula de forma no lineal el descriptor x7 de la geometría del agua cerca de los cationes (Fig. 2b). En el último término, x15 informa sobre la solvatación más lejana sin modular por la identidad del ion. A pesar de su simplicidad, el modelo reducido es preciso para todas las especies de iones monovalentes consideradas aquí (Fig. 2e y Datos ampliados, Fig. 4a). Un modelo de regresión simbólica que se centra en el ensamblaje de LiCl solo muestra que podemos intercambiar menos generalidad por una mayor precisión (Datos ampliados, figuras 4b-d).

A baja temperatura y alta presión, una mezcla líquida de agua y metano se organiza en un hidrato de gas, un sólido parecido al hielo20. En esta transición de fase, las moléculas de agua se ensamblan en una intrincada red cristalina con jaulas regularmente espaciadas y llenas de metano (Fig. 3a). A pesar de la relevancia comercial en el procesamiento del gas natural, el mecanismo de formación de hidratos de gas aún no se comprende completamente, complicado por el carácter de muchos cuerpos del proceso de nucleación y la competencia entre diferentes formas cristalinas20. Estudiar el mecanismo de nucleación es un desafío para los experimentos y, debido a la extrema rareza de los eventos, es imposible en equilibrio MD.

a, Configuraciones moleculares que ilustran el proceso de nucleación extraídas de una simulación MD atomística en solvente explícito. El núcleo se forma en agua, crece y conduce al cristal de clatrato. Las jaulas (líneas) de agua 51262 (azul) y 512 (roja) contienen moléculas de metano (esferas) de colores correspondientes. Las moléculas de metano cerca del núcleo sólido en crecimiento se muestran como esferas verdes y el agua como palos grises. Para mayor claridad, no se muestra el agua a granel. b, Validación del comprometidor aprendido. Correlación cruzada entre el confirmador predicho por la red entrenada y el confirmador obtenido mediante muestreo repetido de configuraciones moleculares en las que el algoritmo no se entrenó (línea gris: identidad). c, Análisis de importancia de los insumos. Las tres coordenadas de entrada más importantes están anotadas como temperatura T, el número de aguas superficiales nw y el número de 51262 cristales nc. d, El modelo mecanicista cuantitativo basado en datos destilado por regresión simbólica revela un cambio en el mecanismo de nucleación. En la ecuación, nw,0 y T0 son el número de referencia de moléculas de agua superficial y la temperatura de referencia, respectivamente, y α, β, γ y δ son constantes numéricas. Superficies iso-committor analíticas para nw,0 = 2, T0 = 270 K, α = 0,0502, β = 3,17, γ = 0,109 K−1, δ = 0,0149 (de izquierda a derecha: amarillo, pB = 1/(1 + e −4); azul, 1/2; verde, 1/(1 + e4)). Los recuadros estructurales ilustran los dos mecanismos competitivos a baja y alta temperatura.

Datos fuente

A las pocas horas de tiempo de cálculo en una sola unidad de procesamiento de gráficos (GPU), el algoritmo extrajo el mecanismo de nucleación de 2225 TP que muestran la formación de clatratos de metano, correspondientes a un total de 445,3 μs de dinámica simulada. Las trayectorias fueron producidas mediante extensas simulaciones de TPS a cuatro temperaturas diferentes y proporcionaron un conjunto de entrenamiento preexistente para nuestro algoritmo21. Describimos configuraciones moleculares utilizando 22 características comúnmente utilizadas para describir procesos de nucleación (Tabla complementaria 1). Consideramos la temperatura T a la que se generó un TP como una característica adicional y entrenamos el modelo de compromiso en las trayectorias acumulativas. Demostramos que el comprometidor aprendido en función de la temperatura es preciso al validar sus predicciones para 160 configuraciones independientes (Fig. 3b). Los modelos generativos construyeron recientemente funciones de distribución a temperaturas no muestreadas22. Al omitir datos en T = 280 K o 285 K en el entrenamiento, mostramos que el comprometidor aprendido interpola y extrapola satisfactoriamente a estados termodinámicos no muestreados (Datos extendidos, Fig. 5).

La temperatura T es el factor más crítico para el resultado de una trayectoria de simulación, seguida por el número nw de moléculas de agua superficial y el número nc de 51262 jaulas, definido por la presencia de 12 pentágonos (512) y dos hexágonos (62) (Fig. .3c). Las tres variables juegan un papel esencial en la teoría clásica de la nucleación homogénea21. La energía libre de activación ΔG para la nucleación está determinada por el tamaño del núcleo en crecimiento, parametrizado por la cantidad de agua superficial y, en el caso de una estructura cristalina, el número de 51262 jaulas. La temperatura determina, a través del grado de sobresaturación, el tamaño del núcleo crítico, la altura de la barrera de energía libre de nucleación y la velocidad.

La regresión simbólica destiló una expresión matemática que revela un cambio dependiente de la temperatura en el mecanismo de nucleación. El mecanismo se cuantifica mediante q (nw, nc, T) (Fig. 3d y Tabla complementaria 1). A bajas temperaturas, el tamaño del núcleo por sí solo decide el crecimiento. A temperaturas más altas, el número de 51262 jaulas de agua gana en importancia, como lo indican las superficies curvas de iso-committor (Fig. 3d). Este modelo mecanicista, generado basándose en datos, revela el cambio del crecimiento amorfo a bajas temperaturas al crecimiento cristalino a temperaturas más altas21,23.

Las proteínas, los ácidos nucleicos y los polímeros pueden autoorganizarse espontáneamente plegándose en estructuras ordenadas. Aplicado a la transición de bobina a cristal de un homopolímero24,25, el algoritmo identificó fácilmente el mecanismo previamente difícil de alcanzar en diferentes niveles de resolución (Datos ampliados, figura 6). A baja resolución, utilizamos un conjunto selecto de 36 características físicas promediadas sobre el polímero. El análisis de atribución seguido de una regresión simbólica representó al comprometidor como una función no lineal de orden orientacional Q6 y energía potencial U sola, lo que resultó altamente predictivo (Fig. 4 y Datos ampliados Fig. 6a). En alta resolución, el aprendizaje profundo produjo una función de compromiso de calidad comparable en un espacio de 384 descriptores generales que representan el entorno local de cada perla de polímero (Datos extendidos, figura 6b) en términos del número de vecinos, el parámetro de orden de enlace local q6 y los coeficientes de conexión cij que miden la correlación entre los entornos locales de las cuentas i y j. De este modo, el algoritmo aprendió representaciones precisas del comprometidor en términos de muchas características generales y pocas características específicas del sistema, y ​​destiló estas últimas en una función compacta y físicamente reveladora de orden de orientación y energía.

a, Representación del mecanismo aprendido. El mapa de calor (barra de color) representa un modelo explícito reducido del comprometidor pB = pF al estado plegado como lo reproduce la expresión \({q}_{{\mathrm{B}}}\left(U,{Q} _{6}\right)=\alpha (U-{U}_{0})+\beta \log \left({Q}_{6}-{Q}_{6,0}\right)+ \gamma\), donde U es la energía potencial total del polímero, Q6 cuantifica su cristalinidad y las constantes numéricas son α = −7.144, β = 3.269, γ = 11.942, U0 = −2.351 y Q6,0 = 0.035. Recuadros: configuraciones moleculares del polímero en pB = 0, 0,5 y 1. Las perlas de polímero se colorean según su valor de q6, desde el blanco (valores bajos) hasta el naranja oscuro (valores altos). b, Validación del comprometidor aprendido. Correlación cruzada entre el confirmador predicho por la red entrenada y el confirmador obtenido mediante muestreo repetido de configuraciones moleculares en las que el algoritmo no se entrenó. El promedio de los confirmadores muestreados (línea azul) y su sd (sombreado en naranja) se calcularon en los contenedores del confirmador aprendido indicados por los pasos verticales. Como referencia, la línea roja indica la identidad.

Datos fuente

Los complejos membrana-proteína desempeñan un papel fundamental en la organización de las células vivas. Aquí investigamos el ensamblaje del homodímero transmembrana del sensor de saturación Mga2 en una bicapa lipídica en la representación cuasi atomística de Martini (Fig. 5a) 26. En extensas simulaciones de MD en equilibrio, se ha observado la asociación espontánea de dos hélices transmembrana de Mga2, pero no se produjo ninguna disociación en aproximadamente 3,6 ms (equivalente a más de 6 meses de cálculos)26.

a, Instantáneas durante un evento de dimerización de Mga2 (de derecha a izquierda). Las hélices transmembrana se muestran como superficies naranjas, las moléculas de lípidos en gris y el agua en cian. b, Autoconsistencia. Recuentos acumulativos del número de transiciones generadas (línea azul) y esperadas (línea discontinua naranja). La curva verde muestra la diferencia acumulada entre los recuentos observados y esperados. c, Validación del comprometidor aprendido. Correlación cruzada entre el confirmador predicho por la red entrenada y el confirmador obtenido mediante muestreo repetido de configuraciones moleculares en las que no se entrenó el modelo de comprometidor. El promedio de los confirmadores muestreados (línea azul) y su sd (sombreado en naranja) se calcularon en los contenedores del confirmador aprendido indicados por los pasos verticales. Como referencia, la línea roja indica la identidad. d, Representación esquemática de las dos coordenadas más relevantes, los contactos interhelicoidales en las posiciones 9 y 22. e, Agrupación jerárquica de todos los TP. Dendrograma en función de la similitud de TP (deformación dinámica del tiempo, consulte 'Ensamblaje de dímero transmembrana de Mga2 en membrana lipídica' en Métodos) calculado en el plano definido por los contactos 9 y 22 (dos grupos principales: azul, naranja). f, g, Densidad de ruta (sombreado gris) para los dos grupos principales en e, calculada en el plano definido por los contactos 9 y 22. Para cada grupo, se muestra un TP representativo desde libre (amarillo) hasta enlazado (rojo). Las isolíneas del comprometidor, como lo predice la regresión simbólica \({q}_{{\mathrm{B}}}({x}_{9},{x}_{22})=-\exp ({ x}_{9}^{2})\log ({x}_{9}-\frac{{x}_{9}}{\log ({x}_{22})})\), se muestran como líneas continuas etiquetadas.

Datos fuente

El muestreo de rutas es naturalmente paralelizable, lo que nos permitió muestrear casi 4000 eventos de disociación en 20 días en una supercomputadora paralela (Fig. 5b). La integración temporal de las trayectorias MD conlleva el mayor coste computacional y sólo es paralelizable hasta un grado limitado. Sin embargo, una sola instancia del algoritmo puede orquestar simultáneamente experimentos virtuales en un número arbitrario de copias del modelo físico (guiando procesos de muestreo paralelos de la cadena de Markov Monte Carlo (MC)) y aprender de todos ellos entrenando sobre los resultados acumulativos. .

Caracterizamos configuraciones moleculares utilizando contactos entre los residuos correspondientes a lo largo de las dos hélices e incluimos, como referencia, una serie de características personalizadas que describen la organización de los lípidos alrededor de las proteínas27 (Datos ampliados, figura 7 y tabla complementaria 2). Validamos el modelo con datos del confirmador para 548 configuraciones moleculares que no se utilizan en el entrenamiento y encontramos que las predicciones son precisas en toda la región de transición entre estados unidos y libres (Fig. 5c).

En una notable reducción de la dimensionalidad, la regresión simbólica logró una representación precisa del comprometidor aprendido como una función simple de solo dos contactos de aminoácidos (Fig. 5d y Datos ampliados, Fig. 8). La regresión simbólica nos proporciona una lista de hipótesis para modelos cuantitativos de la física que gobierna el proceso de ensamblaje molecular (Tabla complementaria 3). Estas hipótesis se clasifican según una combinación de probabilidad estadística (es decir, qué tan bien representan todos los datos disponibles) y su complejidad matemática. Entre las expresiones en el punto culminante del diagrama de Pareto, es decir, con poder predictivo y complejidad comparables, elegimos la que ofrece una interpretación clara en términos de dos mecanismos de ensamblaje en competencia descritos por la formación de contactos que comienzan en un extremo de la hélice u otro. .

Proyectamos todos los TP muestreados en el plano definido por estos dos contactos, calculamos las distancias entre ellos y realizamos una agrupación jerárquica de trayectorias (Fig. 5e). Los TP se organizan en dos grupos principales que revelan dos vías de ensamblaje en competencia con el contacto inicial de la hélice en la parte superior (Fig. 5f) o en la parte inferior (Fig. 5g). Inesperadamente27, la geometría hélice-dímero por sí sola predice el progreso del ensamblaje, lo que implica que el 'solvente' lipídico está implícitamente codificado en los contactos interhelicoidales por pares, a diferencia del solvente agua en la formación de pares iónicos18. Al igual que en el plegamiento de polímeros y la unión de iones, un espacio suficientemente grande de características geométricas generales resultó suficiente para la construcción de confirmadores totalmente predictivos mediante aprendizaje profundo. Este hallazgo es consistente con la teoría de la incrustación28 e implica que el uso de un número pequeño pero suficiente de características generales es tan efectivo como variables colectivas basadas en la intuición física y química.

El muestreo de trayectoria guiado por máquina es general y se puede adaptar inmediatamente para muestrear la dinámica de muchos cuerpos con una noción de "destino probable" similar a la del autor. Este concepto fundamental de la mecánica estadística se extiende desde el juego de ajedrez10 sobre el plegamiento de proteínas3,7 hasta el modelado climático29. El motor de simulación (MD en nuestro caso) se trata como una caja negra y puede ser reemplazado por otros procesos dinámicos, reversibles o no. Tanto el modelo estadístico que define la función de pérdida como la tecnología de aprendizaje automático pueden adaptarse a problemas específicos. Los modelos más sofisticados podrán aprender más con menos datos o incorporar restricciones experimentales. Esquemas de regresión más simples5 pueden reemplazar a las redes neuronales15 cuando el costo de las trayectorias de muestreo limita severamente el volumen de datos de entrenamiento.

Definir los límites de los estados metaestables, como lo requiere nuestro método, no puede ser trivial. Para refinar las definiciones de estado, el marco del comprometidor se puede utilizar de forma iterativa, comenzando desde uno muy conservador, obteniendo una estimación aproximada del comprometidor y luego extendiendo los estados a valores del comprometidor cercanos a 0 y 1. El problema es la presencia de estados metaestables no asignados que dan como resultado largos caminos de transición. En este caso, los enfoques de aprendizaje automático destinados al descubrimiento de estados pueden proporcionar una solución.

Como nuestro método no utiliza fuerzas no físicas para acelerar el muestreo, su aplicabilidad se limita a procesos para los cuales las rutas de transición son lo suficientemente cortas como para que el motor de simulación subyacente pueda muestrear un número razonable de ellas. Esta limitación no es dramática ya que la duración de los TP depende sólo débilmente de la altura de la barrera. Si es necesario, nuestro método también se puede utilizar con superficies potenciales sesgadas para acelerar la dinámica.

La validación experimental del mecanismo de ensamblaje puede provenir de la observación directa, por ejemplo, mediante experimentos con una sola molécula, o por perturbación, por ejemplo, cambiando el estado termodinámico o la composición química. Para el primero, las simulaciones harían predicciones de observables a lo largo de las rutas de montaje. Para este último caso, el efecto de la perturbación en el proceso de montaje tendría que recapitularse en las simulaciones. Por ejemplo, los cambios de temperatura alteraron la velocidad y la vía de nucleación del clatrato21.

El descubrimiento de mecanismos guiados por máquinas30 integra fácilmente los avances en el aprendizaje automático aplicado a campos de fuerza19,31, muestreo32,33,34 y representación molecular19,31,35. El aumento del poder computacional y los avances en la inteligencia artificial simbólica permitirán a los algoritmos destilar descripciones matemáticas cada vez más precisas de los procesos complejos ocultos en datos de alta dimensión36. Como se muestra aquí, el muestreo guiado por máquina y la validación de modelos combinados con la regresión simbólica pueden respaldar el proceso de descubrimiento científico.

El comprometidor pB(x) es la probabilidad de que una trayectoria iniciada en la configuración X con velocidades de Maxwell-Boltzmann alcance el estado (meta)estable B antes de llegar a A. El disparo de trayectoria constituye, por tanto, un proceso de Bernoulli. Esperamos observar que las trayectorias nA y nB terminen en A y B, respectivamente, con probabilidad binomial \(p({n}_{{\mathrm{A}}},{n}_{{\mathrm{B}} }| {{{\bf{x}}}})=\binom{{{n}_{{\mathrm{A}}}+{n}_{{\mathrm{B}}}}}{{ {n}_{{\mathrm{A}}}}}{(1-{p}_{{\mathrm{B}}}({{{\bf{x}}}}))}^{{ n}_{{\mathrm{A}}}}{p}_{{\mathrm{B}}}{({{{\bf{x}}}})}^{{n}_{{\ matemáticas{B}}}}\). Para k puntos de disparo xi, la probabilidad combinada define la probabilidad \({{{\mathcal{L}}}}=\mathop{\prod }\nolimits_{i = 1}^{k}p({n}_{ {\mathrm{A}}}(i),{n}_{{\mathrm{B}}}(i)| {{{{\bf{x}}}}}_{i})\). Aquí ignoramos las correlaciones que surgen en transiciones rápidas dominadas por la inercia para trayectorias disparadas con velocidades iniciales opuestas11,18. Modelamos el confirmador desconocido con una función paramétrica y estimamos sus parámetros w maximizando la probabilidad \({{{\mathcal{L}}}}\) (refs. 5,15). Nos aseguramos de que 0 ≤ pB(x) ≤ 1 escribiendo el comprometidor en términos de una función de activación sigmoidea, \({p}_{{\mathrm{B}}}[q({{{\bf{x}} }}| {{{\bf{w}}}})]=1/(1+\exp [-q({{{\bf{x}}}}| {{{\bf{w}}} })])\). Aquí modelamos la probabilidad logarítmica q(x|w) usando una red neuronal15 y representamos la configuración con un vector x de características. Para N > 2 estados S, la distribución multinomial proporciona un modelo para p(n1, n2, . . . , nN|x), y escribir los confirmadores en los estados S en términos de la función de activación softmax garantiza la normalización, \(\mathop {\sum }\nolimits_{S = 1}^{N}{p}_{S}=1\). La función de pérdida l(w|θ) utilizada en el entrenamiento es el logaritmo negativo de la probabilidad \({{{\mathcal{L}}}}\).

TPS4,12 es un poderoso método MC de cadena de Markov en el espacio de ruta para muestrear TP. El movimiento de disparo bidireccional es una propuesta de movimiento eficiente en TPS4. Consiste en seleccionar aleatoriamente un punto de disparo XSP en el TP χ actual de acuerdo con la probabilidad psel (XSP | χ), dibujar velocidades aleatorias de Maxwell-Boltzmann y propagar dos trayectorias de prueba desde XSP hasta que alcancen cualquiera de los estados. Debido a que una de las trayectorias de prueba se propaga después de invertir primero todos los momentos en el punto inicial, es decir, se propaga hacia atrás en el tiempo, se puede construir un TP continuo si ambas pruebas terminan en estados diferentes. Dado un TP χ, un nuevo TP χ′ generado por disparos bidireccionales se acepta en la cadena de Markov con probabilidad37\({p}_{{{{\rm{acc}}}}}({\chi }^{ {\prime} }| \chi )=\min (1,{p}_{{{{\rm{sel}}}}}({{{{\bf{X}}}}}_{{{ {\rm{SP}}}}}| {\chi }^{{\prime} })/{p}_{{{{\rm{sel}}}}}({{{{\bf{X }}}}}_{{{{\rm{SP}}}}}| \chi ))\). Si se rechaza el nuevo camino, se repite χ.

Conociendo al comprometidor, es posible aumentar la velocidad a la que se generan los TP al sesgar la selección de puntos de disparo hacia el conjunto de estados de transición37, es decir, regiones con alta probabilidad reactiva p(TP|X). Para el caso de dos estados, esto equivale a sesgarse hacia la isosuperficie pB(x) = 1/2 que define los estados de transición con q(x) = 0. Construir un algoritmo que seleccione nuevos puntos de disparo sesgados hacia la mejor estimación actual. Para el conjunto del estado de transición y que aprende iterativamente a mejorar sus conjeturas en función de cada resultado de disparo recién observado, debemos equilibrar la exploración con la explotación. Para este fin, seleccionamos el nuevo punto de disparo X del TP actual χ usando una distribución de Lorentz centrada alrededor del conjunto de estados de transición, \({p}_{{{{\rm{sel}}}}}({{{ \bf{X}}}}| \chi )=1/\mathop{\sum}\limits_{{{{{\bf{x}}}}}^{{\prime} }\in \chi }[ (q{({{{\bf{x}}}})}^{2}+{\gamma }^{2})/(q{({{{{\bf{x}}}}}^ {{\prime} })}^{2}+{\gamma }^{2})],\) donde valores mayores de γ conducen a un aumento de la exploración. La distribución de Lorentz ofrece un equilibrio entre la eficiencia de la producción y la exploración ocasional fuera del estado de transición, que es necesaria para muestrear canales de reacción alternativos.

Con la función de comprometidor aprendida, se puede optimizar la definición de las fronteras estatales. Una definición de estado inicialmente estricta puede suavizarse moviendo los límites hacia afuera, digamos, pB(x) = 0,1 y pB(x) = 0,9. Esta relajación conduce a TP más cortos y acelera el muestreo.

La relación entre el comprometidor y la probabilidad de transición11 nos permite calcular el número esperado de TP generados al disparar desde una configuración X. Validamos el comprometidor aprendido sobre la marcha estimando el número esperado de transiciones antes de disparar desde una configuración y comparando con el resultado del disparo observado. El número esperado de transiciones \({n}_{{{{\rm{TP}}}}}^{\exp }\) calculado sobre una ventana que contiene los k intentos de disparos bidireccionales más recientes4 es \({n }_{{{{\rm{TP}}}}}^{\exp }=\mathop{\sum }\nolimits_{i = 1}^{k}2(1-{p}_{{\mathrm {B}}}({{{{\bf{x}}}}}_{i},i)){p}_{{\mathrm{B}}}({{{{\bf{x} }}}}_{i},i)\), donde pB(xi, i) es la estimación del comprometidor para el punto de disparo de prueba Xi en el paso i antes de observar el resultado del disparo. Iniciamos el aprendizaje cuando lo predicho (\({n}_{{{{\rm{TP}}}}}^{\exp }\)) y realmente generado (\({n}_{{{{\rm {TP}}}}}^{{{\rm{gen}}}}}\)) el número de TP difiere. Definimos un factor de eficiencia, \({\alpha }_{{{{\rm{eff}}}}}=\min (1,{(1-{n}_{{{{\rm{TP}} }}}^{{\mathrm{gen}} }/{n}_{{{{\rm{TP}}}}}^{{\mathrm{exp}} })}^{2})\) , donde un valor de cero indica una predicción perfecta (Datos ampliados, figura 9). Al entrenar sólo cuando es necesario, evitamos el sobreajuste. Aquí utilizamos αeff para escalar la tasa de aprendizaje en el algoritmo de descenso de gradiente. Además, no se realiza ningún entrenamiento si αeff está por debajo de un cierto umbral (especificado más adelante para cada sistema).

Los mecanismos moleculares se pueden describir en diferentes niveles de resolución. Se pueden utilizar muchas características de alta resolución que cuantifican propiedades locales o menos características de baja resolución que miden propiedades globales. Si bien las funciones de alta resolución tienden a estar fácilmente disponibles, la elección de funciones significativas de baja resolución depende de la comprensión física. Centrándonos en eventos moleculares raros, nuestro objetivo era organizar características en una jerarquía de resolución, desde coordenadas cartesianas de posiciones atómicas (la resolución más alta posible) hasta una sola cantidad, el comprometidor.

Diseñamos las redes neuronales en este estudio para alentarlas a aprender la jerarquía de resolución de las características. Las redes neuronales han demostrado la capacidad de aprender características de baja resolución a partir de otras de alta resolución, por ejemplo, cuando se utilizan en el reconocimiento de imágenes. Desde un punto de vista práctico, el ancho de capa de nuestras redes disminuye constantemente, pasando de las entradas a las salidas. Además, a medida que las características aprendidas se vuelven cada vez más globales (y por lo tanto menos redundantes) mientras van a capas más profundas, disminuimos la probabilidad de abandono al ascender en la red. Esto también se refleja en la diferente arquitectura utilizada para la formación del clatrato donde, debido a las características ya de grano grueso y específicas del sistema, utilizamos una red piramidal de alimentación directa comparativamente simple.

En cualquier proceso molecular específico, esperamos que sólo unos pocos de los muchos grados de libertad controlen realmente la dinámica de transición. Identificamos las entradas al modelo de compromiso que tienen el papel más importante en la determinación de su resultado después del entrenamiento. Para este fin, primero calculamos una pérdida de referencia, lref = l(w, θ), sobre el conjunto de entrenamiento no perturbado para comparar con los valores obtenidos perturbando cada entrada una por una38. Luego promediamos la pérdida \(l({{{\bf{w}}}},{\widetilde{{{{\bf{\uptheta }}}}}}_{i})\) sobre ≥100 perturbados conjuntos de entrenamiento \({\widetilde{{{{\bf{\uptheta }}}}}}_{i}\) con valores permutados aleatoriamente de la coordenada de entrada i en la dimensión del lote. La diferencia de pérdida promedio \({{\Delta }}{l}_{i}=\left\langle l({{{\bf{w}}}},{\widetilde{{{{\bf{\uptheta) }}}}}}_{i})\right\rangle -{l}_{{{{\rm{ref}}}}}\) es grande si la iésima entrada influye fuertemente en la salida del modelo entrenado, es decir, es relevante para predecir el cometidor.

En el subconjunto de baja dimensión que consta solo de las entradas z más relevantes (las que tienen el Δli más alto), la regresión simbólica genera expresiones matemáticas compactas que se aproximan al comprometidor completo. Nuestra implementación de regresión simbólica se basa en el paquete Python dcgpy39 y utiliza un algoritmo genético con una estrategia de evolución (N + 1). En cada generación, se generan N nuevas expresiones mediante cambios aleatorios en la estructura matemática de la expresión más adecuada de la generación principal. Posteriormente se utiliza una optimización basada en gradientes para encontrar los mejores parámetros para cada expresión. Luego se elige la expresión más adecuada como padre para la siguiente generación. La idoneidad de cada expresión de prueba pB(z) se mide por \({l}_{{{{\rm{sr}}}}}({{{{\bf{w}}}}}_{{{ {\rm{sr}}}}}| {{{\bf{\uptheta }}}})\equiv -\log {{{\mathcal{L}}}}[{p}_{{\mathrm{ B}}}({{{{\bf{z}}}}}_{{{{\rm{sp}}}}})]+\lambda C\), donde agregamos el término de regularización λC al log-verosimilitud (consulte 'Estimación de máxima verosimilitud de la función de compromiso') para mantener las expresiones simples y evitar el sobreajuste. Aquí λ > 0 y C es una medida de la complejidad de la expresión de prueba, estimada en nuestro caso por el número de operaciones matemáticas.

La regresión simbólica producirá expresiones de diferente complejidad según la fuerza de la regularización. Seleccionamos expresiones con un equilibrio razonable entre simplicidad y precisión utilizando un diagrama de Pareto (Fig. 2c), en el que trazamos la complejidad (medida como el número de operaciones matemáticas) frente a la precisión (medida como la pérdida de datos de validación). . Al escanear un rango de valores de λ, podemos identificar modelos en el frente de Pareto para su posterior análisis.

Investigamos la formación de pares de iones monovalentes en agua para evaluar la capacidad del algoritmo para aprender con precisión el comprometidor de transiciones que están fuertemente influenciadas por los grados de libertad del solvente. Utilizamos seis configuraciones de sistemas diferentes (LiCl, LiI, NaCl, NaI, CsCl y CsI), cada una de las cuales consta de un catión y un anión en agua.

Todas las simulaciones MD se llevaron a cabo en cajas de simulación cúbicas utilizando el campo de fuerza de Joung y Cheatham40 junto con agua TIP3P41. Cada caja de simulación contenía un único par de iones solvatados con 370 moléculas de agua TIP3P. Usamos el motor openMM MD42 para propagar las ecuaciones de movimiento en pasos de tiempo de Δt = 2 fs con un integrador de velocidad Verlet con aleatorización de velocidad43 del paquete openmmtools de Python. Después de un equilibrio inicial de NPT a presión constante P = 1 bar y temperatura constante T = 300 K, todas las simulaciones de producción se realizaron en el conjunto NVT a un volumen constante V y una temperatura constante de T = 300 K. La fricción se estableció en 1 ps-1. Las interacciones no enlazadas se calcularon utilizando un esquema de Ewald de malla de partículas44 con un límite de espacio real de 1 nm y una tolerancia de error de 0,0005. Cada simulación de TPS (que consta de simulaciones MD y entrenamiento de redes neuronales) se llevó a cabo en medio nodo utilizando una unidad central de procesamiento (CPU) Xeon Gold 6248 junto con una GPU RTX5000. En TPS, los estados asociados y disociados se definieron de acuerdo con las distancias interiónicas (consulte la Tabla complementaria 4 para conocer los valores de cada especie iónica).

El cometidor de una configuración es invariante bajo traslaciones y rotaciones globales en ausencia de campos externos, y además es invariante con respecto a permutaciones de partículas idénticas. Por lo tanto, optamos por transformar el sistema de coordenadas del espacio cartesiano a una representación que incorpore las simetrías físicas del comprometidor. Para lograr una transformación casi sin pérdidas, utilizamos la distancia interiónica para describir el soluto y adaptamos funciones de simetría para describir la configuración del disolvente45. Las funciones de simetría se desarrollaron originalmente para describir estructuras moleculares en potenciales de redes neuronales19,46, pero también se han utilizado con éxito para detectar y describir diferentes fases del hielo en simulaciones atomísticas47. Estas funciones describen el entorno que rodea a un átomo central sumando todas las partículas idénticas a una distancia radial determinada. El tipo de función de simetría \({G}_{i}^{2}\) cuantifica la densidad de las moléculas de disolvente alrededor de un átomo de soluto i en una capa centrada en rs

donde la suma abarca todos los átomos de disolvente j de un tipo de átomo específico, rij es la distancia entre el átomo central i y el átomo j, y η controla el ancho de la capa. La función fc(r) es un límite de Fermi definido como:

lo que asegura que la contribución de los átomos distantes del solvente desaparezca. El parámetro escalar αc controla la inclinación del corte. El tipo de función de simetría \({G}_{i}^{5}\) analiza adicionalmente la distribución angular del disolvente alrededor del átomo central i.

donde la suma abarca todos los pares de átomos del disolvente distintos, ϑijk es el ángulo formado entre los dos átomos del disolvente y el átomo central del soluto, el parámetro ζ es un número par que controla la nitidez de la distribución angular, y λ = ±1 establece el ubicación del mínimo con respecto a ϑijk en π y 0, respectivamente. Consulte la Tabla complementaria 5 para conocer las combinaciones de parámetros utilizadas. Escalamos todas las entradas para que se encuentren aproximadamente en el rango [0, 1] para aumentar la estabilidad numérica del entrenamiento. En particular, normalizamos las funciones de simetría dividiéndolas por el número promedio esperado de átomos (o pares de átomos) para una distribución isotrópica en el volumen de sondeo.

Las funciones de simetría del tipo G2 cuentan el número de átomos de disolvente en el volumen de sondeo; la constante de normalización \({\langle N[{G}_{i}^{2}]\rangle }_{{{{\rm{iso}}}}}\) es, por lo tanto, el número esperado de átomos en la volumen de sondeo \({V}_{{{{\rm{sonda}}}}}^{(2)}\)

donde ρN es la densidad numérica promedio del tipo de átomo de solvente analizado. El volumen de sondeo exacto para el tipo G2 se puede aproximar como

para η pequeño y rcut > rs.

Las funciones de tipo G5 incluyen un término angular adicional y cuentan el número de pares de átomos de disolvente ubicados en lados opuestos del átomo de soluto central. El número esperado de pares \({\langle {N}_{{{{\rm{pairs}}}}}\rangle }_{{{{\rm{iso}}}}}\) se puede calcular a partir de el número esperado de átomos en el volumen sondeado \({\langle {N}_{{{{\rm{átomos}}}}}\rangle }_{{{{\rm{iso}}}}}\) como \({\langle {N}_{{{{\rm{átomos}}}}}\rangle }_{{{{\rm{iso}}}}}({\langle {N}_{{ {{\rm{átomos}}}}}\rangle }_{{{{\rm{iso}}}}}-1)/2\). Esta expresión solo es exacta para valores enteros de \({\langle {N}_{{{{\rm{átomos}}}}}\rangle }_{{{{\rm{iso}}}}}\) e incluso puede volverse negativo si \({\langle {N}_{{{{\rm{átomos}}}}}\rangle }_{{{{\rm{iso}}}}} < 1\). Por lo tanto, utilizamos una aproximación que se garantiza que no será negativa.

El número esperado de átomos \({\langle {N}_{{{{\rm{átomos}}}}}\rangle }_{{{{\rm{iso}}}}}\) se puede calcular a partir de el volumen que se sondea para un átomo de soluto fijo y con un átomo de disolvente fijo

Con la expectativa de que la mayoría de los grados de libertad del sistema no controlen la transición, diseñamos redes neuronales que filtran progresivamente entradas irrelevantes y construyen una función altamente no lineal de las restantes. Probamos tres arquitecturas de redes neuronales piramidales diferentes, 'ResNet I', 'ResNet II' y 'SNN', donde los nombres que contienen 'ResNet' indican el uso de unidades residuales48,49 y 'SNN' es una arquitectura de red neuronal autonormalizada50 (ver Tablas complementarias 6–8 para las arquitecturas exactas utilizadas). La arquitectura con mejor rendimiento es 'ResNet I' (consulte el archivo de datos complementario 1 para comparar el rendimiento de las diferentes arquitecturas para todos los sistemas iónicos). ResNet Utilicé una pila piramidal de cinco unidades residuales de preactivación, cada una con cuatro capas ocultas. El número de unidades ocultas por capa se reduce mediante un factor constante \(f={(10/{n}_{{\mathrm{in}}})}^{1/4}\) después de cada bloque de unidades residuales y disminuye de nin = 221 en la primera unidad a 10 en la última. Además, después de cada bloque residual se aplica una eliminación de 0,1fi, donde i es el índice unitario residual que oscila entre 0 y 4. La optimización de los pesos de la red se realiza utilizando el descenso de gradiente de Adam51. Para todas las arquitecturas, el entrenamiento se realizó después de cada tercer paso de TPS MC durante una época con una tasa de aprendizaje de lr = αeff10-3, si lr ≥ 10-4. El factor de eficiencia esperado αeff se calculó en una ventana de k = 100 pasos de TPS. Realizamos todo el aprendizaje profundo con código escrito personalizado basado en Keras52. Las simulaciones de TPS se llevaron a cabo utilizando una versión personalizada de openpathsampling53,54 junto con nuestro propio módulo Python.

Para el entrenamiento de transferencia, la última capa con una sola neurona (es decir, el predictor logarítmico) de un modelo originalmente entrenado en LiCl se aleatorizó y todos los demás pesos se mantuvieron fijos durante el entrenamiento posterior en los datos de disparo para las otras especies iónicas ( LiI, NaCl, NaI, CsCl y CsI). El entrenamiento se realizó utilizando el optimizador Adam con una tasa de aprendizaje de lr = 2,5 × 10−5. La pérdida de prueba se calculó después de cada época en el 20% de los datos utilizados como conjunto de prueba. El entrenamiento se detuvo cuando no se observó ninguna disminución en la pérdida de prueba durante más de 1000 épocas. Luego se utilizó el modelo con la mejor pérdida de prueba.

Para la extrapolación en el espacio químico (Datos ampliados, figura 3), configuramos una red neuronal de múltiples iones de arquitectura 'ResNet I'. El modelo fue entrenado en los resultados del disparo para diferentes pares de especies iónicas simultáneamente, como se especifica. Utilizó las coordenadas del conjunto 'SF shortranged' junto con los parámetros ϵ y σ Lennard-Jones del campo de fuerza para distinguir las diferentes especies iónicas. El entrenamiento se realizó con el optimizador Adam (lr = 10-3) utilizando el 10% de los datos como conjunto de prueba. El entrenamiento finalizaba si la pérdida de prueba no disminuía durante 1000 épocas y luego se utilizaba el modelo con la mejor pérdida de prueba.

Seleccionamos las siete coordenadas más relevantes identificadas por la red neuronal de iones múltiples como entradas para las regresiones simbólicas de iones múltiples (Tablas complementarias 9 a 12). Usamos entre 3 y 7 de estas coordenadas más relevantes para ejecuciones de regresión simbólica independientes usando los valores de regularización λ = 0,001, λ = 0,0001 y λ = 0,00001. Luego seleccionamos la expresión informada en la Fig. 2d usando el diagrama de Pareto en la Fig. 2c.

También seleccionamos las cinco coordenadas más relevantes identificadas a partir de una red neuronal entrenada en LiCl para ejecuciones de regresión simbólica (Datos extendidos, figuras 4b-d). Regularizamos las expresiones producidas penalizando el número total de operaciones matemáticas elementales con λ = 10−6 y λ = 10−7.

Las contribuciones de cada átomo al comprometidor en un sistema X particular (Fig. 1b) se calcularon como la magnitud del gradiente de la coordenada de reacción q (x) con respecto a sus coordenadas cartesianas. Todas las magnitudes de gradiente se escalaron con la masa atómica inversa.

El algoritmo de aprendizaje se aplicó a un conjunto de datos TPS existente de nucleación de clatrato de metano producido inicialmente para la referencia. 21. Contiene datos para simulaciones realizadas a cuatro temperaturas diferentes T = 270 K, 275 K, 280 K y 285 K (consulte la Tabla complementaria 14 para obtener más detalles). Se realizaron nuevas simulaciones para obtener los valores de compromiso muestreados utilizados en la validación. Todas las simulaciones del confirmador se realizaron con OpenMM 7.1.142 en GPU NVIDIA GeForce GTX TITAN 1080Ti, disparando entre 6 y 18 trayectorias por configuración utilizando el mismo protocolo de simulación que en la ref. 21.

Utilizamos 22 características diferentes para describir el tamaño, la cristalinidad, la estructura y la composición del núcleo cristalino de hidrato de metano en crecimiento (Tabla complementaria 1). Además de las características que describen configuraciones moleculares, utilizamos la temperatura como entrada para las redes neuronales y la regresión simbólica. En una red piramidal de avance con 9 capas, reducimos el número de unidades por capa de 23 en la entrada a una en la última capa (Tabla complementaria 15). La red se entrenó con el optimizador Adam con una tasa de aprendizaje lr = 10-3 en los datos TPS existentes para todas las temperaturas, dejando fuera el 10% de los puntos de disparo como datos de prueba. Detuvimos el entrenamiento después de que la pérdida en el conjunto de prueba no disminuyó durante 10,000 épocas y utilizamos el modelo con la mejor pérdida de prueba. Todo el entrenamiento de la red neuronal se realizó en una GPU RTX6000. Usamos las tres coordenadas más relevantes como entradas para ejecuciones de regresión simbólica con una penalización en el número total de operaciones matemáticas elementales usando λ = 10−5.

Aplicamos nuestro algoritmo de aprendizaje automático a los datos de disparo existentes de cristalización de polímeros24,25. Usamos dos conjuntos diferentes de características para describir la transición, un conjunto de 35 características de baja resolución (de grano grueso) que también se usó en trabajos anteriores y un conjunto de características de alta resolución que describen cada perla de polímero por sí sola. Las características de baja resolución contienen una serie de medidas globales como la energía potencial U y los parámetros de orden de enlace de Steinhardt Q4 y Q6, descripciones del entorno local de partículas de polímero seleccionadas, varias medidas que describen la estructura del polímero mediante el conteo de cadenas y bucles y algunas distancias seleccionadas (consulte la Tabla complementaria 16 para obtener una lista exhaustiva). El conjunto de características de alta resolución consta del número de conexiones, vecinos y los parámetros de orden de enlace de Lechner-Dellago Steinhardt promediados por vecinos55 para cada perla de polímero, es decir, cada configuración corresponde a un vector de características con 3 × 128 = 384 entradas.

Tanto para la descripción de alta resolución como para la de baja resolución, utilizamos redes neuronales en forma de pirámide (Tablas complementarias 17 y 18). En ambos casos, el entrenamiento se realizó utilizando el método de descenso de gradiente de Adam con una tasa de aprendizaje lr = 10-3 utilizando el 20% de los datos como datos de prueba. Los modelos se guardaron y se calculó la pérdida de prueba después de cada época. El entrenamiento se detuvo si la pérdida de la prueba no disminuía durante 10.000 épocas. Luego se utilizó el modelo con la menor pérdida de prueba como modelo entrenado final. Todo el entrenamiento de la red neuronal se realizó en una GPU RTX6000.

Utilizamos entre dos y cinco de las cinco características de baja resolución más relevantes como entradas en ejecuciones de regresión simbólica (Tablas complementarias 19 a 22). Regularizamos penalizando el número de operaciones matemáticas elementales con λ = 10−2, 10−3, 10−4, 10−5 y 10−6.

Utilizamos el campo de fuerza de Martini de grano grueso (v2.2)56,57,58,59 para describir el ensamblaje del homodímero transmembrana alfa-helicoidal Mga2. Todas las simulaciones de MD se llevaron a cabo con GROMACS v4.6.760,61,62,63 con un paso de tiempo de integración de Δt = 0,02 ps, utilizando una caja de simulación cúbica que contiene dos hélices alfa idénticas de 30 aminoácidos de longitud en una membrana lipídica hecha de 313 moléculas de 1-palmitoil-2-oleoil-glicero-3-fosfocolina (POPC). La membrana abarca la caja en el plano x-y y se solvató con agua (5348 perlas de agua) e iones NaCl correspondientes a una concentración de 150 mM (58 Na+, 60 Cl-). Se aplicó una temperatura de referencia de T = 300 K utilizando el termostato v-rescale64 con una constante de acoplamiento de 1 ps por separado en la proteína, la membrana y el disolvente. Se aplicó una presión de 1 bar por separado en el plano x-y y en z utilizando un barostato semiisotrópico Parrinello-Rahman65 con una constante de tiempo de 12 ps y una compresibilidad de 3 × 10-4 bar-1. Cada simulación MD se llevó a cabo en un único nodo informático con dos CPU E5-2680-v3 y 64 GB de memoria. Todo el entrenamiento de la red neuronal se realizó en una GPU RTX6000.

Para describir el ensamblaje del homodímero Mga2, utilizamos 28 distancias interhelicoidales por pares entre las perlas principales de las dos hélices junto con el número total de contactos interhelicoidales, la distancia entre los centros de masa de las hélices y una serie de características hechas a mano que describen el Organización de lípidos alrededor de las dos hélices (Tabla complementaria 2). Para garantizar que todas las entradas de la red se encuentren aproximadamente en [0, 1], utilizamos la función sigmoidal \(f(r)={(1-(r/{R}_{0})^{6})}/{ (1-(r/{R}_{0})^{12})}\) con R0 = 2 nm para todas las distancias por pares, mientras que escalamos todas las características de los lípidos utilizando los valores mínimos y máximos tomados a lo largo de la transición. Los estados ensamblado y desensamblado se definen como configuraciones con ≥130 contactos interhelicoidales y con distancias entre centro de masa hélice-hélice dCoM ≥ 3 nm, respectivamente.

La red neuronal utilizada para ajustar el committor se implementó usando Keras52 y consistió en una parte piramidal inicial de 3 capas en la que el número de unidades disminuye de las 36 entradas a 6 en la última capa usando un factor constante de (6/36)1 /2 seguido de 6 unidades residuales48,49, cada una con 4 capas y 6 neuronas por capa (Tabla complementaria 23). Se aplica una abandono de 0,01 a las entradas y la red se entrena utilizando el protocolo de descenso de gradiente Adam con una tasa de aprendizaje de lr = 0,0001 (ref. 51).

Para investigar el mecanismo de ensamblaje de Mga2, realizamos un muestreo guiado por máquina en paralelo en múltiples nodos de un grupo de computadoras de alto rendimiento. Ejecutamos 500 cadenas TPS independientes guiadas por el modelo de compromiso actual. Las simulaciones de 500 TPS se inicializaron con TP iniciales aleatorios. La red neuronal utilizada para seleccionar los puntos de disparo iniciales se entrenó en intentos de disparo preliminares (8.044 disparos independientes desde 1.160 puntos diferentes). Después de dos rondas (dos pasos en cada una de las 500 cadenas de TPS independientes), actualizamos el modelo de compromiso entrenando con todos los datos nuevos. Volvimos a entrenar después de la sexta ronda. No se requirió más capacitación, como lo indican los números consistentes de recuentos esperados y observados de TP. Realizamos otras 14 rondas para las 500 cadenas de TPS para recolectar TP (Fig. 5b). La selección del punto de disparo, la configuración del TPS y el entrenamiento de la red neuronal se automatizaron completamente en código Python utilizando MDAnalysis66,67, numpy68 y nuestro paquete Python personalizado.

El análisis de importancia de los insumos reveló que el número total de contactos ncontacts es el insumo más importante (Datos ampliados, figura 7). Sin embargo, ninguna expresión generada por la regresión simbólica en función de ncontactos únicamente fue precisa para reproducir al autor del compromiso. Es probable que la red entrenada utilice ncontacts solo como un interruptor binario para distinguir los dos regímenes diferentes: cerca de los estados vinculados o no vinculados. Al restringir el análisis de importancia de entrada a puntos de entrenamiento cercanos al estado libre, encontramos que la red utiliza varios contactos interhelicoidales que aproximadamente siguen un patrón helicoidal (Datos extendidos, figura 7). Realizamos una regresión simbólica en todas las combinaciones posibles realizadas por una, dos o tres de las siete coordenadas de entrada más importantes (Tabla complementaria 3). Las mejores expresiones en términos de pérdida se seleccionaron utilizando datos del confirmador de validación que no se habían utilizado durante la optimización. Este conjunto de validación constaba de datos del confirmador para 516 configuraciones con 30 disparos de prueba cada una y 32 configuraciones con 10 disparos de prueba.

Para evaluar la variabilidad en los mecanismos de reacción observados, realizamos una agrupación jerárquica de todos los TP proyectados en el plano definido por los contactos 9 y 22, que ingresan a la parametrización más precisa generada por la regresión simbólica. Luego utilizamos la distorsión dinámica del tiempo69 para calcular la similitud por pares entre todos los TP para la agrupación, que realizamos utilizando el módulo de agrupación scipy70,71. Los gráficos de densidad de rutas (Fig. 5f, g) se histogramaron de acuerdo con la cantidad de rutas, no con la cantidad de configuraciones, es decir, el contador de cada celda visitada por una ruta en particular se incrementó en uno para esta ruta.

Los datos y archivos del conjunto de entrenamiento para configurar simulaciones de dinámica molecular para el ensamblaje de LiCl se incluyen en la cápsula Code Ocean72. Los datos para reproducir este estudio para todos los sistemas restantes (todos los iones restantes, polímero, clatrato y dímero transmembrana MGA2) están disponibles públicamente en un repositorio de Zenodo73. Los datos originales se proporcionan con este documento.

Una versión ejecutable del código 'Inteligencia artificial para el descubrimiento de mecanismos moleculares' (AIMMD) está disponible en la cápsula del software Code Ocean72. El código AIMMD también está disponible en https://github.com/bio-phys/aimmd como repositorio git.

Pena-Francesch, A., Jung, H., Demirel, MC & Sitti, M. Materiales biosintéticos autorreparantes para máquinas blandas. Nat. Madre. 19, 1230-1235 (2020).

Artículo de Google Scholar

Van Driessche, AES et al. Mecanismos de nucleación molecular y estrategias de control para la selección de polimorfos cristalinos. Naturaleza 556, 89–94 (2018).

Artículo de Google Scholar

Chung, HS, Piana-Agostinetti, S., Shaw, DE y Eaton, WA Origen estructural de la difusión lenta en el plegamiento de proteínas. Ciencia 349, 1504-1510 (2015).

Artículo de Google Scholar

Dellago, C., Bolhuis, PG y Chandler, D. Muestreo eficiente de rutas de transición: aplicación a reordenamientos de grupos de Lennard-Jones. J. química. Física. 108, 9236–9245 (1998).

Artículo de Google Scholar

Peters, B. & Trout, BL Obtención de coordenadas de reacción mediante maximización de probabilidad. J. química. Física. 125, 054108 (2006).

Artículo de Google Scholar

Bolhuis, PG, Dellago, C. y Chandler, D. Coordenadas de reacción de isomerización biomolecular. Proc. Acad. Nacional. Ciencia. Estados Unidos 97, 5877–5882 (2000).

Artículo de Google Scholar

Best, RB y Hummer, G. Coordenadas y tasas de reacción de rutas de transición. Proc. Acad. Nacional. Ciencia. Estados Unidos 102, 6732–6737 (2005).

Artículo de Google Scholar

Berezhkovskii, AM & Szabo, A. Difusión a lo largo de la coordenada de reacción de probabilidad de división/compromiso. J. Física. Química. B 117, 13115-13119 (2013).

Artículo de Google Scholar

E, W. & Vanden-Eijnden, E. Hacia una teoría de los caminos de transición. J. estadística. Física. 123, 503 (2006).

Artículo MathSciNet MATEMÁTICAS Google Scholar

Krivov, SV Reducción óptima de la dimensionalidad de dinámicas complejas: el juego de ajedrez como difusión en un paisaje de energía libre. Física. Rev. E 84, 011135 (2011).

Artículo de Google Scholar

Hummer, G. De caminos de transición a estados de transición y coeficientes de tasa. J. química. Física. 120, 516–523 (2003).

Artículo de Google Scholar

Bolhuis, PG, Chandler, D., Dellago, C. & Geissler, PL Muestreo de caminos de transición: arrojar cuerdas sobre pasos de montaña accidentados, en la oscuridad. Año. Rev. Phys. Química. 53, 291–318 (2002).

Artículo de Google Scholar

LeCun, Y., Bengio, Y. y Hinton, G. Aprendizaje profundo. Naturaleza 521, 436–444 (2015).

Artículo de Google Scholar

Schmidhuber, J. Aprendizaje profundo en redes neuronales: una descripción general. Red neuronal. 61, 85-117 (2015).

Artículo de Google Scholar

Ma, A. & Dinner, AR Método automático para identificar coordenadas de reacción en sistemas complejos. J. Física. Química. B 109, 6769–6779 (2005).

Artículo de Google Scholar

Vanden-Eijnden, E., Venturoli, M., Ciccotti, G. y Elber, R. Sobre los supuestos subyacentes a los hitos. J. química. Física. 129, 174102 (2008).

Artículo de Google Scholar

Schmidt, M. & Lipson, H. Destilación de leyes naturales de forma libre a partir de datos experimentales. Ciencia 324, 81–85 (2009).

Artículo de Google Scholar

Ballard, AJ y Dellago, C. Hacia el mecanismo de disociación iónica en el agua. J. Física. Química. B 116, 13490–13497 (2012).

Artículo de Google Scholar

Behler, J. & Parrinello, M. Representación generalizada de redes neuronales de superficies de energía potencial de alta dimensión. Física. Rev. Lett. 98, 146401 (2007).

Artículo de Google Scholar

Walsh, MR, Koh, CA, Sloan, ED, Sum, AK y Wu, DT Simulaciones de microsegundos de nucleación y crecimiento espontáneos de hidrato de metano. Ciencia 326, 1095–1098 (2009).

Artículo de Google Scholar

Arjun, Berendsen, TA y Bolhuis, PG Visión atomística imparcial de los mecanismos de nucleación competitivos de los hidratos de metano. Proc. Acad. Nacional. Ciencia. Estados Unidos 116, 19305–19310 (2019).

Wang, Y., Herron, L. y Tiwary, P. De los datos al ruido y a los datos para mezclar la física en diferentes temperaturas con inteligencia artificial generativa. Proc. Acad. Nacional. Ciencia. Estados Unidos 119, e2203656119 (2022).

Artículo de Google Scholar

Jacobson, LC, Hujo, W. & Molinero, V. Precursores amorfos en la nucleación de hidratos de clatrato. Mermelada. Química. Soc. 132, 11806–11811 (2010).

Artículo de Google Scholar

Leitold, C. & Dellago, C. Mecanismo de plegado de una cadena de polímero con atracciones de corto alcance. J. química. Física. 141, 134901 (2014).

Artículo de Google Scholar

Leitold, C., Lechner, W. y Dellago, C. Una coordenada de reacción de cuerda para el plegamiento de una cadena de polímero. J. Física. Condensa. Asunto 27, 194126 (2015).

Artículo de Google Scholar

Covino, R. y col. Un sensor eucariota para la saturación de lípidos de membrana. Mol. Celda 63, 49–59 (2016).

Artículo de Google Scholar

Chiavazzo, E. et al. Exploración de la dinámica de mapas intrínsecos en busca de paisajes de energía libre efectivos e inexplorados. Proc. Acad. Nacional. Ciencia. Estados Unidos 114, E5494–E5503 (2017).

Artículo de Google Scholar

Bittracher, A. y col. Múltiples de transición de sistemas metaestables complejos. J. Ciencia no lineal. 28, 471–512 (2018).

Artículo MathSciNet MATEMÁTICAS Google Scholar

Lucente, D., Duffner, S., Herbert, C., Rolland, J. y Bouchet, F. Aprendizaje automático de funciones de compromiso para predecir eventos climáticos de alto impacto. Preimpresión en arXiv https://doi.org/10.48550/arXiv.1910.11736 (2019).

Wang, Y., Lamim Ribeiro, JM & Tiwary, P. Enfoques de aprendizaje automático para analizar y mejorar simulaciones de dinámica molecular. actual. Opinión. Estructura. Biol. 61, 139-145 (2020).

Artículo de Google Scholar

Noé, F., Tkatchenko, A., Müller, K.-R. & Clementi, C. Aprendizaje automático para simulación molecular. Año. Rev. Phys. Química. 71, 361–390 (2020).

Artículo de Google Scholar

Noé, F., Olsson, S., Köhler, J. & Wu, H. Generadores de Boltzmann: muestreo de estados de equilibrio de sistemas de muchos cuerpos con aprendizaje profundo. Ciencia 365, eaaw1147 (2019).

Rogal, J., Schneider, E. y Tuckerman, ME Variables colectivas de ruta basadas en redes neuronales para un muestreo mejorado de transformaciones de fase. Física. Rev. Lett. 123, 245701 (2019).

Artículo de Google Scholar

Sidky, H., Chen, W. y Ferguson, AL Aprendizaje automático para el descubrimiento colectivo de variables y muestreo mejorado en simulación biomolecular. Mol. Física. 118, e1737742 (2020).

Artículo de Google Scholar

Bartók, AP et al. El aprendizaje automático unifica el modelado de materiales y moléculas. Ciencia. Adv. 3, e1701816 (2017).

Artículo de Google Scholar

Udrescu, S.-M. & Tegmark, M. AI Feynman: un método inspirado en la física para la regresión simbólica. Ciencia. Adv. 6, fácil2631 (2020).

Artículo de Google Scholar

Jung, H., Okazaki, K.-i y Hummer, G. Muestreo de rutas de transición de eventos raros disparando desde arriba. J. química. Física. 147, 152716 (2017).

Artículo de Google Scholar

Kemp, SJ, Zaradic, P. y Hansen, F. Un enfoque para determinar la importancia y el significado relativos de los parámetros de entrada en redes neuronales artificiales. Ecológico. Modelo. 204, 326–334 (2007).

Artículo de Google Scholar

Izzo, D. & Biscani, F. dcgp: programación genética cartesiana diferenciable simplificada. J. Software de código abierto. 5, 2290 (2020).

Artículo de Google Scholar

Joung, IS y Cheatham, TE Determinación de parámetros de iones monovalentes alcalinos y haluro para su uso en simulaciones biomoleculares explícitamente solvatadas. J. Física. Química. B 112, 9020–9041 (2008).

Artículo de Google Scholar

Jorgensen, WL, Chandrasekhar, J., Madura, JD, Impey, RW y Klein, ML Comparación de funciones potenciales simples para simular agua líquida. J. química. Física. 79, 926–935 (1983).

Artículo de Google Scholar

Eastman, P. y col. OpenMM 7: rápido desarrollo de algoritmos de alto rendimiento para dinámica molecular. Computación PLoS. Biol. 13, e1005659 (2017).

Sivak, DA, Chodera, JD & Crooks, GE El reescalado por pasos de tiempo recupera propiedades dinámicas de tiempo continuo para la integración de Langevin en tiempo discreto de sistemas sin equilibrio. J. Física. Química. B 118, 6466–6474 (2014).

Essmann, U. et al. Un método de Ewald de malla de partículas suave. J. química. Física. 103, 8577–8593 (1995).

Artículo de Google Scholar

Behler, J. Funciones de simetría centradas en el átomo para la construcción de potenciales de redes neuronales de alta dimensión. J. química. Física. 134, 074106 (2011).

Artículo de Google Scholar

Behler, J. Representación de superficies de energía potencial mediante potenciales de redes neuronales de alta dimensión. J. Física. Condensa. Asunto 26, 183001 (2014).

Artículo de Google Scholar

Geiger, P. & Dellago, C. Redes neuronales para la detección de estructuras locales en sistemas polimórficos. J. química. Física. 139, 164105 (2013).

Artículo de Google Scholar

He, K., Zhang, X., Ren, S. & Sun, J. Aprendizaje residual profundo para el reconocimiento de imágenes. En la Conferencia IEEE de 2016 sobre visión por computadora y reconocimiento de patrones 770–778 (IEEE, 2016).

Él, K., Zhang, X., Ren, S. y Sun, J. Mapeos de identidad en redes residuales profundas. En Computer Vision—ECCV 2016, Lecture Notes in Computer Science (eds Leibe, B. et al.) 630–645 (Springer, 2016).

Klambauer, G., Unterthiner, T., Mayr, A. y Hochreiter, S. Redes neuronales autonormalizadas. En Proc. 31ª Conferencia Internacional sobre Sistemas de Procesamiento de Información Neuronal 972–981 (Curran Associates, 2017).

Kingma, DP & Ba, J. Adam: un método de optimización estocástica. Preimpresión en arXiv https://doi.org/10.48550/arXiv.1412.6980 (2017).

Chollet, F.Keras. https://github.com/fchollet/keras (2015).

Swenson, DWH, Prinz, J.-H., Noe, F., Chodera, JD & Bolhuis, PG openpathsampling: un marco de Python para simulaciones de muestreo de rutas. 1. Conceptos básicos. J. química. Computación teórica. 15, 813–836 (2019).

Artículo de Google Scholar

Swenson, DWH, Prinz, J.-H., Noe, F., Chodera, JD & Bolhuis, PG openpathsampling: un marco de Python para simulaciones de muestreo de rutas. 2. Crear y personalizar conjuntos de rutas y esquemas de muestra. J. química. Computación teórica. 15, 837–856 (2019).

Artículo de Google Scholar

Lechner, W. y Dellago, C. Determinación precisa de estructuras cristalinas basada en parámetros de orden de enlaces locales promediados. J. química. Física. 129, 114707 (2008).

Artículo de Google Scholar

Marrink, SJ, de Vries, AH y Mark, AE Modelo de grano grueso para simulaciones semicuantitativas de lípidos. J. Física. Química. B 108, 750–760 (2004).

Artículo de Google Scholar

Marrink, SJ, Risselada, HJ, Yefimov, S., Tieleman, DP y de Vries, AH El campo de fuerza MARTINI: modelo de grano grueso para simulaciones biomoleculares. J. Física. Química. B 111, 7812–7824 (2007).

Artículo de Google Scholar

Monticelli, L. y col. El campo de fuerza de grano grueso MARTINI: extensión a las proteínas. J. química. Computación teórica. 4, 819–834 (2008).

Artículo de Google Scholar

de Jong, DH et al. Parámetros mejorados para el campo de fuerza de proteína de grano grueso de martini. J. química. Computación teórica. 9, 687–697 (2013).

Artículo de Google Scholar

Berendsen, H., van der Spoel, D. y van Drunen, R. GROMACS: una implementación de dinámica molecular paralela de transmisión de mensajes. Computadora. Física. Comunitario. 91, 43–56 (1995).

Artículo de Google Scholar

Hess, B., Kutzner, C., van der Spoel, D. y Lindahl, E. GROMACS 4: algoritmos para una simulación molecular escalable, con carga equilibrada y altamente eficiente. J. química. Computación teórica. 4, 435–447 (2008).

Artículo de Google Scholar

Pronk, S. y col. GROMACS 4.5: un conjunto de herramientas de simulación molecular de código abierto, altamente paralelo y de alto rendimiento. Bioinformática 29, 845–854 (2013).

Artículo de Google Scholar

Abraham, MJ y cols. GROMACS: simulaciones moleculares de alto rendimiento mediante paralelismo multinivel desde portátiles hasta supercomputadoras. SoftwareX 1–2, 19–25 (2015).

Artículo de Google Scholar

Bussi, G., Donadio, D. y Parrinello, M. Muestreo canónico mediante reescalado de velocidad. J. química. Física. 126, 014101 (2007).

Artículo de Google Scholar

Parrinello, M. & Rahman, A. Transiciones polimórficas en monocristales: un nuevo método de dinámica molecular. J. Aplica. Física. 52, 7182–7190 (1981).

Artículo de Google Scholar

Michaud-Agrawal, N., Denning, EJ, Woolf, TB y Beckstein, O. MDAnalysis: un conjunto de herramientas para el análisis de simulaciones de dinámica molecular. J. Computación. Química. 32, 2319–2327 (2011).

Artículo de Google Scholar

Gowers, R. y col. MDAnalysis: un paquete Python para el análisis rápido de simulaciones de dinámica molecular. Actas de las conferencias Python in Science (eds Benthall, S. & Rostrup, S.) 98–105 (2016); https://doi.org/10.25080/issn.2575-9752

Harris, CR y cols. Programación de matrices con NumPy. Naturaleza 585, 357–362 (2020).

Artículo de Google Scholar

Meert, W., Hendrickx, K. y Van Craenendonck, T. Wannesm/dtaidistance v2.0.0. Zenodo https://zenodo.org/record/3276100 (2020)

Virtanen, P. y col. SciPy 1.0: algoritmos fundamentales para la computación científica en Python. Nat. Métodos 17, 261–272 (2020).

Artículo de Google Scholar

Müllner, D. Algoritmos modernos de agrupamiento aglomerativo y jerárquico. Preimpresión en arXiv https://doi.org/10.48550/arXiv.1109.2378 (2011).

Jung, H. y col. Muestreo de rutas guiado por máquina para descubrir mecanismos de autoorganización molecular (cápsula de software). Código Océano https://doi.org/10.24433/CO.7949737.v1 (2023).

Jung, H. y col. Muestreo de rutas guiado por máquina para descubrir mecanismos de autoorganización molecular (datos de capacitación y validación). Zenodo https://doi.org/10.5281/zenodo.7704969 (2023).

Descargar referencias

Agradecemos a FE Blanc por sus útiles comentarios y a la comunidad de openpathsampling, en particular a D. Swenson, por las discusiones y el soporte técnico. HJ, RC y GH agradecen el apoyo de la Sociedad Max Planck. RC agradece el apoyo del Instituto de Estudios Avanzados de Frankfurt. RC y GH reconocen el apoyo del programa LOEWE CMMS del estado de Hesse, el Centro de Investigación Colaborativa 1507 de la Fundación Alemana de Investigación y la Escuela Internacional de Investigación Max Planck en Biofísica Celular. AA y PGB reconocen el apoyo del programa CSER de la Organización Holandesa para la Investigación Científica (NWO) y de Shell Global Solutions International. BV, CL y CD agradecen el apoyo del Fondo Austriaco para la Ciencia (FWF) a través del SFB TACO, subvención número F 81-N. Los financiadores no tuvieron ningún papel en el diseño del estudio, la recopilación y análisis de datos, la decisión de publicar o la preparación del manuscrito.

Financiamiento de acceso abierto proporcionado por la Sociedad Max Planck.

Estos autores contribuyeron igualmente: Hendrik Jung, Roberto Covino.

Departamento de Biofísica Teórica, Instituto Max Planck de Biofísica, Frankfurt am Main, Alemania

Hendrik Jung y Gerhard Hummer

Instituto de Estudios Avanzados de Frankfurt, Frankfurt am Main, Alemania

Roberto Covino

Instituto van 't Hoff de Ciencias Moleculares, Universidad de Ámsterdam, Ámsterdam, Países Bajos

A. Arjun y Peter G. Bolhuis

Facultad de Física, Universidad de Viena, Viena, Austria

Christian Leitold y Christoph Dellago

Instituto de Biofísica, Universidad Goethe de Frankfurt, Frankfurt am Main, Alemania

Gerhard Hummer

También puedes buscar este autor en PubMed Google Scholar.

También puedes buscar este autor en PubMed Google Scholar.

También puedes buscar este autor en PubMed Google Scholar.

También puedes buscar este autor en PubMed Google Scholar.

También puedes buscar este autor en PubMed Google Scholar.

También puedes buscar este autor en PubMed Google Scholar.

También puedes buscar este autor en PubMed Google Scholar.

HJ desarrolló el código de aprendizaje automático, realizó simulaciones de ensamblaje de iones y proteínas y analizó los datos junto con RC y GHAA realizó las simulaciones de hidratos de gas y analizó los datos junto con HJ, RC y PGBCL realizaron las simulaciones de plegamiento de polímeros y analizaron los datos juntos. con HJ, RC y CDGH concibieron el estudio. HJ, RC y GH diseñaron la investigación. HJ, RC y GH escribieron el manuscrito con aportaciones de todos los autores. RC, PGB, CD y GH coordinaron el proyecto.

Correspondencia a Gerhard Hummer.

Los autores declaran no tener conflictos de intereses.

Editora de manejo principal: Kaitlin McCardle, en colaboración con el equipo de Nature Computational Science.

Nota del editor Springer Nature se mantiene neutral con respecto a reclamos jurisdiccionales en mapas publicados y afiliaciones institucionales.

Función de autocorrelación normalizada del tiempo de transición tTP en función del número de pasos de Monte Carlo (MC) para simulaciones de TPS sobre diferentes especies iónicas, en las que los puntos de disparo fueron seleccionados por el algoritmo (curvas roja y verde, corridas independientes) o por selección aleatoria uniforme (curva naranja y azul, ejecuciones independientes). Cada función de autocorrelación se calculó sobre cadenas MC con una longitud de 100.000 pasos cada una. Los tiempos de descorrelación de las cadenas de Markov, definidos como los puntos donde las funciones de autocorrelación normalizadas alcanzan el valor e−1, están marcados con líneas verticales de puntos del mismo color que la función de autocorrelación respectiva. El muestreo guiado por máquina conduce a tiempos de descorrelación de cuatro pasos de MC en todos los casos, mientras que para la selección aleatoria uniforme este valor oscila entre 20 y 60 pasos dependiendo de la especie iónica.

Datos fuente

Gráficos de correlación cruzada de comitadores para modelos obtenidos mediante entrenamiento por transferencia en diferentes especies iónicas con un número variable de puntos de disparo de entrenamiento. El modelo original se entrenó en LiCl y es el mismo en todos los casos (consulte la Fig. 1 para ver un gráfico de correlación cruzada del confirmador del modelo original que predice LiCl). El modelo es de arquitectura 'ResNet I' y utiliza las entradas 'SF longranged II' (consulte las Tablas complementarias 5 y 6). El entrenamiento de transferencia se realizó aleatorizando la última capa del modelo original con una sola neurona y luego entrenando en el número especificado de puntos de disparo con datos para cada sistema iónico por separado. El promedio de los confirmadores muestreados (línea azul) +/- un SD (sombreado en naranja) se calculó agrupando el confirmador muestreado en el rango del confirmador aprendido indicado por los pasos verticales. Como referencia, la línea roja indica la identidad.

Datos fuente

Se muestran las correlaciones cruzadas de los comprometidores muestreados (eje x) y los comprometidores previstos (eje y) para modelos entrenados en diferentes subconjuntos de los seis pares de iones diferentes. (Fila 1) Entrenamiento en los seis pares de iones (Li+ Cl−, Li+ I−, Na+ Cl−, Na+ I−, Cs+ Cl− y Cs+ I−) simultáneamente. (Fila 2) Modelo entrenado sin datos para Li+, es decir, excluyendo Li+ Cl− y Li+ I−. (Fila 3) Modelo entrenado sin datos para Na+, es decir, excluyendo Na+ Cl− y Na+ I−. (Fila 4) Modelo entrenado sin datos para Cl−, es decir, excluyendo Li+ Cl−, Na+ Cl− y Cs+ Cl−. (Fila 5) Modelo entrenado sin datos para I−, es decir, excluyendo Li+ I−, Na+ I− y Cs+ I−. Todos los modelos son de arquitectura 'ResNet I' (consulte la Tabla complementaria 6). Para todos los modelos, las diferentes especies iónicas se distinguieron agregando los parámetros ϵ y σ Lennard-Jones del campo de fuerza como descriptores adicionales al conjunto 'SF de corto alcance' (consulte la Tabla complementaria 5), ​​que se usó para describir el solvente alrededor del par de iones. El promedio de los confirmadores muestreados (línea azul) +/- un SD (sombreado en naranja) se calculó agrupando el confirmador muestreado en el rango del confirmador aprendido indicado por los pasos verticales. Como referencia, la línea roja indica la identidad.

Datos fuente

a, Validación del modelo reducido de iones múltiples para cuatro pares de iones adicionales. El promedio de los confirmadores muestreados (línea azul) y su desviación estándar (sombreado en naranja) se calculan agrupando a lo largo del confirmador previsto (línea roja: identidad). b, Coordenadas de entrada más importantes que determinan el comprometidor entrenado en simulaciones de asociación de un solo par iónico, en este caso LiCl (consulte también la Tabla complementaria 13 para obtener una lista de las diez coordenadas más relevantes). c, Modelos reducidos q0 y q1 que describen la asociación de Li+ Cl− en agua obtenidos mediante regresión simbólica con regularización estricta (λ = 10−6) y suave (λ = 10−7), respectivamente. Tenga en cuenta que el primer modelo no depende de los grados de libertad del agua (consulte nuevamente la Tabla complementaria 13 para obtener una descripción de las coordenadas). d, Gráficos de correlación cruzada entre los datos del comprometidor no entrenado y las predicciones de regresión simbólica como validaciones independientes de la precisión de q0 y q1. El promedio de los confirmadores muestreados (línea azul) +/- un SD (sombreado en naranja) se calcula agrupando a lo largo del confirmador previsto (línea roja: identidad).

Datos fuente

Correlación cruzada entre el comprometidor aprendido y el comprometido obtenido mediante muestreo repetido en configuraciones no entrenadas para dos modelos que están entrenados en solo tres de las cuatro temperaturas disponibles en el conjunto de entrenamiento. (Izquierda) Modelo de compromiso entrenado solo con datos para T=270 K, 275 K y 285 K para evaluar la capacidad del modelo para interpolar a T=280 K. (Derecha) Modelo entrenado con datos para T=270 K, 275 K y 280 K para evaluar la capacidad del modelo para extrapolar a T = 285 K. La línea roja representa la identidad.

Datos fuente

a, Correlación cruzada del confirmador en configuraciones no entrenadas. b, Relevancia de entrada para redes neuronales entrenadas con datos de nucleación de polímeros utilizando características de baja y alta resolución, respectivamente. En a, el promedio de los confirmadores muestreados (línea azul) y su desviación estándar (sombreado en naranja) se calcularon agrupando a lo largo del confirmador previsto (línea roja: identidad). En b, las barras azules muestran la media de n repeticiones independientes del análisis de importancia de entrada (n = 100 para funciones de alta resolución y n = 250 para funciones de baja resolución), y las barras de error indican +/- una DE.

Datos fuente

Análisis de importancia de entrada utilizando todos los puntos de entrenamiento (panel superior) y un subconjunto con ncontacts < 0,01 (panel inferior), correspondiente a puntos de entrenamiento cercanos al estado independiente. La altura de cada barra es la media de 50 análisis independientes (n = 50), mientras que las barras indican +/- una DE. Todos los valores están normalizados a la mayor importancia en cada conjunto.

Datos fuente

Gráfico de correlación cruzada del comitador para la expresión de regresión simbólica \({q}_{B}({x}_{9},{x}_{2}2)=-\exp ({x}_{9}^{2) })\log ({x}_{9}-\frac{{x}_{9}}{\log ({x}_{22})})\) en datos del confirmador de validación no capacitado del ensamblaje transmembrana Mga2. La expresión es una función del contacto interhelicoidal 9 (x9) en la parte superior de las dos hélices y del contacto 22 (x22) en la parte inferior (consulte las Tablas complementarias 2 y 3). El promedio de los confirmadores muestreados (línea azul) y su desviación estándar (sombreado en naranja) se calculan agrupando a lo largo del confirmador previsto (línea roja: identidad).

Datos fuente

Iteraciones de entrenamiento para el ensamblaje Li+ Cl−. La línea azul muestra la tasa de aprendizaje calculada a partir del factor de eficiencia en cada paso, las cruces naranjas muestran cuándo ocurrió realmente el entrenamiento. El recuadro muestra la pérdida de entrenamiento por punto de tiro para cada entrenamiento. Sólo se muestran los primeros 26.000 escalones de Montecarlo. El modelo utiliza las entradas 'SF longranged II' y tiene la arquitectura 'ResNet I' (consulte las Tablas complementarias 5 y 6).

Datos fuente

Tablas complementarias 1 a 23.

Eficiencia de las simulaciones de TPS de iones: número de transiciones aceptadas, generadas y esperadas por paso de Monte Carlo para todas las simulaciones de TPS asistidas por aprendizaje automático en todas las especies iónicas. La diferencia que se muestra en la última columna es entre las transiciones esperadas y generadas, es decir, (esperado − generado)/pasos. Los valores para la fase inicial (final) se calcularon sobre los primeros (últimos) pasos de 1000 MC, el valor para la simulación completa es sobre los 100.000 pasos de MC.

Datos de fuente estadística para las figuras 1b-e.

Datos de fuente estadística para la Fig. 2a, c, e.

Datos de fuente estadística para la Fig. 3b,c.

Datos de fuente estadística para la Fig. 4a, b.

Datos de fuente estadística para las figuras 5b, c, e – g.

Fuente de datos estadísticos para datos extendidos Fig. 1.

Fuente de datos estadísticos para datos extendidos Fig. 2.

Fuente de datos estadísticos para datos extendidos Fig. 3.

Datos de fuente estadística para datos extendidos Fig. 4a, b, d.

Fuente de datos estadísticos para datos extendidos Fig. 5.

Fuente de datos estadísticos para datos extendidos Fig. 6.

Fuente de datos estadísticos para datos extendidos Fig. 7.

Fuente de datos estadísticos para datos extendidos Fig. 8.

Fuente de datos estadísticos para datos extendidos Fig. 9.

Acceso Abierto Este artículo está bajo una Licencia Internacional Creative Commons Attribution 4.0, que permite el uso, compartir, adaptación, distribución y reproducción en cualquier medio o formato, siempre y cuando se dé el crédito apropiado a los autores originales y a la fuente. proporcione un enlace a la licencia Creative Commons e indique si se realizaron cambios. Las imágenes u otro material de terceros en este artículo están incluidos en la licencia Creative Commons del artículo, a menos que se indique lo contrario en una línea de crédito al material. Si el material no está incluido en la licencia Creative Commons del artículo y su uso previsto no está permitido por la normativa legal o excede el uso permitido, deberá obtener permiso directamente del titular de los derechos de autor. Para ver una copia de esta licencia, visite http://creativecommons.org/licenses/by/4.0/.

Reimpresiones y permisos

Jung, H., Covino, R., Arjun, A. et al. Muestreo de rutas guiado por máquina para descubrir mecanismos de autoorganización molecular. Nat Comput Sci 3, 334–345 (2023). https://doi.org/10.1038/s43588-023-00428-z

Descargar cita

Recibido: 15 de febrero de 2023

Aceptado: 10 de marzo de 2023

Publicado: 24 de abril de 2023

Fecha de emisión: abril de 2023

DOI: https://doi.org/10.1038/s43588-023-00428-z

Cualquier persona con la que comparta el siguiente enlace podrá leer este contenido:

Lo sentimos, actualmente no hay un enlace para compartir disponible para este artículo.

Proporcionado por la iniciativa de intercambio de contenidos Springer Nature SharedIt