Resultados y discusión
El evento de autoionización se investigó utilizando simulaciones ab initio RETIS como se describe en Materiales y Métodos. Para las simulaciones RETIS, utilizamos un parámetro de orden de distancia geométrica relativamente simple, λ, como se ilustra en la Fig. 1: Cuando el sistema consiste en sólo especies de H2O, λ es la mayor distancia de enlace covalente O-H, y cuando el sistema contiene especies de OH- y H3O+, λ se toma como la distancia más corta entre el oxígeno en OH- y los átomos de hidrógeno en H3O+. En lo sucesivo, nos referiremos al átomo de oxígeno utilizado para el parámetro de orden como Oλ. El tipo de especie (OH-, H2O o H3O+) se identificó asignando a cada hidrógeno un único enlace que lo conecta con el oxígeno más cercano. Nótese que la definición del parámetro de orden no requiere un umbral para definir un enlace químico ni restringe el parámetro de orden a moléculas de agua específicas durante la duración de la simulación. Esto significa que calculamos la tasa de disociación de cualquier molécula de agua en el sistema en lugar de un único enlace O-H objetivo o molécula de agua.
Condiciones de inicio y variables colectivas locales. (A) Distribuciones reactivas (rλc,λr(ξ)) y no reactivas (uλc,λr(ξ)) para ξ={w4,na} y λc=1,16 Å y λr=2,0 Å. Para su visualización, las distribuciones representadas están normalizadas . Los insertos superior y derecho muestran las proyecciones unidimensionales de las distribuciones. Puede verse una clara separación de las dos distribuciones a lo largo de la coordenada w4, lo que indica que las trayectorias reactivas están más comprimidas en comparación con las no reactivas. Además, el átomo de oxígeno utilizado en el cálculo del parámetro de orden (Oλ) acepta de media un mayor número de enlaces de hidrógeno en las trayectorias reactivas, en comparación con las no reactivas. (B) Instantánea ilustrativa de una trayectoria reactiva en la que Oλ se muestra en azul. Los cuatro átomos de oxígeno circundantes que se utilizan para el cálculo del parámetro de orden tetraédrico q se muestran en naranja. El hilo de agua se destaca con una línea amarilla (y esferas grises transparentes) y se indica el parámetro de ángulo qcos. En esta instantánea, el hilo de agua está comprimido, q muestra una desviación de una estructura tetraédrica, qcos indica que tres átomos de oxígeno se alinean en el hilo, y Oλ acepta tres enlaces de hidrógeno y dona uno (mostrado con líneas verdes).
Si consideramos la coordenada q, observamos que rλc,λr se desplaza hacia valores q más bajos en comparación con uλc,λr, lo que indica que una distorsión de una disposición tetraédrica alrededor de la especie de agua que se disocia también puede iniciar el evento. Este hallazgo es algo sorprendente ya que en algunas otras reacciones químicas en fase acuosa se encontró el efecto contrario (31). Se pueden extraer conclusiones similares para la distribución de ξ=(w4,qcos). Aquí, hay un pico a lo largo de la coordenada qcos para la distribución reactiva más cercana a una disposición lineal de las moléculas de agua. En la Fig. 4B mostramos una instantánea representativa, obtenida al principio (después de 3 fs) de una trayectoria reactiva. En general, los resultados mostrados en la Fig. 3 informan de que una compresión del hilo de agua (medida por w4) y la hipercoordinación (medida por na) o la distorsión (medida por q y qcos) son condiciones de iniciación necesarias para la autoionización. Sin embargo, no son condiciones suficientes, como muestran los valores de TAλc,λr en la Fig. 3B: todavía el 60% de las trayectorias que parten del rango ideal de parámetros ξ no consiguen establecer un salto de protón concertado.
El aprendizaje automático (ML) aplicado a los datos de muestreo de trayectorias (33, 34) es un enfoque prometedor para encontrar variables colectivas importantes que la intuición humana puede pasar por alto fácilmente. Para explorar esta posibilidad, construimos modelos ML para predecir el resultado de las trayectorias dado el estado del sistema de agua al principio de las mismas. Nos centramos en el mismo rango que en el análisis del poder predictivo y utilizamos el estado del sistema, cuando se alcanza por primera vez λ>1,15 Å, para predecir el resultado. Utilizamos varias técnicas de ML en las que cada conjunto de trayectorias impares se incluyó en la calibración y los conjuntos de trayectorias pares se utilizaron para el conjunto de pruebas. Una división alternativa en la que los datos dentro de cada conjunto de trayectorias se dividieron uniformemente en dos dio resultados similares. Además, como las distribuciones muy sesgadas son difíciles de tratar con ML, omitimos además la reponderación de los conjuntos de datos con los pesos estadísticos de los conjuntos de trayectorias correspondientes. Sin embargo, aplicamos las técnicas de ML como un enfoque cualitativo para encontrar nuevos parámetros que pudieran ser probados cuantitativamente dentro del método de poder predictivo (19).
Además, para evitar un riesgo potencial de sobreinterpretación optamos por restringir la complejidad del proceso de decisión de ML e impusimos un máximo de cuatro parámetros de orden al calcular TAλc,λr. Por ejemplo, se obtuvieron excelentes rendimientos predictivos (>90%) utilizando las máquinas de gradiente-boosting basadas en conjuntos (35, 36). Sin embargo, la interpretación del modelo es problemática, ya que se utiliza un conjunto de 100-150 árboles de decisión profundos (añadidos en una secuencia). Aunque el rendimiento mejora, aumenta la posibilidad de sobreajuste con correlaciones accidentales. Por ello, nos hemos limitado a los modelos de decisión basados en un solo árbol de clasificación y regresión (CART) (20). La restricción a cuatro parámetros de orden para la función TAλc,λr se basa en razones similares. Si se añaden más parámetros se obtienen matrices más dispersas que representan las distribuciones reactivas/no reactivas y, como resultado, la integración numérica para calcular el solapamiento entre estas distribuciones se vuelve muy sensible al tamaño de la bandeja y podría subestimar el solapamiento debido a que las bandejas están vacías por una estadística insuficiente.
Consideramos 138 variables colectivas que consisten en las distancias oxígeno-oxígeno; las distancias oxígeno-hidrógeno para las moléculas de agua inicialmente unidas; todos los ángulos formados por Oλ y sus cuatro vecinos de oxígeno más cercanos; y los parámetros de orden de Steinhardt de los órdenes 3, 4 y 6 (32) (ver Materiales y Métodos para más detalles). Además, se añadieron los parámetros de orden ya considerados. La Fig. 5A muestra el árbol de decisión resultante. Sorprendentemente, de todos los parámetros de entrada, el parámetro w4 se encuentra en la cima del árbol de decisión y es la variable más importante, medida por la reducción del error de clasificación atribuido a cada variable en cada división del árbol de decisión (20) (Apéndice SI, Fig. S9). También el ordenamiento tetraédrico y el número de enlaces de hidrógeno aceptados aparecen en el árbol de decisión. Para describir el primer efecto, el enfoque ML priorizó el parámetro de ordenación q4 de Steinhardt por encima del parámetro q similar utilizado previamente por nosotros. Algunas distancias que también aparecen en el árbol de decisión como d25, la distancia entre Oλ y su 25º oxígeno más cercano, se deben muy probablemente a correlaciones accidentales causadas por el tamaño limitado del conjunto de datos. Esto se comprueba inspeccionando la importancia de esta variable: d25 no aparece entre las 20 variables más importantes (Apéndice SI, Fig. S9), y, de hecho, otras variables similares (p. ej., d24) ocupan un lugar más alto, aunque con poca importancia. Un parámetro más importante e intuitivo que sugiere el enfoque ML es λ2, la distancia OH entre el oxígeno más cercano a Oλ y su hidrógeno con el mayor enlace intramolecular. Al volver a calcular la capacidad de predicción utilizando los parámetros del árbol ML (Fig. 5B) no se obtuvieron rendimientos más elevados que la combinación w4, q, na y qcos, pero debería concebirse como igualmente buena, teniendo en cuenta las incertidumbres estadísticas.
Resultados del análisis de aprendizaje automático. (A) Árbol de clasificación y regresión para predecir el resultado de las trayectorias iniciadas. Aquí, consideramos varias variables colectivas adicionales (descripción en Materiales y Métodos), pero sólo un pequeño subconjunto es finalmente necesario para construir el árbol: w4, q4 , λ2 (la longitud del enlace de hidrógeno estirado en la molécula de agua más cercana a la especie Oλ), di (la distancia del Oλ al ith oxígeno más cercano), y d¯i (la distancia media considerando los i oxígenos más cercanos). La notación de los nodos se explica con el nodo independiente de la esquina superior izquierda. Este árbol predice que las trayectorias sean reactivas, es decir, que alcancen un λ≥2, o no reactivas basándose en las variables colectivas obtenidas en el fotograma de las trayectorias cuando λ es primero ≥1,15. Los nodos que predicen trayectorias reactivas están coloreados en azul (clase 1) mientras que los nodos que predicen trayectorias no reactivas están coloreados en verde (clase 0). Obsérvese que los porcentajes de la parte inferior de los cuadros no reflejan las fracciones físicamente correctas, ya que los conjuntos de trayectorias no se reponen utilizando sus pesos estadísticos. Las reglas son representaciones textuales del recorrido del árbol; por ejemplo, la regla 5 (que predice trayectorias reactivas) puede expresarse como w4≥7,6 y λ2≥1,1. Estas reglas dan diferentes condiciones de iniciación, y se enumeran en el Apéndice SI, Tabla S1, para la fila inferior de nodos. (B) El poder predictivo y la probabilidad de cruce en función de λr para λc=1,16 Å y diferentes combinaciones de variables colectivas. Aquí comparamos el poder predictivo utilizando las variables colectivas que identificamos con las variables marcadas como importantes por el análisis de aprendizaje automático. (C) Distribuciones reactivas (rλc,λr(ξ)) y no reactivas (uλc,λr(ξ)) para ξ={λ2,d¯2} y λc=1,16 Å y λr=2,0 Å. Para su visualización, las distribuciones representadas están normalizadas. Los insertos superior y derecho muestran las proyecciones unidimensionales de las distribuciones.