Blog

Dec 11, 2023

Maschine

Nature Computational Science Band 3, Seiten 334–345 (2023)Diesen Artikel zitieren 8339 Zugriffe 8 Zitate 44 Altmetrische Metriken Einzelheiten zur molekularen Selbstorganisation, angetrieben durch konzertierte Vielteilchen

Nature Computational Science Band 3, Seiten 334–345 (2023)Diesen Artikel zitieren

8339 Zugriffe

8 Zitate

44 Altmetrisch

Details zu den Metriken

Die durch konzertierte Vielteilcheninteraktionen angetriebene molekulare Selbstorganisation erzeugt die geordneten Strukturen, die sowohl unbelebte als auch lebende Materie definieren. Hier präsentieren wir einen autonomen Pfad-Sampling-Algorithmus, der Deep Learning und Übergangspfadtheorie integriert, um den Mechanismus molekularer Selbstorganisationsphänomene zu entdecken. Der Algorithmus nutzt das Ergebnis neu initiierter Trajektorien, um quantitative mechanistische Modelle zu erstellen, zu validieren und – falls erforderlich – zu aktualisieren. Zum Abschluss des Lernzyklus leiten die Modelle die Probenahme, um die Probenahme seltener Montageereignisse zu verbessern. Die symbolische Regression verdichtet den erlernten Mechanismus in einer für den Menschen interpretierbaren Form im Hinblick auf relevante physikalische Observablen. Angewandt auf die Ionenassoziation in Lösung, die Bildung von Gashydratkristallen, die Polymerfaltung und den Membran-Protein-Zusammenbau erfassen wir die Vielteilchen-Lösungsmittelbewegungen, die den Zusammenbauprozess steuern, identifizieren die Variablen der klassischen Nukleationstheorie und decken den Faltungsmechanismus auf verschiedenen Ebenen auf Auflösung und offenbaren konkurrierende Montagepfade. Die mechanistischen Beschreibungen sind auf thermodynamische Zustände und den chemischen Raum übertragbar.

Das Verständnis, wie generische, aber subtil orchestrierte Wechselwirkungen bei der Bildung komplexer Strukturen zusammenwirken, ist der Schlüssel zur Steuerung der molekularen Selbstorganisation1,2. Als Computerexperimente versprechen uns Simulationen der Molekulardynamik (MD) atomar detaillierte und unvoreingenommene Einblicke in Selbstorganisationsprozesse3. Die meisten kollektiven Selbstorganisationsprozesse sind jedoch seltene Ereignisse, die auf Zeitskalen ablaufen, die um viele Größenordnungen länger sind als die schnellen molekularen Bewegungen, die den MD-Integrationsschritt begrenzen. Das System verbringt die meiste Zeit in metastabilen Zuständen, und die seltenen und schnellen stochastischen Übergänge zwischen Zuständen werden in unvoreingenommenen MD-Simulationen, wenn überhaupt, selten gelöst. Diese Übergangspfade (TPs) sind die ganz besonderen Trajektoriensegmente, die den Reorganisationsprozess erfassen. Das Erlernen molekularer Mechanismen aus Simulationen erfordert die Konzentration der Rechenleistung auf die Probenahme von TPs4 und die Destillation quantitativer Modelle daraus5. Aufgrund der hohen Dimensionalität des Konfigurationsraums sind sowohl die Probenentnahme als auch die Informationsextraktion in der Praxis äußerst anspruchsvoll. Unser Algorithmus bewältigt beide Herausforderungen gleichzeitig. Es erstellt autonom und gleichzeitig quantitative mechanistische Modelle komplexer molekularer Ereignisse, validiert die Modelle im laufenden Betrieb und nutzt sie, um die Probenentnahme im Vergleich zur regulären MD um Größenordnungen zu beschleunigen.

Die statistische Mechanik bietet einen allgemeinen Rahmen, um niedrigdimensionale mechanistische Modelle von Selbstorganisationsereignissen zu erhalten. In diesem Artikel konzentrieren wir uns auf Systeme, die sich zwischen zwei Zuständen A und B (zusammengebaut bzw. zerlegt) neu organisieren, aber eine Verallgemeinerung auf eine beliebige Anzahl von Zuständen ist einfach. Jedes TP, das die beiden Zustände verbindet, enthält eine Folge von Snapshots, die das System während seiner Reorganisation erfassen. Folglich ist das Transition Path Ensemble (TPE) der Mechanismus mit der höchsten Auflösung. Da der Übergang praktisch stochastisch ist, erfordert die Quantifizierung seines Mechanismus einen probabilistischen Rahmen. Wir definieren den Committor pS(x) als die Wahrscheinlichkeit, dass eine Trajektorie zuerst in den Zustand S eintritt, mit S = A bzw. B, wobei x ein Vektor von Merkmalen ist, die den Startpunkt X im Konfigurationsraum darstellen, und pA(x) + pB(x) = 1 für ergodische Dynamik. Der Committor pB berichtet über den Fortschritt der Reaktion A → B und sagt das Schicksal der Flugbahn auf Markovsche Weise voraus6,7, was sie zur idealen Reaktionskoordinate8,9 macht. Beim Schachspiel kann man sich den Committor beispielsweise als die Wahrscheinlichkeit vorstellen, dass Schwarz in wiederholten Partien für gegebene Anfangspositionen auf dem Brett gewinnt10. Die Mindestanforderungen für Anwendungen über molekulare Simulationen hinaus bestehen darin, (1) dass eine Größe vorhanden ist, die einem Committor ähnelt, und (2) dass die Dynamik des Systems wiederholt abgetastet werden kann, zumindest in Vorwärtsrichtung. Die Wahrscheinlichkeit verschiedener möglicher Ereignisse (A, B, …) sollte daher zumindest teilweise kodiert sein (und somit im Hinblick auf den momentanen Zustand X des Systems erlernbar sein) und die Dynamik des Systems sollte vorzugsweise einer wiederholten Abtastung zugänglich sein durch effiziente Computersimulation. Wenn man jedoch wiederholt ein reales System mit zufriedenstellender Kontrolle über die Anfangsbedingungen vorbereiten kann, kann man mithilfe unseres Rahmenwerks lernen, das wahrscheinliche Schicksal dieses Systems anhand der beobachteten und kontrollierten Anfangsbedingungen vorherzusagen.

Das Abtasten von TPs für seltene Ereignisse und das Erlernen der zugehörigen Committor-Funktion pB(x) sind zwei herausragende und untrennbar miteinander verbundene Herausforderungen. Angesichts der Tatsache, dass TPs in einem sehr hochdimensionalen Raum äußerst selten sind, ist eine uninformierte Suche zwecklos. Allerdings konvergieren TPs in der Nähe von Übergangszuständen7, wo die Trajektorie mit pA(x) = pB(x) = 1/2 am wenigsten festgelegt ist. Für die Markovsche und zeitreversible Dynamik P(TP|x) erfüllt die Wahrscheinlichkeit, dass eine durch x verlaufende Trajektorie ein TP ist, P(TP|x) = 2pB(x)(1 − pB(x)), d. h , bestimmt der Committer die Wahrscheinlichkeit, einen TP11 abzutasten. Die Herausforderungen der Informationsextraktion und der Probenentnahme sind somit miteinander verknüpft.

Um diese doppelten Herausforderungen zu bewältigen, haben wir einen iterativen Algorithmus entwickelt, der im Sinne des Reinforcement Learning auf der Transition Path Theory9 und dem Transition Path Sampling (TPS)12 aufbaut. Der Algorithmus lernt durch wiederholte Durchführung virtueller Experimente den Beteiligten seltener Ereignisse in komplexen Vielteilchensystemen und nutzt den Wissensgewinn, um die Probenahme von TPs zu verbessern (Abb. 1a). In jedem Experiment wählt der Algorithmus einen Punkt Um ein detailliertes Gleichgewicht der TPS sicherzustellen, wählt der Algorithmus Strukturen aus dem aktuellen Übergangspfad aus, um Schussbewegungen mit neu gezeichneten Maxwell-Boltzmann-Anfangsgeschwindigkeiten zu initiieren. Nach wiederholten Schüssen von verschiedenen Punkten Nur wenn die Vorhersage schlecht ist, trainiert der Algorithmus das Modell anhand der Ergebnisse aller virtuellen Experimente neu, wodurch eine Überanpassung verhindert wird. Mit zunehmender Vorhersagekraft des mechanistischen Modells wird der Algorithmus bei der Abtastung von TPs effizienter, indem er Anfangspunkte in der Nähe von Übergangszuständen wählt, d. h. gemäß P(TP|x).

a, Mechanismuslernen durch Pfadabtastung. Die Methode iteriert zwischen der Abtastung von Übergangspfaden aus einer Konfiguration x zwischen den metastabilen Zuständen A und B (links) und dem Erlernen des Committors pB(x) (rechts). Eine neuronale Netzwerkfunktion aus molekularen Merkmalen (x1 bis x4) modelliert den Committor. Der Log-Prädiktor, der die letzte Schicht bildet, wird nicht angezeigt. Bei der Konvergenz destilliert die symbolische Regression einen interpretierbaren Ausdruck, der den molekularen Mechanismus anhand ausgewählter Merkmale (x1, x2) und numerischer Konstanten (a, c) quantifiziert, die durch mathematische Operationen (hier: +, −, ×, exp) verbunden sind. b, Schnappschüsse entlang eines TP, die die Bildung eines LiCl-Ionenpaars (von rechts nach links) in einer atomistischen MD-Simulation zeigen. Wasser wird als Stäbchen dargestellt, Li+ als kleine Kugel und Cl− als große Kugel. Atome werden entsprechend ihrem Beitrag zum Reaktionsfortschritt von niedrig (blau) nach hoch (rot) gefärbt, quantifiziert durch ihren Beitrag zum Gradienten der Reaktionskoordinate q(x|w). c, Selbstkonsistenz. Zählt die generierte (blaue Linie) und erwartete (orange gestrichelte Linie) Anzahl von Übergangsereignissen. Die grüne Linie zeigt die kumulative Differenz zwischen der beobachteten und der erwarteten Anzahl. Der Einschub zeigt eine Vergrößerung der ersten 1.000 Iterationen. d, Validierung des gelernten Committors. Kreuzkorrelation zwischen dem vom trainierten Netzwerk vorhergesagten Committor und dem Committor, der durch wiederholte Stichproben aus molekularen Konfigurationen erhalten wurde, für die das Committor-Modell nicht trainiert wurde. Der Durchschnitt der Stichproben-Committors (blaue Linie) und ihr SD (orange schattiert) wurden in den durch die vertikalen Schritte angegebenen Bins des erlernten Committors berechnet. Als Referenz zeigt die rote Linie die Identität an. e, Übertragbarkeit des gelernten Committors. Darstellung des Transferlernens und Kreuzkorrelationen zwischen ausgewählten Committoren für die NaCl- und NaI-Ionenpaarung und Vorhersagen des Committors aus einem Modell, das auf Daten für LiCl trainiert und durch Transferlernen unter Verwendung von jeweils nur 1.000 zusätzlichen Schussergebnissen angepasst wurde. Farben und SD (angezeigt durch orange Schattierung) sind wie in d.

Quelldaten

Der Algorithmus lernt aus seinen wiederholten Versuchen, indem er Deep Learning5,13,14 auf selbstkonsistente Weise nutzt. Wir modellieren den Committor pB(x) = 1/(1 + e−q(x|w)) mit einem neuronalen Netzwerk15 q(x|w) der Gewichte w. Bei dieser Wahl ist q eine invertierbare Funktion von pB, und wir können beide Funktionen als Reaktionskoordinaten betrachten. Bei jedem Versuch, ein TP zu generieren, propagiert der Algorithmus zwei Trajektorien, eine zeitlich vorwärts und eine rückwärts, indem er MD-Simulationen durchführt4. Im Erfolgsfall geht zunächst eine Trajektorie in den Zustand A und die andere in den Zustand B über und bildet so gemeinsam ein neues TP. Bei diesem Bernoulli-Münzwurfprozess lernt der Algorithmus sowohl aus Erfolgen als auch aus Misserfolgen. Die negative Log-Likelihood5 für k Versuche definiert die Verlustfunktion \(l({{{\bf{w}}}}| {{{\bf{\uptheta }}}})=\mathop{\sum }\nolimits_ {i = 1}^{k}\log (1+{{\mathrm{e}}}^{{s}_{i}q({{{{\bf{x}}}}}_{i }| {{{\bf{w}}}})})\), wobei si = 1, wenn die von Xi ausgehende Flugbahn i zuerst in A eintritt, und si = −1, wenn sie zuerst in B eintritt. Der Trainingssatz θ enthält die k Merkmalsvektoren xi, die den Schusspunkten Xi und den Ergebnissen si zugeordnet sind. Indem wir das Netzwerk q(x|w) trainieren, um den Verlust l(w|θ) zu minimieren, erhalten wir eine Maximum-Likelihood-Schätzung des Committors5, die differenzierbar ist und anspruchsvolle Analysemethoden ermöglicht16.

Anschließend verwenden wir die symbolische Regression17, um den molekularen Mechanismus in eine für den Menschen interpretierbare Form zu verdichten (Abb. 1a) und physikalische Erkenntnisse zu gewinnen. Zunächst identifiziert eine Attributanalyse des trainierten Netzwerks eine kleine Teilmenge z der Eingabekoordinaten x, die das Ergebnis der Netzwerkvorhersage dominieren. Dann destilliert die symbolische Regression explizite mathematische Ausdrücke qsr(z|wsr) mithilfe eines genetischen Algorithmus, der sowohl Funktions- als auch Parameterräume durchsucht, um den Verlust l(wsr|θ) auf dem Trainingssatz θ unabhängig vom vorherigen neuronalen Netzwerktraining zu minimieren, wobei Der Index „sr“ weist auf eine symbolische Regression hin. Die resultierenden analytischen Ausdrücke liefern uns eine Liste von Hypothesen für quantitative Modelle der Physik, die den molekularen Aufbauprozess steuert. Zur weiteren Untersuchung werden diese Hypothesen nach einer Kombination aus statistischer Wahrscheinlichkeit (d. h. wie gut sie alle verfügbaren Daten berücksichtigen) und mathematischer Komplexität eingestuft.

Die Bildung von Ionenpaaren in Wasser ist ein paradigmatischer Aufbauprozess, der durch Vielteilchenwechselwirkungen im umgebenden Lösungsmittelmedium gesteuert wird. Obwohl MD den Prozess effizient simulieren kann, stellt die kollektive Reorganisation von Wassermolekülen bis heute eine Herausforderung für die Formulierung quantitativer mechanistischer Modelle dar18 (Abb. 1b).

Der Algorithmus lernte schnell, die Bildung von Ionenpaaren optimal abzutasten. Für Lithium- (Li+) und Chlorid- (Cl−) Ionenpaare in Wasser (Abb. 1b, c) verwendete das Netzwerk den interionischen Abstand rLiCl und 220 molekulare Merkmale x1, …, x220, die die Winkelanordnung von Wassersauerstoff- und Wasserstoffatomen beschreiben in einem bestimmten Abstand von jedem Ion19. Diese Koordinaten liefern eine allgemeine Darstellung des Systems, die in Bezug auf physikalische Symmetrien und den Austausch von Atombezeichnungen invariant ist. Nach den ersten 500 Iterationen stimmen die vorhergesagte und die beobachtete Anzahl von TPs überein (Abb. 1c). Die Probenahme ist etwa zehnmal schneller als bei herkömmlichem TPS (Extended Data Abb. 1). Wir stellen fest, dass diese Beschleunigung vollständig durch eine Verbesserung der Effizienz beim Abtasten neuer Übergangspfade und ohne Verzerrung der zugrunde liegenden Dynamik erreicht wird. Wir haben die erlernte Committor-Funktion weiter validiert, indem wir ihre Vorhersagen anhand unabhängiger Simulationen überprüft haben. Aus 1.000 Konfigurationen, die im Training nicht verwendet wurden, haben wir jeweils 500 unabhängige Simulationen initiiert und den pB des erfassten Committors als den Anteil der Trajektorien geschätzt, die zuerst in den ungebundenen Zustand eintreten. Prognostizierte und in die Stichprobe einbezogene Committoren stimmen quantitativ überein (Abb. 1d).

Entgegen einer allgemeinen Sorge bei Modellen des maschinellen Lernens ist der erlernte Mechanismus allgemein und beschreibt mit minimalen Anpassungen den Aufbau chemisch unterschiedlicher Ionenspezies. Wir führten Transferlernen auf fünf weiteren Systemen durch, indem wir Änderungen nur in der letzten linearen Schicht des trainierten Netzwerks zuließen, die ein einzelnes Neuron enthielt (Abb. 1e und Extended Data Abb. 2). Eine sehr begrenzte Anzahl neuer simulierter Übergänge reicht aus, um das Netzwerk, das den LiCl-Committor enthält, so anzupassen, dass der Committor für LiI, NaCl, NaI, CsCl und CsI korrekt vorhergesagt wird.

Wir haben auch Multiionenmodelle gebaut, die sich über den gesamten chemischen Raum erstrecken. Als Reporter für Ionengröße und Energie haben wir die Parameter Partikelgröße σ und Dispersionsenergie ϵ des Lennard-Jones-Potentials in die Merkmalsvektoren x einbezogen. Wir fanden heraus, dass Modelle, die auf der kombinierten TP-Statistik für verschiedene Ionenpaarkombinationen trainiert wurden, im chemischen Raum inter- und extrapolieren können und so vernünftige Vorhersagen für den Assoziationsmechanismus von Ionenspezies treffen, auf die sie nicht trainiert wurden (Extended Data, Abb. 3).

Lösungsmittelumlagerungen spielen eine entscheidende Rolle bei der Bestimmung der Ionenanordnung. Die Attributionsanalyse für ein Modell, das auf der Anordnung von LiCl, LiI, NaCl, NaI, CsCl und CsI trainiert wurde, identifizierte gleichzeitig den interionischen Abstand Rion und die Lennard-Jones-Parameter als entscheidend (Abb. 2a). Darüber hinaus steuern die Symmetriefunktionen, die die Geometrie der Wassermoleküle um das Kation herum beschreiben, den Aufbaumechanismus. Als wichtigstes der 176 molekularen Merkmale, die zur Beschreibung des Lösungsmittels verwendet werden, quantifiziert x7 die Sauerstoffanisotropie in einem radialen Abstand von 0,1 nm von den Kationen (Abb. 2b). Für eine erfolgreiche Ionenpaarbildung müssen diese Wassermoleküle in der inneren Schale Platz für das ankommende Anion schaffen. Die Bedeutung der Wasserumlagerung innerhalb der Schale steht im Einklang mit einer visuellen Analyse, die Atome in einem TP entsprechend ihrem Beitrag zum Committorgradienten hervorhebt (Abb. 1b).

a, Eingaberelevanz für alle 179 Eingabekoordinaten, die für Deep Learning verwendet werden. Die ersten 176 beschreiben die Geometrie von Wassermolekülen um Kationen und Anionen. Die übrigen sind der interionische Abstand rion und die Lennard-Jones-Parameter, wobei σ die Ionengröße ist. b, Schematische Darstellung der wichtigsten Lösungsmittelumorientierung. Die Symmetriefunktion x7 gibt die Geometrie der Wassersauerstoffatome (O, in Blau) bei 0,1 nm um das Kation (in Rosa) an (siehe Kasten für die Definition von x7, wobei rij und rik die Abstände zwischen dem zentralen Kation i und Sauerstoff sind). Atome j und k, und ϑijk ist der Winkel, den das zentrale Kation und zwei Sauerstoffatome bilden). c, Pareto-Diagramm aller Modelle, destilliert durch symbolische Regression. Jeder Punkt entspricht einem alternativen Modell qsr(z|wsr), gefärbt entsprechend der Anzahl der verwendeten Eingabekoordinaten (Nin). Das rote Kreuz kennzeichnet das optimale Modell am Knie der Pareto-Front. d, Multiionenmodell aus symbolischer Regression, das den Zusammenbaumechanismus von LiCl, LiI, NaCl, NaI, CsCl und CsI in Wasser beschreibt. Das Modell, \(q\left({r}_{{\mathrm{ion}}},\sigma ,{x}_{7},{x}_{15};{\sigma }_{w} \right)\) ist eine Funktion des interionischen Abstands (Rion- und Ionengröße σ in Einheiten des Wassergrößenparameters σw = 0,315 nm) und der Geometrie des Wassers um die Kationen (x7 und x15). e, Validierung des Multiionen-Assemblierungsmodells durch Kreuzkorrelation zwischen ungeschulten Proben-Committern und der Vorhersage für jede Ionenspezies separat, hier gezeigt für LiCl und CsCl (siehe erweiterte Daten Abb. 4a für alle übrigen Spezies). Der Durchschnitt der Stichproben-Committors (blaue Linie) und ihr SD (orange schattiert) wurden in den durch die vertikalen Schritte angegebenen Bins des erlernten Committors berechnet. Als Referenz zeigt die rote Linie die Identität an.

Quelldaten

Die symbolische Regression liefert ein quantitatives und interpretierbares Mehrionenmodell des Montagemechanismus. In unabhängigen symbolischen Regressionen variierten wir die Anzahl der Eingaben und den Komplexitätsnachteil. Anschließend haben wir Modelle in einem Pareto-Diagramm ausgewählt (Abb. 2c). Modelle am Knie der Pareto-Front bieten gute Kompromisse zwischen der Modellqualität, gemessen am Verlust, und der Modellkomplexität, gemessen an der Anzahl der mathematischen Operationen.

Das destillierte Multiionenmodell ist interpretierbar und bietet physikalische Einblicke in die Anordnung einwertiger Ionen in Wasser (Abb. 2d). Im führenden Term in q wird ein skalierter Ionengrößenparameter σ vom interionischen Abstand subtrahiert, was mit der physikalischen Intuition übereinstimmt. Im zweiten Term moduliert die Ionengröße nichtlinear den Deskriptor x7 der Wassergeometrie in der Nähe der Kationen (Abb. 2b). Im letzten Term berichtet x15 über eine weiter entfernte Solvatation, die nicht durch die Ionenidentität moduliert wird. Trotz seiner Einfachheit ist das reduzierte Modell für alle hier betrachteten monovalenten Ionenarten genau (Abb. 2e und erweiterte Daten Abb. 4a). Ein symbolisches Regressionsmodell, das sich auf den Aufbau von LiCl konzentriert, zeigt nur, dass wir weniger Allgemeingültigkeit gegen höhere Genauigkeit eintauschen können (Extended Data Abb. 4b – d).

Bei niedriger Temperatur und hohem Druck organisiert sich ein flüssiges Gemisch aus Wasser und Methan zu einem Gashydrat, einem eisähnlichen Feststoff20. Bei diesem Phasenübergang fügen sich Wassermoleküle zu einem komplizierten Kristallgitter mit regelmäßig beabstandeten, mit Methan gefüllten Käfigen zusammen (Abb. 3a). Trotz der kommerziellen Relevanz bei der Erdgasverarbeitung ist der Mechanismus der Gashydratbildung noch immer unvollständig verstanden, was durch den Vielteilchencharakter des Keimbildungsprozesses und die Konkurrenz zwischen verschiedenen Kristallformen erschwert wird20. Die Untersuchung des Keimbildungsmechanismus ist für Experimente eine Herausforderung und aufgrund der außerordentlichen Seltenheit der Ereignisse im MD-Gleichgewicht unmöglich.

a, Molekülkonfigurationen, die den Keimbildungsprozess veranschaulichen, extrahiert aus einer atomistischen MD-Simulation in explizitem Lösungsmittel. Der Kern bildet sich im Wasser, wächst und führt zum Clathratkristall. 51262 (blau) und 512 (rot) Wasserkäfige (Linien) enthalten entsprechend gefärbte Methanmoleküle (Kugeln). Methanmoleküle in der Nähe des wachsenden festen Kerns werden als grüne Kugeln und Wasser als graue Stäbchen dargestellt. Der Übersichtlichkeit halber ist das Grundwasser nicht dargestellt. b, Validierung des erlernten Committors. Kreuzkorrelation zwischen dem vom trainierten Netzwerk vorhergesagten Committor und dem Committor, der durch wiederholte Stichproben aus molekularen Konfigurationen erhalten wurde, auf denen der Algorithmus nicht trainiert hat (graue Linie: Identität). c, Analyse der Eingabewichtigkeit. Die drei wichtigsten Eingabekoordinaten sind mit der Temperatur T, der Anzahl der Oberflächengewässer nw und der Anzahl der 51262 Kristalle nc versehen. d: Ein durch symbolische Regression destilliertes datengesteuertes quantitatives mechanistisches Modell zeigt einen Wechsel im Keimbildungsmechanismus. In der Gleichung sind nw,0 und T0 die Referenzzahl der Oberflächenwassermoleküle bzw. die Referenztemperatur und α, β, γ und δ sind numerische Konstanten. Analytische Iso-Committor-Oberflächen für nw,0 = 2, T0 = 270 K, α = 0,0502, β = 3,17, γ = 0,109 K−1, δ = 0,0149 (von links nach rechts: gelb, pB = 1/(1 + e −4); Blau, 1/2; Grün, 1/(1 + e4)). Die strukturellen Einschübe veranschaulichen die beiden konkurrierenden Mechanismen bei niedriger und hoher Temperatur.

Quelldaten

Innerhalb weniger Stunden Rechenzeit auf einer einzelnen Grafikverarbeitungseinheit (GPU) extrahierte der Algorithmus den Keimbildungsmechanismus aus 2.225 TPs, der die Bildung von Methanclathraten zeigt, was einer simulierten Dynamik von insgesamt 445,3 μs entspricht. Die Trajektorien wurden durch umfangreiche TPS-Simulationen bei vier verschiedenen Temperaturen erstellt und stellten einen bereits vorhandenen Trainingssatz für unseren Algorithmus dar21. Wir haben molekulare Konfigurationen anhand von 22 Merkmalen beschrieben, die üblicherweise zur Beschreibung von Keimbildungsprozessen verwendet werden (Ergänzungstabelle 1). Wir haben die Temperatur T, bei der ein TP erzeugt wurde, als zusätzliches Merkmal betrachtet und das Committor-Modell auf die kumulativen Trajektorien trainiert. Wir haben gezeigt, dass der gelernte Committor als Funktion der Temperatur genau ist, indem wir seine Vorhersagen für 160 unabhängige Konfigurationen validiert haben (Abb. 3b). Generative Modelle konstruierten kürzlich Verteilungsfunktionen bei Temperaturen, die nicht erfasst wurden22. Indem wir Daten bei T = 280 K oder 285 K im Training weglassen, zeigen wir, dass der gelernte Committor zufriedenstellend interpoliert und auf nicht abgetastete thermodynamische Zustände extrapoliert (Erweiterte Daten, Abb. 5).

Die Temperatur T ist der kritischste Faktor für das Ergebnis einer Simulationstrajektorie, gefolgt von der Anzahl nw der Oberflächenwassermoleküle und der Anzahl nc von 51262 Käfigen, definiert durch das Vorhandensein von 12 Fünfecken (512) und zwei Sechsecken (62) (Abb . 3c). Alle drei Variablen spielen eine wesentliche Rolle in der klassischen Theorie der homogenen Keimbildung21. Die freie Aktivierungsenergie ΔG für die Keimbildung wird durch die Größe des wachsenden Keims bestimmt, parametrisiert durch die Menge an Oberflächenwasser und – im Falle einer kristallinen Struktur – durch die Anzahl der 51262 Käfige. Die Temperatur bestimmt über den Grad der Übersättigung die Größe des kritischen Kerns, die Höhe der Barriere für die freie Energie der Keimbildung und die Rate.

Die symbolische Regression destillierte einen mathematischen Ausdruck, der einen temperaturabhängigen Wechsel im Keimbildungsmechanismus offenbarte. Der Mechanismus wird durch q (nw, nc, T) quantifiziert (Abb. 3d und Ergänzungstabelle 1). Bei niedrigen Temperaturen entscheidet allein die Größe des Kerns über das Wachstum. Bei höheren Temperaturen gewinnt die Zahl von 51262 Wasserkäfigen an Bedeutung, was durch gekrümmte Iso-Committor-Oberflächen angezeigt wird (Abb. 3d). Dieses datengesteuert generierte mechanistische Modell zeigt den Übergang vom amorphen Wachstum bei niedrigen Temperaturen zum kristallinen Wachstum bei höheren Temperaturen21,23.

Proteine, Nukleinsäuren und Polymere können sich spontan selbst organisieren, indem sie sich in geordnete Strukturen falten. Angewandt auf den Knäuel-zu-Kristall-Übergang eines Homopolymers24,25 identifizierte der Algorithmus problemlos den bisher schwer fassbaren Mechanismus auf verschiedenen Auflösungsebenen (Extended Data Abb. 6). Bei niedriger Auflösung verwendeten wir einen ausgewählten Satz von 36 physikalischen Eigenschaften, gemittelt über das Polymer. Die Attributionsanalyse mit anschließender symbolischer Regression stellte den Committor als eine nichtlineare Funktion der Orientierungsordnung Q6 und der potentiellen Energie U allein dar, was sich als äußerst prädiktiv erwies (Abb. 4 und erweiterte Daten Abb. 6a). Bei hoher Auflösung erzeugte Deep Learning eine Committor-Funktion von vergleichbarer Qualität in einem Raum von 384 allgemeinen Deskriptoren, die die lokale Umgebung jedes Polymerkügelchens (Extended Data Abb. 6b) in Bezug auf die Anzahl der Nachbarn und den lokalen Bindungsordnungsparameter q6 darstellten und die Verbindungskoeffizienten cij, die die Korrelation zwischen den lokalen Umgebungen der Perlen i und j messen. Der Algorithmus lernte so genaue Committor-Darstellungen sowohl im Hinblick auf viele allgemeine als auch auf wenige systemspezifische Merkmale und destillierte letztere in eine kompakte und physikalisch aufschlussreiche Funktion der Orientierungsordnung und Energie.

a, Darstellung des erlernten Mechanismus. Die Wärmekarte (Farbbalken) stellt ein reduziertes explizites Modell des Committors pB = pF auf den gefalteten Zustand dar, wie durch den Ausdruck \({q}_{{\mathrm{B}}}\left(U,{Q}) reproduziert. _{6}\right)=\alpha (U-{U}_{0})+\beta \log \left({Q}_{6}-{Q}_{6,0}\right)+ \gamma\), wobei U die gesamte potentielle Energie des Polymers ist, Q6 seine Kristallinität quantifiziert und die numerischen Konstanten α = −7,144, β = 3,269, γ = 11,942, U0 = −2,351 und Q6,0 = 0,035 sind. Einschübe: Molekülkonfigurationen des Polymers bei pB = 0, 0,5 und 1. Polymerkügelchen werden entsprechend ihrem q6-Wert gefärbt, von weiß (niedrige Werte) bis dunkelorange (hohe Werte). b, Validierung des erlernten Committors. Kreuzkorrelation zwischen dem vom trainierten Netzwerk vorhergesagten Committor und dem Committor, der durch wiederholte Stichproben aus molekularen Konfigurationen erhalten wird, für die der Algorithmus nicht trainiert hat. Der Durchschnitt der Stichproben-Committors (blaue Linie) und ihr SD (orange schattiert) wurden in den durch die vertikalen Schritte angegebenen Bins des erlernten Committors berechnet. Als Referenz zeigt die rote Linie die Identität an.

Quelldaten

Membran-Protein-Komplexe spielen eine grundlegende Rolle bei der Organisation lebender Zellen. Hier untersuchten wir den Aufbau des Transmembranhomodimers des Sättigungssensors Mga2 in einer Lipiddoppelschicht in der quasi-atomistischen Martini-Darstellung (Abb. 5a)26. In umfangreichen MD-Gleichgewichtssimulationen wurde die spontane Assoziation zweier Mga2-Transmembranhelices beobachtet, es kam jedoch innerhalb von etwa 3,6 ms zu keiner Dissoziation (entspricht einer Berechnungszeit von mehr als 6 Monaten)26.

a, Schnappschüsse während eines Mga2-Dimerisierungsereignisses (von rechts nach links). Die Transmembranhelices sind als orangefarbene Flächen dargestellt, die Lipidmoleküle in Grau und Wasser in Cyan. b, Selbstkonsistenz. Kumulierte Anzahl der generierten (blaue Linie) und erwarteten (orange gestrichelte Linie) Anzahl von Übergängen. Die grüne Kurve zeigt die kumulative Differenz zwischen den beobachteten und den erwarteten Zählungen. c, Validierung des erlernten Committors. Kreuzkorrelation zwischen dem vom trainierten Netzwerk vorhergesagten Committor und dem Committor, der durch wiederholte Stichproben aus molekularen Konfigurationen erhalten wurde, für die das Committor-Modell nicht trainiert wurde. Der Durchschnitt der Stichproben-Committors (blaue Linie) und ihr SD (orange schattiert) wurden in den durch die vertikalen Schritte angegebenen Bins des erlernten Committors berechnet. Als Referenz zeigt die rote Linie die Identität an. d, Schematische Darstellung der beiden wichtigsten Koordinaten, der interhelikalen Kontakte an den Positionen 9 und 22. e, Hierarchische Clusterung aller TPs. Dendrogramm als Funktion der TP-Ähnlichkeit (dynamische Zeitverzerrung, siehe „Mga2-Transmembran-Dimer-Anordnung in Lipidmembran“ in Methoden), berechnet in der durch die Kontakte 9 und 22 definierten Ebene (zwei Hauptcluster: blau, orange). f,g, Pfaddichte (graue Schattierung) für die beiden Hauptcluster in e, berechnet in der durch die Kontakte 9 und 22 definierten Ebene. Für jeden Cluster wird ein repräsentatives TP von ungebunden (gelb) bis gebunden (rot) angezeigt. Die Isolinien des Committors, wie durch die symbolische Regression \({q}_{{\mathrm{B}}}({x}_{9},{x}_{22})=-\exp ({ x}_{9}^{2})\log ({x}_{9}-\frac{{x}_{9}}{\log ({x}_{22})})\), werden als beschriftete durchgezogene Linien dargestellt.

Quelldaten

Die Pfadabtastung ist von Natur aus parallelisierbar, was es uns ermöglichte, in 20 Tagen fast 4.000 Dissoziationsereignisse auf einem parallelen Supercomputer abzutasten (Abb. 5b). Die Zeitintegration von MD-Trajektorien verursacht den höchsten Rechenaufwand und ist nur bedingt parallelisierbar. Eine einzelne Instanz des Algorithmus kann jedoch gleichzeitig virtuelle Experimente an einer beliebigen Anzahl von Kopien des physischen Modells orchestrieren (durch Steuerung paralleler Markov-Ketten-Monte-Carlo-(MC)-Stichprobenprozesse) und durch Training der kumulativen Ergebnisse aus allen lernen .

Wir haben molekulare Konfigurationen mithilfe von Kontakten zwischen entsprechenden Resten entlang der beiden Helices dargestellt und als Referenz eine Reihe von handgefertigten Merkmalen hinzugefügt, die die Organisation von Lipiden um die Proteine ​​herum beschreiben (Erweiterte Daten, Abb. 7 und Ergänzungstabelle 2). Wir validierten das Modell anhand von Committor-Daten für 548 Molekülkonfigurationen, die nicht im Training verwendet wurden, und stellten fest, dass die Vorhersagen über den gesamten Übergangsbereich zwischen gebundenen und ungebundenen Zuständen genau waren (Abb. 5c).

In einer bemerkenswerten Reduzierung der Dimensionalität erreichte die symbolische Regression eine genaue Darstellung des erlernten Committors als einfache Funktion von nur zwei Aminosäurekontakten (Abb. 5d und erweiterte Daten, Abb. 8). Die symbolische Regression liefert uns eine Liste von Hypothesen für quantitative Modelle der Physik, die den Molekülassemblierungsprozess regeln (Ergänzungstabelle 3). Diese Hypothesen werden nach einer Kombination aus statistischer Wahrscheinlichkeit (d. h. wie gut sie alle verfügbaren Daten berücksichtigen) und ihrer mathematischen Komplexität eingestuft. Unter den Ausdrücken am Ende des Pareto-Diagramms, das heißt mit vergleichbarer Vorhersagekraft und Komplexität, haben wir denjenigen ausgewählt, der eine klare Interpretation im Sinne zweier konkurrierender Montagemechanismen bietet, die durch die Bildung von Kontakten beschrieben werden, die an dem einen oder anderen Helix-Terminus beginnen .

Wir haben alle abgetasteten TPs auf die durch diese beiden Kontakte definierte Ebene projiziert, die Abstände zwischen ihnen berechnet und eine hierarchische Trajektorienclusterung durchgeführt (Abb. 5e). TPs organisieren sich in zwei Hauptclustern, die zwei konkurrierende Aufbauwege offenbaren, wobei der anfängliche Helixkontakt oben (Abb. 5f) oder unten (Abb. 5g) liegt. Unerwarteterweise27 sagt allein die Helix-Dimer-Geometrie den Fortschritt des Zusammenbaus voraus, was bedeutet, dass das Lipid-„Lösungsmittel“ im Gegensatz zum Wasserlösungsmittel bei der Ionenpaarbildung implizit in den interhelikalen paarweisen Kontakten kodiert ist18. Wie bei der Polymerfaltung und Ionenbindung erwies sich somit ein ausreichend großer Raum allgemeiner geometrischer Merkmale als ausreichend für die Konstruktion vollständig prädiktiver Committoren durch Deep Learning. Dieser Befund steht im Einklang mit der Einbettungstheorie28 und impliziert, dass die Verwendung einer kleinen, aber ausreichenden Anzahl allgemeiner Merkmale genauso effektiv ist wie kollektive Variablen, die auf physikalischer und chemischer Intuition basieren.

Die maschinengesteuerte Trajektorienerfassung ist allgemein und kann sofort angepasst werden, um Vielkörperdynamiken mit einer Vorstellung vom „wahrscheinlichen Schicksal“ ähnlich dem Committor abzutasten. Dieses grundlegende Konzept der statistischen Mechanik reicht vom Schachspiel10 über die Proteinfaltung3,7 bis zur Klimamodellierung29. Die Simulationsmaschine – in unserem Fall MD – wird als Blackbox behandelt und kann durch andere dynamische Prozesse ersetzt werden, ob reversibel oder nicht. Sowohl das statistische Modell, das die Verlustfunktion definiert, als auch die Technologie des maschinellen Lernens können auf spezifische Probleme zugeschnitten werden. Anspruchsvollere Modelle werden in der Lage sein, aus weniger Daten mehr zu lernen oder experimentelle Einschränkungen zu berücksichtigen. Einfachere Regressionsschemata5 können neuronale Netze15 ersetzen, wenn die Kosten für die Abtastung von Trajektorien das Volumen der Trainingsdaten stark einschränken.

Die Definition der Grenzen der metastabilen Zustände, wie sie unsere Methode erfordert, kann nicht trivial sein. Um die Zustandsdefinitionen zu verfeinern, kann das Committor-Framework auf iterative Weise verwendet werden. Man beginnt mit einer sehr konservativen Methode, erhält eine grobe Schätzung des Committors und erweitert dann die Zustände auf Werte des Committors nahe 0 und 1. A verwandt Das Problem ist das Vorhandensein nicht zugewiesener metastabiler Zustände, die zu langen Übergangspfaden führen. Hier können maschinelle Lernansätze, die auf die Zustandserkennung abzielen, eine Lösung bieten.

Da unsere Methode keine unphysikalischen Kräfte zur Beschleunigung der Abtastung nutzt, ist ihre Anwendbarkeit auf Prozesse beschränkt, für die die Übergangspfade kurz genug sind, damit die zugrunde liegende Simulationsmaschine eine angemessene Anzahl davon abtasten kann. Diese Einschränkung ist nicht dramatisch, da die Dauer der TPs nur schwach von der Barrierenhöhe abhängt. Bei Bedarf kann unsere Methode auch mit voreingenommenen Potentialflächen eingesetzt werden, um die Dynamik zu beschleunigen.

Die experimentelle Validierung des Aufbaumechanismus kann durch direkte Beobachtung, beispielsweise durch Einzelmolekülexperimente, oder durch Störung, beispielsweise durch Änderung des thermodynamischen Zustands oder der chemischen Zusammensetzung, erfolgen. Im ersteren Fall würden die Simulationen Vorhersagen über Observable entlang der Montagerouten treffen. Für Letzteres müsste die Auswirkung der Störung auf den Montageprozess in den Simulationen rekapituliert werden. Beispielsweise veränderten Temperaturänderungen die Geschwindigkeit und den Weg der Clathrat-Keimbildung21.

Die maschinengesteuerte Mechanismuserkennung30 integriert problemlos Fortschritte im maschinellen Lernen, die auf Kraftfelder19,31, Probenahme32,33,34 und molekulare Darstellung19,31,35 angewendet werden. Zunehmende Rechenleistung und Fortschritte in der symbolischen künstlichen Intelligenz werden es Algorithmen ermöglichen, immer genauere mathematische Beschreibungen der komplexen Prozesse zu destillieren, die in hochdimensionalen Daten verborgen sind36. Wie hier gezeigt, können maschinengestützte Probenahme und Modellvalidierung in Kombination mit symbolischer Regression den wissenschaftlichen Entdeckungsprozess unterstützen.

Der Committor pB(x) ist die Wahrscheinlichkeit, dass eine bei Konfiguration Wir erwarten, dass nA- und nB-Trajektorien mit der Binomialwahrscheinlichkeit \(p({n}_{{\mathrm{A}}},{n}_{{\mathrm{B}}) in A und B enden. }| {{{\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}_{{\ mathrm{B}}}}\). Für k Schusspunkte xi definiert die kombinierte Wahrscheinlichkeit die Wahrscheinlichkeit \({{{\mathcal{L}}}}=\mathop{\prod }\nolimits_{i = 1}^{k}p({n}_{ {\mathrm{A}}}(i),{n}_{{\mathrm{B}}}(i)| {{{{\bf{x}}}}}_{i})\). Hier ignorieren wir die Korrelationen, die bei schnellen trägheitsdominierten Übergängen für Flugbahnen auftreten, die mit entgegengesetzten Anfangsgeschwindigkeiten abgeschossen werden11,18. Wir modellieren den unbekannten Committor mit einer parametrischen Funktion und schätzen seine Parameter w durch Maximieren der Wahrscheinlichkeit \({{{\mathcal{L}}}}\) (Lit. 5,15). Wir stellen sicher, dass 0 ≤ pB(x) ≤ 1, indem wir den Committor in Form einer sigmoidalen Aktivierungsfunktion schreiben, \({p}_{{\mathrm{B}}}[q({{{\bf{x}} }}| {{{\bf{w}}}})]=1/(1+\exp [-q({{\bf{x}}}}| {{{\bf{w}}} })])\). Hier modellieren wir die Log-Wahrscheinlichkeit q(x|w) mithilfe eines neuronalen Netzwerks15 und stellen die Konfiguration mit einem Vektor x von Merkmalen dar. Für N > 2 Zustände S stellt die Multinomialverteilung ein Modell für p(n1, n2, . . . , nN|x) bereit, und das Schreiben der Committoren in Zustände S im Sinne der Softmax-Aktivierungsfunktion gewährleistet die Normalisierung, \(\mathop {\sum }\nolimits_{S = 1}^{N}{p}_{S}=1\). Die im Training verwendete Verlustfunktion l(w|θ) ist der negative Logarithmus der Wahrscheinlichkeit \({{{\mathcal{L}}}}\).

TPS4,12 ist eine leistungsstarke Markov-Ketten-MC-Methode im Pfadraum zum Abtasten von TPs. Der Zwei-Wege-Schusszug ist ein effizienter Vorschlagzug in TPS4. Es besteht aus der zufälligen Auswahl eines Schusspunkts XSP auf dem aktuellen TP χ gemäß der Wahrscheinlichkeit psel(XSP|χ), dem Zeichnen zufälliger Maxwell-Boltzmann-Geschwindigkeiten und der Ausbreitung zweier Versuchstrajektorien von XSP, bis sie einen der Zustände erreichen. Da sich eine der Versuchsbahnen ausbreitet, nachdem zunächst alle Impulse am Startpunkt umgekehrt wurden, also in der Zeit rückwärts, kann ein kontinuierlicher TP konstruiert werden, wenn beide Versuche in unterschiedlichen Zuständen enden. Bei gegebenem TP χ wird ein neues TP χ′, das durch Zwei-Wege-Schießen erzeugt wird, mit der Wahrscheinlichkeit37\({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 ))\). Wenn der neue Pfad abgelehnt wird, wird χ wiederholt.

Wenn man den Committor kennt, ist es möglich, die Rate, mit der TPs erzeugt werden, zu erhöhen, indem man die Auswahl der Aufnahmepunkte auf das Übergangszustandsensemble37 ausrichtet, d. h. auf Regionen mit hoher reaktiver Wahrscheinlichkeit p(TP|X). Für den Zwei-Zustands-Fall entspricht dies einer Ausrichtung auf die pB(x) = 1/2-Isofläche, die die Übergangszustände mit q(x) = 0 definiert. Konstruieren eines Algorithmus, der neue Aufnahmepunkte auswählt, die auf die aktuell beste Schätzung ausgerichtet sind Für das Übergangszustandsensemble, das iterativ lernt, seine Schätzung basierend auf jedem neu beobachteten Schießergebnis zu verbessern, müssen wir Erkundung und Ausbeutung in Einklang bringen. Zu diesem Zweck wählen wir den neuen Aufnahmepunkt \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})],\), wobei größere Werte von γ zu einer Zunahme der Erkundung führen. Die Lorentz-Verteilung bietet einen Kompromiss zwischen Produktionseffizienz und der gelegentlichen Erkundung außerhalb des Übergangszustands, die zur Erprobung alternativer Reaktionskanäle erforderlich ist.

Mit der erlernten Committor-Funktion kann man die Definition der Staatsgrenzen optimieren. Eine anfänglich enge Zustandsdefinition kann aufgeweicht werden, indem die Grenzen nach außen verschoben werden, beispielsweise auf pB(x) = 0,1 und pB(x) = 0,9. Diese Lockerung führt zu kürzeren TPs und beschleunigt die Probenahme.

Die Beziehung zwischen dem Committor und der Übergangswahrscheinlichkeit11 ermöglicht es uns, die erwartete Anzahl von TPs zu berechnen, die durch das Schießen aus einer Konfiguration es mit dem beobachteten Schießergebnis. Die erwartete Anzahl von Übergängen \({n}_{{{{\rm{TP}}}}}^{\exp }\), berechnet über ein Fenster, das die k letzten Zwei-Wege-Schießversuche4 enthält, beträgt \({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)\), wobei pB(xi, i) die Committor-Schätzung für den Probeschießpunkt Xi in Schritt i vor der Beobachtung des Schießergebnisses ist. Wir beginnen mit dem Lernen, wenn das vorhergesagte (\({n}_{{{{\rm{TP}}}}}^{\exp }\)) und das tatsächlich generierte (\({n}_{{{{\rm {TP}}}}}^{{{\rm{gen}}}}}\)) Anzahl der TPs unterschiedlich. Wir definieren einen Effizienzfaktor, \({\alpha }_{{{{\rm{eff}}}}}=\min (1,{(1-{n}_{{{{\rm{TP}} }}}^{{\mathrm{gen}} }/{n}_{{{{\rm{TP}}}}}^{{\mathrm{exp}} })}^{2})\) , wobei ein Wert von Null eine perfekte Vorhersage anzeigt (Extended Data Abb. 9). Indem wir nur bei Bedarf trainieren, vermeiden wir eine Überanpassung. Hier haben wir αeff verwendet, um die Lernrate im Gradientenabstiegsalgorithmus zu skalieren. Darüber hinaus findet kein Training statt, wenn αeff unter einem bestimmten Schwellenwert liegt (weiter unten für jedes System angegeben).

Molekulare Mechanismen können mit unterschiedlichen Auflösungsstufen beschrieben werden. Man kann viele hochauflösende Features verwenden, die lokale Eigenschaften quantifizieren, oder weniger niedrig aufgelöste Features, die globale Eigenschaften messen. Während Funktionen mit hoher Auflösung in der Regel leicht verfügbar sind, hängt die Auswahl sinnvoller Funktionen mit niedriger Auflösung vom physikalischen Verständnis ab. Mit einem Fokus auf seltene molekulare Ereignisse wollten wir die Merkmale in einer Auflösungshierarchie anordnen, die von kartesischen Koordinaten atomarer Positionen – der höchstmöglichen Auflösung – bis hin zu einer einzigen Größe, dem Committor, reicht.

Wir haben die neuronalen Netze in dieser Studie entworfen, um sie zum Erlernen der Auflösungshierarchie von Features zu ermutigen. Neuronale Netze haben gezeigt, dass sie in der Lage sind, Merkmale mit niedriger Auflösung von Merkmalen mit hoher Auflösung zu lernen, beispielsweise bei der Verwendung in der Bilderkennung. Aus praktischer Sicht nimmt die Schichtbreite unserer Netzwerke ständig ab und bewegt sich von der Eingabe zur Ausgabe. Darüber hinaus verringern wir die Dropout-Wahrscheinlichkeit, wenn wir uns im Netzwerk nach oben bewegen, da die erlernten Funktionen immer globaler (und daher weniger redundant) werden, während sie in tiefere Schichten vordringen. Dies spiegelt sich auch in der unterschiedlichen Architektur der Clathratbildung wider, bei der wir aufgrund der bereits grobkörnigen und systemspezifischen Merkmale ein vergleichsweise einfaches Pyramiden-Feed-Forward-Netzwerk verwendeten.

Wir gehen davon aus, dass bei jedem spezifischen molekularen Prozess nur wenige der vielen Freiheitsgrade tatsächlich die Übergangsdynamik steuern. Wir identifizieren die Eingaben in das Committor-Modell, die nach dem Training die größte Rolle bei der Bestimmung seiner Ausgabe spielen. Zu diesem Zweck berechnen wir zunächst einen Referenzverlust, lref = l(w, θ), über den ungestörten Trainingssatz, um ihn mit den Werten zu vergleichen, die wir erhalten, indem wir jeden Eingang einzeln stören38. Anschließend mitteln wir den Verlust \(l({{\bf{w}}}},{\widetilde{{{{\bf{\uptheta }}}}}}_{i})\) über ≥100 Störungen Trainingsmengen \({\widetilde{{{{\bf{\uptheta }}}}}}_{i}\) mit zufällig permutierten Werten der Eingabekoordinate i in der Batch-Dimension. Die durchschnittliche Verlustdifferenz \({{\Delta }}{l}_{i}=\left\langle l({{{\bf{w}}}},{\widetilde{{{{\bf{\uptheta }}}}}}_{i})\right\rangle -{l}_{{{{\rm{ref}}}}}\) ist groß, wenn die i-te Eingabe die Ausgabe des trainierten Modells stark beeinflusst, das heißt, es ist relevant für die Vorhersage des Committors.

In der niedrigdimensionalen Teilmenge, die nur aus den relevantesten Eingaben z besteht (diejenigen mit dem höchsten Δli), generiert die symbolische Regression kompakte mathematische Ausdrücke, die sich dem vollständigen Committor annähern. Unsere Implementierung der symbolischen Regression basiert auf dem Python-Paket dcgpy39 und verwendet einen genetischen Algorithmus mit einer (N + 1)-Evolutionsstrategie. In jeder Generation werden N neue Ausdrücke durch zufällige Änderungen an der mathematischen Struktur des geeignetsten Ausdrucks der übergeordneten Generation generiert. Anschließend wird eine gradientenbasierte Optimierung verwendet, um für jeden Ausdruck die besten Parameter zu finden. Der geeignetste Ausdruck wird dann als übergeordneter Ausdruck für die nächste Generation ausgewählt. Die Fitness jedes Versuchsausdrucks pB(z) wird durch \({l}_{{{{\rm{sr}}}}}({{{{\bf{w}}}}}_{{{ {\rm{sr}}}}}| {{{\bf{\uptheta }}}})\equiv -\log {{{\mathcal{L}}}}[{p}_{{\mathrm{ B}}}({{{{\bf{z}}}}}_{{{\rm{sp}}}}})]+\lambda C\), wobei wir den Regularisierungsterm λC zum hinzugefügt haben Log-Likelihood (siehe „Maximum-Likelihood-Schätzung der Committor-Funktion“), um Ausdrücke einfach zu halten und eine Überanpassung zu vermeiden. Hier ist λ > 0 und C ist ein Maß für die Komplexität des Versuchsausdrucks, in unserem Fall geschätzt durch die Anzahl der mathematischen Operationen.

Die symbolische Regression führt je nach Regularisierungsstärke zu Ausdrücken unterschiedlicher Komplexität. Wir wählen Ausdrücke mit einem angemessenen Kompromiss zwischen Einfachheit und Genauigkeit mithilfe eines Pareto-Diagramms (Abb. 2c) aus, in dem wir die Komplexität (gemessen als Anzahl der mathematischen Operationen) gegen die Genauigkeit (gemessen als Verlust der Validierungsdaten) auftragen. . Durch Scannen eines Bereichs von λ-Werten können wir Modelle an der Pareto-Front für die weitere Analyse identifizieren.

Wir haben die Bildung einwertiger Ionenpaare in Wasser untersucht, um die Fähigkeit des Algorithmus zu beurteilen, den Committor für Übergänge, die stark von den Freiheitsgraden des Lösungsmittels beeinflusst werden, genau zu lernen. Wir verwendeten sechs verschiedene Systemaufbauten (LiCl, LiI, NaCl, NaI, CsCl und CsI), die jeweils aus einem Kation und einem Anion in Wasser bestanden.

Alle MD-Simulationen wurden in kubischen Simulationsboxen unter Verwendung des Joung- und Cheatham-Kraftfelds40 zusammen mit TIP3P41-Wasser durchgeführt. Jede Simulationsbox enthielt ein einzelnes Ionenpaar, solvatisiert mit 370 TIP3P-Wassermolekülen. Wir haben die openMM MD-Engine42 verwendet, um die Bewegungsgleichungen in Zeitschritten von Δt = 2 fs mit einem Geschwindigkeits-Verlet-Integrator mit Geschwindigkeits-Randomisierung43 aus dem Python-Paket openmmtools fortzupflanzen. Nach einer anfänglichen NPT-Äquilibrierung bei konstantem Druck P = 1 bar und konstanter Temperatur T = 300 K wurden alle Produktionssimulationen im NVT-Ensemble bei einem konstanten Volumen V und einer konstanten Temperatur von T = 300 K durchgeführt. Die Reibung wurde auf 1 gesetzt ps−1. Nichtgebundene Wechselwirkungen wurden unter Verwendung eines Partikelnetz-Ewald-Schemas44 mit einem Realraum-Grenzwert von 1 nm und einer Fehlertoleranz von 0,0005 berechnet. Jede TPS-Simulation (bestehend aus MD-Simulationen und neuronalem Netzwerktraining) wurde auf einem halben Knoten mit einer Xeon Gold 6248-Zentraleinheit (CPU) in Verbindung mit einer RTX5000-GPU durchgeführt. In TPS wurden die assoziierten und dissoziierten Zustände anhand der interionischen Abstände definiert (die Werte für jede Ionenart finden Sie in der Ergänzungstabelle 4).

Der Committor einer Konfiguration ist invariant unter globalen Translationen und Rotationen ohne externe Felder und darüber hinaus invariant in Bezug auf Permutationen identischer Teilchen. Wir haben uns daher entschieden, die Systemkoordinaten aus dem kartesischen Raum in eine Darstellung umzuwandeln, die die physikalischen Symmetrien des Committors berücksichtigt. Um eine nahezu verlustfreie Transformation zu erreichen, verwendeten wir den interionischen Abstand zur Beschreibung des gelösten Stoffes und passten Symmetriefunktionen zur Beschreibung der Lösungsmittelkonfiguration an45. Symmetriefunktionen wurden ursprünglich entwickelt, um molekulare Strukturen in neuronalen Netzwerkpotentialen zu beschreiben19,46, wurden aber auch erfolgreich zur Erkennung und Beschreibung verschiedener Eisphasen in atomistischen Simulationen eingesetzt47. Diese Funktionen beschreiben die Umgebung eines Zentralatoms, indem sie alle identischen Teilchen in einem bestimmten radialen Abstand summieren. Die Symmetriefunktion vom Typ \({G}_{i}^{2}\) quantifiziert die Dichte von Lösungsmittelmolekülen um ein gelöstes Atom i in einer bei rs zentrierten Schale

wobei die Summe über alle Lösungsmittelatome j einer bestimmten Atomart läuft, rij der Abstand zwischen Zentralatom i und Atom j ist und η die Breite der Schale steuert. Die Funktion fc(r) ist ein Fermi-Grenzwert, definiert als:

was dafür sorgt, dass der Beitrag entfernter Lösungsmittelatome verschwindet. Der Skalarparameter αc steuert die Steilheit des Cut-Offs. Die Symmetriefunktion vom Typ \({G}_{i}^{5}\) untersucht zusätzlich die Winkelverteilung des Lösungsmittels um das Zentralatom i

wobei die Summe über alle unterschiedlichen Lösungsmittelatompaare läuft, ϑijk der Winkel ist, der zwischen den beiden Lösungsmittelatomen und dem zentralen gelösten Atom aufgespannt wird, der Parameter ζ eine gerade Zahl ist, die die Schärfe der Winkelverteilung steuert, und λ = ±1 legt fest Ort des Minimums in Bezug auf ϑijk bei π bzw. 0. Die verwendeten Parameterkombinationen finden Sie in der Ergänzungstabelle 5. Wir haben alle Eingaben so skaliert, dass sie ungefähr im Bereich [0, 1] liegen, um die numerische Stabilität des Trainings zu erhöhen. Insbesondere haben wir die Symmetriefunktionen normalisiert, indem wir sie durch die erwartete durchschnittliche Anzahl von Atomen (oder Atompaaren) für eine isotrope Verteilung im Untersuchungsvolumen dividiert haben.

Die Symmetriefunktionen vom Typ G2 zählen die Anzahl der Lösungsmittelatome im Untersuchungsvolumen; die Normalisierungskonstante \({\langle N[{G}_{i}^{2}]\rangle }_{{{{\rm{iso}}}}}\) ist daher die erwartete Anzahl von Atomen in der Sondierungsvolumen \({V}_{{{{\rm{probe}}}}}^{(2)}\)

wobei ρN die durchschnittliche Zahlendichte des untersuchten Lösungsmittelatomtyps ist. Das genaue Sondierungsvolumen für den G2-Typ kann angenähert werden als:

für kleines η und rcut > rs.

Die Funktionen vom Typ G5 enthalten einen zusätzlichen Winkelterm und zählen die Anzahl der Lösungsmittelatompaare, die sich auf gegenüberliegenden Seiten des zentralen gelösten Atoms befinden. Die erwartete Anzahl von Paaren \({\langle {N}_{{{{\rm{pairs}}}}}\rangle }_{{{{\rm{iso}}}}}\) kann berechnet werden aus die erwartete Anzahl von Atomen im untersuchten Volumen \({\langle {N}_{{{{\rm{Atome}}}}}\rangle }_{{{{\rm{iso}}}}}\) als \({\langle {N}_{{{{\rm{Atome}}}}}\rangle }_{{{\rm{iso}}}}}({\langle {N}_{{ {{\rm{Atome}}}}}\rangle }_{{{\rm{iso}}}}}-1)/2\). Dieser Ausdruck ist nur für ganzzahlige Werte von \({\langle {N}_{{{{\rm{Atome}}}}}\rangle }_{{{{\rm{iso}}}}}\) exakt. und kann sogar negativ werden, wenn \({\langle {N}_{{{{\rm{Atome}}}}}\rangle }_{{{{\rm{iso}}}}} < 1\). Wir haben daher eine Näherung verwendet, die garantiert nicht negativ ist

Die erwartete Anzahl der Atome \({\langle {N}_{{{{\rm{Atome}}}}}\rangle }_{{{{\rm{iso}}}}}\) kann berechnet werden aus das Volumen, das nach einem festen gelösten Atom und einem festen Lösungsmittelatom untersucht wird

Mit der Erwartung, dass die meisten Freiheitsgrade des Systems den Übergang nicht steuern, haben wir neuronale Netze entworfen, die nach und nach irrelevante Eingaben herausfiltern und aus den verbleibenden eine hochgradig nichtlineare Funktion erstellen. Wir haben drei verschiedene Pyramidenarchitekturen neuronaler Netzwerke getestet: „ResNet I“, „ResNet II“ und „SNN“, wobei Namen, die „ResNet“ enthalten, auf die Verwendung von Resteinheiten hinweisen48,49 und „SNN“ eine selbstnormalisierende neuronale Netzwerkarchitektur50 ist (siehe Ergänzungstabellen). 6–8 für die genauen verwendeten Architekturen). Die Architektur mit der besten Leistung ist „ResNet I“ (siehe Zusatzdatendatei 1 für einen Leistungsvergleich der verschiedenen Architekturen für alle Ionensysteme). ResNet Ich habe einen Pyramidenstapel aus fünf Resteinheiten vor der Aktivierung mit jeweils vier verborgenen Schichten verwendet. Die Anzahl der verborgenen Einheiten pro Schicht wird nach jedem verbleibenden Einheitsblock um einen konstanten Faktor \(f={(10/{n}_{{\mathrm{in}}})}^{1/4}\) reduziert sinkt von nin = 221 in der ersten Einheit auf 10 in der letzten. Zusätzlich wird nach jedem Restblock ein Dropout von 0,1fi angewendet, wobei i der Resteinheitsindex im Bereich von 0 bis 4 ist. Die Optimierung der Netzwerkgewichte erfolgt mithilfe des Adam-Gradientenabstiegs51. Für alle Architekturen wurde das Training nach jedem dritten TPS MC-Schritt für eine Epoche mit einer Lernrate von lr = αeff10−3 durchgeführt, wenn lr ≥ 10−4. Der erwartete Effizienzfaktor αeff wurde über ein Fenster von k = 100 TPS-Schritten berechnet. Wir haben das gesamte Deep Learning mit individuell geschriebenem Code auf Basis von Keras52 durchgeführt. Die TPS-Simulationen wurden mit einer angepassten Version von openpathsampling53,54 zusammen mit unserem eigenen Python-Modul durchgeführt.

Für das Transfertraining wurde die letzte Schicht mit einem einzelnen Neuron (d. h. der Log-Prädiktor) eines ursprünglich auf LiCl trainierten Modells randomisiert und alle anderen Gewichte wurden während des anschließenden Trainings auf den Aufnahmedaten für die anderen Ionenspezies unverändert gehalten ( LiI, NaCl, NaI, CsCl und CsI). Das Training wurde mit dem Adam-Optimierer mit einer Lernrate von lr = 2,5 × 10−5 durchgeführt. Der Testverlust wurde nach jeder Epoche für 20 % der als Testsatz verwendeten Daten berechnet. Das Training wurde abgebrochen, wenn über mehr als 1.000 Epochen kein Rückgang des Testverlusts beobachtet wurde. Anschließend wurde das Modell mit dem besten Testverlust verwendet.

Für die Extrapolation im chemischen Raum (Extended Data Abb. 3) haben wir ein multiionisches neuronales Netzwerk der Architektur „ResNet I“ aufgebaut. Das Modell wurde wie angegeben auf die Schießergebnisse für verschiedene Paare ionischer Spezies gleichzeitig trainiert. Es verwendete die Koordinaten aus dem Satz „SF shortranged“ zusammen mit den Lennard-Jones-Parametern ϵ und σ des Kraftfelds, um die verschiedenen Ionenarten zu unterscheiden. Das Training wurde mit dem Adam-Optimierer (lr = 10−3) durchgeführt, wobei 10 % der Daten als Testsatz verwendet wurden. Das Training wurde abgebrochen, wenn der Testverlust 1.000 Epochen lang nicht abnahm, und dann wurde das Modell mit dem besten Testverlust verwendet.

Wir haben die sieben relevantesten Koordinaten, die vom neuronalen Netzwerk mit mehreren Ionen identifiziert wurden, als Eingaben für die symbolischen Regressionen mit mehreren Ionen ausgewählt (Ergänzungstabellen 9–12). Wir haben zwischen 3 und 7 dieser relevantesten Koordinaten für unabhängige symbolische Regressionsläufe unter Verwendung der Regularisierungswerte λ = 0,001, λ = 0,0001 und λ = 0,00001 verwendet. Anschließend haben wir den in Abb. 2d dargestellten Ausdruck mithilfe des Pareto-Diagramms in Abb. 2c ausgewählt.

Wir haben auch die fünf relevantesten Koordinaten ausgewählt, die aus einem auf LiCl trainierten neuronalen Netzwerk für symbolische Regressionsläufe identifiziert wurden (Extended Data Abb. 4b – d). Wir haben die erzeugten Ausdrücke reguliert, indem wir die Gesamtzahl der elementaren mathematischen Operationen mit λ = 10−6 und λ = 10−7 bestraft haben.

Die Beiträge jedes Atoms zum Committor in einem bestimmten System X (Abb. 1b) wurden als Größe des Gradienten der Reaktionskoordinate q(x) in Bezug auf ihre kartesischen Koordinaten berechnet. Alle Gradientengrößen wurden mit der inversen Atommasse skaliert.

Der Lernalgorithmus wurde auf einen vorhandenen TPS-Datensatz der Methan-Clathrat-Keimbildung angewendet, der ursprünglich für Referenz erstellt wurde. 21. Es enthält Daten für Simulationen, die bei vier verschiedenen Temperaturen T = 270 K, 275 K, 280 K und 285 K durchgeführt wurden (Einzelheiten siehe Ergänzungstabelle 14). Neue Simulationen wurden durchgeführt, um die in der Validierung verwendeten Stichproben-Committor-Werte zu erhalten. Alle Committor-Simulationen wurden mit OpenMM 7.1.142 auf NVIDIA GeForce GTX TITAN 1080Ti-GPUs durchgeführt, wobei zwischen 6 und 18 Flugbahnen pro Konfiguration unter Verwendung des gleichen Simulationsprotokolls wie in Ref. aufgenommen wurden. 21.

Wir haben 22 verschiedene Merkmale verwendet, um Größe, Kristallinität, Struktur und Zusammensetzung des wachsenden Methanhydrat-Kristallkerns zu beschreiben (Ergänzungstabelle 1). Zusätzlich zu den Merkmalen, die molekulare Konfigurationen beschreiben, verwendeten wir die Temperatur als Eingabe für die neuronalen Netze und die symbolische Regression. In einem Pyramiden-Feed-Forward-Netzwerk mit 9 Schichten haben wir die Anzahl der Einheiten pro Schicht von 23 am Eingang auf eine in der letzten Schicht reduziert (Ergänzungstabelle 15). Das Netzwerk wurde mit dem Adam-Optimierer mit der Lernrate lr = 10−3 auf den vorhandenen TPS-Daten für alle Temperaturen trainiert, wobei 10 % der Aufnahmepunkte als Testdaten weggelassen wurden. Wir stoppten das Training, nachdem der Verlust des Testsatzes 10.000 Epochen lang nicht abgenommen hatte, und verwendeten das Modell mit dem besten Testverlust. Das gesamte neuronale Netzwerktraining wurde auf einer RTX6000-GPU durchgeführt. Wir verwendeten die drei relevantesten Koordinaten als Eingaben für symbolische Regressionsläufe mit einem Abzug auf die Gesamtzahl der elementaren mathematischen Operationen unter Verwendung von λ = 10−5.

Wir haben unseren Algorithmus für maschinelles Lernen auf vorhandene Aufnahmedaten der Polymerkristallisation angewendet24,25. Wir haben zwei verschiedene Merkmalssätze verwendet, um den Übergang zu beschreiben: einen Satz von 35 niedrig aufgelösten (grobkörnigen) Merkmalen, der auch in früheren Arbeiten verwendet wurde, und einen Satz hochauflösender Merkmale, die jede Polymerperle einzeln beschreiben. Die niedrig aufgelösten Merkmale enthalten eine Reihe globaler Maße wie die potentielle Energie U und die Steinhardt-Bindungsordnungsparameter Q4 und Q6, Beschreibungen der lokalen Umgebung ausgewählter Polymerpartikel, verschiedene Maße, die die Struktur des Polymers durch Zählen von Ketten beschreiben und Schleifen und einige ausgewählte Entfernungen (eine vollständige Liste finden Sie in der Ergänzungstabelle 16). Der hochauflösende Merkmalssatz besteht aus der Anzahl der Verbindungen, Nachbarn und den nachbarschaftsgemittelten Lechner-Dellago-Steinhardt-Bindungsordnungsparametern55 für jede Polymerperle, d. h. jede Konfiguration entspricht einem Merkmalsvektor mit 3 × 128 = 384 Einträgen.

Sowohl für die Beschreibung mit hoher als auch mit niedriger Auflösung verwendeten wir pyramidenförmige neuronale Netze (Ergänzungstabellen 17 und 18). In beiden Fällen wurde das Training mit der Adam-Gradientenabstiegsmethode mit einer Lernrate lr = 10−3 durchgeführt, wobei 20 % der Daten als Testdaten verwendet wurden. Die Modelle wurden gespeichert und der Testverlust nach jeder Epoche berechnet. Das Training wurde abgebrochen, wenn der Testverlust 10.000 Epochen lang nicht abnahm. Das Modell mit dem geringsten Testverlust wurde dann als endgültiges trainiertes Modell verwendet. Das gesamte neuronale Netzwerktraining wurde auf einer RTX6000-GPU durchgeführt.

Wir haben zwischen zwei und fünf der fünf relevantesten Merkmale mit niedriger Auflösung als Eingaben in symbolischen Regressionsläufen verwendet (Ergänzungstabellen 19–22). Wir haben reguliert, indem wir die Anzahl der elementaren mathematischen Operationen mit λ = 10−2, 10−3, 10−4, 10−5 und 10−6 bestraft haben.

Wir verwendeten das grobkörnige Martini-Kraftfeld (v2.2)56,57,58,59, um den Aufbau des alpha-helikalen Transmembran-Homodimers Mga2 zu beschreiben. Alle MD-Simulationen wurden mit GROMACS v4.6.760,61,62,63 mit einem Integrationszeitschritt von Δt = 0,02 ps durchgeführt, wobei eine kubische Simulationsbox verwendet wurde, die die beiden identischen 30 Aminosäuren langen Alpha-Helices in einer hergestellten Lipidmembran enthielt aus 313 1-Palmitoyl-2-oleoyl-glycero-3-phosphocholin (POPC)-Molekülen. Die Membran überspannt den Kasten in der x-y-Ebene und wurde mit Wasser (5.348 Wasserkügelchen) und NaCl-Ionen entsprechend einer Konzentration von 150 mM (58 Na+, 60 Cl−) solvatisiert. Eine Referenztemperatur von T = 300 K wurde unter Verwendung des V-Rescale-Thermostats64 mit einer Kopplungskonstante von 1 ps getrennt für das Protein, die Membran und das Lösungsmittel erzwungen. Ein Druck von 1 bar wurde separat in der x-y-Ebene und in z unter Verwendung eines semiisotropen Parrinello-Rahman-Barostaten65 mit einer Zeitkonstante von 12 ps und einer Kompressibilität von 3 × 10−4 bar−1 ausgeübt. Jede MD-Simulation wurde auf einem einzelnen Rechenknoten mit zwei E5-2680-v3-CPUs und 64 GB Speicher durchgeführt. Das gesamte neuronale Netzwerktraining wurde auf einer RTX6000-GPU durchgeführt.

Um den Aufbau des Mga2-Homodimers zu beschreiben, verwendeten wir 28 interhelikale paarweise Abstände zwischen den Grundgerüstkügelchen der beiden Helices zusammen mit der Gesamtzahl der interhelikalen Kontakte, dem Abstand zwischen den Helix-Massenzentren und einer Reihe von handgeschneiderten Merkmalen, die das beschreiben Organisation der Lipide um die beiden Helices (Ergänzungstabelle 2). Um sicherzustellen, dass alle Netzwerkeingaben ungefähr in [0, 1] liegen, haben wir die Sigmoidfunktion \(f(r)={(1-(r/{R}_{0})^{6})}/{ verwendet. (1-(r/{R}_{0})^{12})}\) mit R0 = 2 nm für alle paarweisen Abstände, während wir alle Lipidmerkmale unter Verwendung der minimalen und maximalen Werte entlang des Übergangs skaliert haben. Die zusammengesetzten und zerlegten Zustände sind als Konfigurationen mit ≥130 interhelikalen Kontakten bzw. mit Helix-Helix-Schwerpunktabständen dCoM ≥ 3 nm definiert.

Das zur Anpassung des Committors verwendete neuronale Netzwerk wurde mit Keras52 implementiert und bestand aus einem anfänglichen dreischichtigen Pyramidenteil, in dem die Anzahl der Einheiten von 36 Eingaben auf 6 in der letzten Schicht unter Verwendung eines konstanten Faktors von (6/36)1 abnimmt /2 gefolgt von 6 Resteinheiten48,49, jeweils mit 4 Schichten und 6 Neuronen pro Schicht (Ergänzungstabelle 23). Auf die Eingaben wird ein Dropout von 0,01 angewendet und das Netzwerk wird mithilfe des Adam-Gradientenabstiegsprotokolls mit einer Lernrate von lr = 0,0001 trainiert (Ref. 51).

Um den Zusammenbaumechanismus von Mga2 zu untersuchen, führten wir parallel eine maschinengesteuerte Probenahme auf mehreren Knoten eines Hochleistungscomputerclusters durch. Wir haben 500 unabhängige TPS-Ketten nach dem aktuellen Committor-Modell betrieben. Die 500 TPS-Simulationen wurden mit zufälligen Anfangs-TPs initialisiert. Das zur Auswahl der ersten Schusspunkte verwendete neuronale Netzwerk wurde anhand vorläufiger Schussversuche (8.044 unabhängige Schüsse von 1.160 verschiedenen Punkten) trainiert. Nach zwei Runden (zwei Schritte in jeder der 500 unabhängigen TPS-Ketten) haben wir das Committor-Modell aktualisiert, indem wir auf allen neuen Daten trainiert haben. Nach der sechsten Runde haben wir erneut trainiert. Es war keine weitere Schulung erforderlich, wie aus der konsistenten Anzahl der erwarteten und beobachteten TP-Zählungen hervorgeht. Wir haben weitere 14 Runden für alle 500 TPS-Ketten durchgeführt, um TPs zu sammeln (Abb. 5b). Die Auswahl des Aufnahmepunkts, die TPS-Einrichtung und das Training des neuronalen Netzwerks wurden mithilfe von MDAnalysis66,67, numpy68 und unserem benutzerdefinierten Python-Paket vollständig im Python-Code automatisiert.

Die Input-Wichtigkeitsanalyse ergab, dass die Gesamtzahl der Kontakte ncontacts der wichtigste Input war (Extended Data Abb. 7). Allerdings konnte kein durch symbolische Regression als Funktion von ncontacts allein erzeugter Ausdruck den Committor genau reproduzieren. Es ist wahrscheinlich, dass ncontacts vom trainierten Netzwerk nur als binärer Schalter verwendet wird, um die beiden verschiedenen Regime zu unterscheiden – nahe an den gebundenen oder an den ungebundenen Zuständen. Indem wir die Eingabewichtigkeitsanalyse auf Trainingspunkte in der Nähe des ungebundenen Zustands beschränkten, stellten wir fest, dass das Netzwerk verschiedene interhelikale Kontakte verwendet, die annähernd ein helikales Muster nachzeichnen (Erweiterte Daten, Abb. 7). Wir haben eine symbolische Regression für alle möglichen Kombinationen durchgeführt, die aus einer, zwei oder drei der sieben wichtigsten Eingabekoordinaten bestehen (Ergänzungstabelle 3). Die besten Ausdrücke im Hinblick auf den Verlust wurden anhand von Validierungs-Committor-Daten ausgewählt, die während der Optimierung nicht verwendet wurden. Dieser Validierungssatz bestand aus Committor-Daten für 516 Konfigurationen mit jeweils 30 Testschüssen und 32 Konfigurationen mit 10 Testschüssen.

Um die Variabilität der beobachteten Reaktionsmechanismen abzuschätzen, führten wir eine hierarchische Gruppierung aller TPs durch, die in die durch die Kontakte 9 und 22 definierte Ebene projiziert werden und in die genaueste Parametrisierung eingehen, die durch symbolische Regression erzeugt wird. Anschließend haben wir dynamisches Time Warping69 verwendet, um die paarweise Ähnlichkeit zwischen allen TPs für das Clustering zu berechnen, was wir mit dem Scipy-Clustering-Modul70,71 durchgeführt haben. Die Pfaddichtediagramme (Abb. 5f, g) wurden entsprechend der Anzahl der Pfade und nicht der Anzahl der Konfigurationen histogrammiert, d. h. der Zähler jeder von einem bestimmten Pfad besuchten Zelle wurde für diesen Pfad um eins erhöht.

Trainingssatzdaten und Dateien zum Einrichten von Molekulardynamiksimulationen für den Aufbau von LiCl sind in der Code Ocean-Kapsel72 enthalten. Daten zur Reproduktion dieser Studie für alle verbleibenden Systeme (alle verbleibenden Ionen, Polymer, Clathrat und MGA2-Transmembrandimer) sind in einem Zenodo-Repository73 öffentlich verfügbar. Quelldaten werden mit diesem Dokument bereitgestellt.

Eine ausführbare Version des Codes „Artificial Intelligence for Molecular Mechanism Discovery“ (AIMMD) ist in der Code Ocean-Softwarekapsel72 verfügbar. Der AIMMD-Code ist auch unter https://github.com/bio-phys/aimmd als Git-Repository verfügbar.

Pena-Francesch, A., Jung, H., Demirel, MC & Sitti, M. Biosynthetische selbstheilende Materialien für weiche Maschinen. Nat. Mater. 19, 1230–1235 (2020).

Artikel Google Scholar

Van Driessche, AES et al. Molekulare Keimbildungsmechanismen und Kontrollstrategien für die Auswahl von Kristallpolymorphen. Natur 556, 89–94 (2018).

Artikel Google Scholar

Chung, HS, Piana-Agostinetti, S., Shaw, DE & Eaton, WA Struktureller Ursprung der langsamen Diffusion bei der Proteinfaltung. Science 349, 1504–1510 (2015).

Artikel Google Scholar

Dellago, C., Bolhuis, PG & Chandler, D. Effiziente Übergangspfad-Probenahme: Anwendung auf Lennard-Jones-Cluster-Umordnungen. J. Chem. Physik. 108, 9236–9245 (1998).

Artikel Google Scholar

Peters, B. & Trout, BL Erhalten von Reaktionskoordinaten durch Wahrscheinlichkeitsmaximierung. J. Chem. Physik. 125, 054108 (2006).

Artikel Google Scholar

Bolhuis, PG, Dellago, C. & Chandler, D. Reaktionskoordinaten der biomolekularen Isomerisierung. Proz. Natl Acad. Wissenschaft. USA 97, 5877–5882 (2000).

Artikel Google Scholar

Best, RB & Hummer, G. Reaktionskoordinaten und -raten von Übergangspfaden. Proz. Natl Acad. Wissenschaft. USA 102, 6732–6737 (2005).

Artikel Google Scholar

Berezhkovskii, AM & Szabo, A. Diffusion entlang der Spaltungs-/Commitment-Wahrscheinlichkeits-Reaktionskoordinate. J. Phys. Chem. B 117, 13115–13119 (2013).

Artikel Google Scholar

E, W. & Vanden-Eijnden, E. Auf dem Weg zu einer Theorie der Übergangspfade. J. Stat. Physik. 123, 503 (2006).

Artikel MathSciNet MATH Google Scholar

Krivov, SV Optimale Dimensionsreduktion komplexer Dynamiken: das Schachspiel als Diffusion in einer freien Energielandschaft. Physik. Rev. E 84, 011135 (2011).

Artikel Google Scholar

Hummer, G. Von Übergangspfaden zu Übergangszuständen und Ratenkoeffizienten. J. Chem. Physik. 120, 516–523 (2003).

Artikel Google Scholar

Bolhuis, PG, Chandler, D., Dellago, C. & Geissler, PL Probenahme von Übergangspfaden: Werfen von Seilen über raue Bergpässe, im Dunkeln. Annu. Rev. Phys. Chem. 53, 291–318 (2002).

Artikel Google Scholar

LeCun, Y., Bengio, Y. & Hinton, G. Deep Learning. Natur 521, 436–444 (2015).

Artikel Google Scholar

Schmidhuber, J. Deep Learning in neuronalen Netzen: ein Überblick. Neuronales Netz. 61, 85–117 (2015).

Artikel Google Scholar

Ma, A. & Dinner, AR Automatische Methode zur Identifizierung von Reaktionskoordinaten in komplexen Systemen. J. Phys. Chem. B 109, 6769–6779 (2005).

Artikel Google Scholar

Vanden-Eijnden, E., Venturoli, M., Ciccotti, G. & Elber, R. Zu den Annahmen, die der Meilensteinbildung zugrunde liegen. J. Chem. Physik. 129, 174102 (2008).

Artikel Google Scholar

Schmidt, M. & Lipson, H. Destillation freier Naturgesetze aus experimentellen Daten. Wissenschaft 324, 81–85 (2009).

Artikel Google Scholar

Ballard, AJ & Dellago, C. Auf dem Weg zum Mechanismus der Ionendissoziation in Wasser. J. Phys. Chem. B 116, 13490–13497 (2012).

Artikel Google Scholar

Behler, J. & Parrinello, M. Verallgemeinerte neuronale Netzwerkdarstellung hochdimensionaler Potentialenergieoberflächen. Physik. Rev. Lett. 98, 146401 (2007).

Artikel Google Scholar

Walsh, MR, Koh, CA, Sloan, ED, Sum, AK & Wu, DT Mikrosekundensimulationen der spontanen Keimbildung und des Wachstums von Methanhydrat. Science 326, 1095–1098 (2009).

Artikel Google Scholar

Arjun, Berendsen, TA & Bolhuis, PG Unvoreingenommene atomistische Einblicke in die konkurrierenden Keimbildungsmechanismen von Methanhydraten. Proz. Natl Acad. Wissenschaft. USA 116, 19305–19310 (2019).

Wang, Y., Herron, L. & Tiwary, P. Von Daten über Rauschen bis hin zu Daten zum Mischen von Physik über Temperaturen hinweg mit generativer künstlicher Intelligenz. Proz. Natl Acad. Wissenschaft. USA 119, e2203656119 (2022).

Artikel Google Scholar

Jacobson, LC, Hujo, W. & Molinero, V. Amorphe Vorläufer bei der Keimbildung von Clathrathydraten. Marmelade. Chem. Soc. 132, 11806–11811 (2010).

Artikel Google Scholar

Leitold, C. & Dellago, C. Faltungsmechanismus einer Polymerkette mit Anziehungskraft im Nahbereich. J. Chem. Physik. 141, 134901 (2014).

Artikel Google Scholar

Leitold, C., Lechner, W. & Dellago, C. Eine String-Reaktionskoordinate für die Faltung einer Polymerkette. J. Phys. Kondensiert. Matter 27, 194126 (2015).

Artikel Google Scholar

Covino, R. et al. Ein eukaryotischer Sensor für die Membranlipidsättigung. Mol. Zelle 63, 49–59 (2016).

Artikel Google Scholar

Chiavazzo, E. et al. Erkundung der intrinsischen Kartendynamik für unbekannte effektive Landschaften mit freier Energie. Proz. Natl Acad. Wissenschaft. USA 114, E5494–E5503 (2017).

Artikel Google Scholar

Bittracher, A. et al. Übergangsmannigfaltigkeiten komplexer metastabiler Systeme. J. Nichtlineare Wissenschaft. 28, 471–512 (2018).

Artikel MathSciNet MATH Google Scholar

Lucente, D., Duffner, S., Herbert, C., Rolland, J. & Bouchet, F. Maschinelles Lernen von Committor-Funktionen zur Vorhersage von Klimaereignissen mit großen Auswirkungen. Vorabdruck bei arXiv https://doi.org/10.48550/arXiv.1910.11736 (2019).

Wang, Y., Lamim Ribeiro, JM & Tiwary, P. Ansätze des maschinellen Lernens zur Analyse und Verbesserung von Molekulardynamiksimulationen. Curr. Meinung. Struktur. Biol. 61, 139–145 (2020).

Artikel Google Scholar

Noé, F., Tkatchenko, A., Müller, K.-R. & Clementi, C. Maschinelles Lernen für die molekulare Simulation. Annu. Rev. Phys. Chem. 71, 361–390 (2020).

Artikel Google Scholar

Noé, F., Olsson, S., Köhler, J. & Wu, H. Boltzmann-Generatoren: Probenahme von Gleichgewichtszuständen von Vielteilchensystemen mit Deep Learning. Wissenschaft 365, eaaw1147 (2019).

Rogal, J., Schneider, E. & Tuckerman, ME Auf neuronalen Netzwerken basierende kollektive Pfadvariablen für eine verbesserte Abtastung von Phasentransformationen. Physik. Rev. Lett. 123, 245701 (2019).

Artikel Google Scholar

Sidky, H., Chen, W. & Ferguson, AL Maschinelles Lernen für die Entdeckung kollektiver Variablen und verbesserte Probenentnahme in der biomolekularen Simulation. Mol. Physik. 118, e1737742 (2020).

Artikel Google Scholar

Bartók, AP et al. Maschinelles Lernen vereinheitlicht die Modellierung von Materialien und Molekülen. Wissenschaft. Adv. 3, e1701816 (2017).

Artikel Google Scholar

Udrescu, S.-M. & Tegmark, M. AI Feynman: eine von der Physik inspirierte Methode zur symbolischen Regression. Wissenschaft. Adv. 6, eaay2631 (2020).

Artikel Google Scholar

Jung, H., Okazaki, K.-i & Hummer, G. Übergangspfad-Sampling seltener Ereignisse durch Aufnahme von oben. J. Chem. Physik. 147, 152716 (2017).

Artikel Google Scholar

Kemp, SJ, Zaradic, P. & Hansen, F. Ein Ansatz zur Bestimmung der relativen Wichtigkeit und Signifikanz von Eingabeparametern in künstlichen neuronalen Netzen. Ökologisch. Modell. 204, 326–334 (2007).

Artikel Google Scholar

Izzo, D. & Biscani, F. dcgp: Differenzierbare kartesische genetische Programmierung leicht gemacht. J. Open Source Softw. 5, 2290 (2020).

Artikel Google Scholar

Joung, IS & Cheatham, TE Bestimmung der Parameter monovalenter Alkali- und Halogenidionen zur Verwendung in explizit solvatisierten biomolekularen Simulationen. J. Phys. Chem. B 112, 9020–9041 (2008).

Artikel Google Scholar

Jorgensen, WL, Chandrasekhar, J., Madura, JD, Impey, RW & Klein, ML Vergleich einfacher Potentialfunktionen zur Simulation von flüssigem Wasser. J. Chem. Physik. 79, 926–935 (1983).

Artikel Google Scholar

Eastman, P. et al. OpenMM 7: schnelle Entwicklung von Hochleistungsalgorithmen für die Molekulardynamik. PLoS Comput. Biol. 13, e1005659 (2017).

Sivak, DA, Chodera, JD & Crooks, GE Zeitschritt-Neuskalierung stellt zeitkontinuierliche dynamische Eigenschaften für die zeitdiskrete Langevin-Integration von Nichtgleichgewichtssystemen wieder her. J. Phys. Chem. B 118, 6466–6474 (2014).

Essmann, U. et al. Eine Ewald-Methode mit glattem Partikelnetz. J. Chem. Physik. 103, 8577–8593 (1995).

Artikel Google Scholar

Behler, J. Atomzentrierte Symmetriefunktionen zur Konstruktion hochdimensionaler neuronaler Netzwerkpotentiale. J. Chem. Physik. 134, 074106 (2011).

Artikel Google Scholar

Behler, J. Darstellung potentieller Energieoberflächen durch hochdimensionale neuronale Netzwerkpotentiale. J. Phys. Kondensiert. Materie 26, 183001 (2014).

Artikel Google Scholar

Geiger, P. & Dellago, C. Neuronale Netze zur lokalen Strukturerkennung in polymorphen Systemen. J. Chem. Physik. 139, 164105 (2013).

Artikel Google Scholar

He, K., Zhang, X., Ren, S. & Sun, J. Deep Residual Learning für die Bilderkennung. Im Jahr 2016 IEEE-Konferenz zu Computer Vision und Mustererkennung 770–778 (IEEE, 2016).

He, K., Zhang, X., Ren, S. & Sun, J. Identitätszuordnungen in tiefen Restnetzwerken. In Computer Vision – ECCV 2016, Lecture Notes in Computer Science (Hrsg. Leibe, B. et al.) 630–645 (Springer, 2016).

Klambauer, G., Unterthiner, T., Mayr, A. & Hochreiter, S. Selbstnormalisierende neuronale Netze. In Proc. 31. Internationale Konferenz über neuronale Informationsverarbeitungssysteme 972–981 (Curran Associates, 2017).

Kingma, DP & Ba, J. Adam: eine Methode zur stochastischen Optimierung. Vorabdruck bei 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: ein Python-Framework für Path-Sampling-Simulationen. 1. Grundlagen. J. Chem. Theorieberechnung. 15, 813–836 (2019).

Artikel Google Scholar

Swenson, DWH, Prinz, J.-H., Noe, F., Chodera, JD & Bolhuis, PG openpathsampling: ein Python-Framework für Path-Sampling-Simulationen. 2. Erstellen und Anpassen von Pfadensembles und Beispielschemata. J. Chem. Theorieberechnung. 15, 837–856 (2019).

Artikel Google Scholar

Lechner, W. & Dellago, C. Genaue Bestimmung von Kristallstrukturen basierend auf gemittelten Parametern der lokalen Bindungsordnung. J. Chem. Physik. 129, 114707 (2008).

Artikel Google Scholar

Marrink, SJ, de Vries, AH & Mark, AE Grobkörniges Modell für semiquantitative Lipidsimulationen. J. Phys. Chem. B 108, 750–760 (2004).

Artikel Google Scholar

Marrink, SJ, Risselada, HJ, Yefimov, S., Tieleman, DP & de Vries, AH Das MARTINI-Kraftfeld: grobkörniges Modell für biomolekulare Simulationen. J. Phys. Chem. B 111, 7812–7824 (2007).

Artikel Google Scholar

Monticelli, L. et al. Das grobkörnige MARTINI-Kraftfeld: Erweiterung auf Proteine. J. Chem. Theorieberechnung. 4, 819–834 (2008).

Artikel Google Scholar

de Jong, DH et al. Verbesserte Parameter für das grobkörnige Martini-Protein-Kraftfeld. J. Chem. Theorieberechnung. 9, 687–697 (2013).

Artikel Google Scholar

Berendsen, H., van der Spoel, D. & van Drunen, R. GROMACS: eine Implementierung der parallelen Molekulardynamik zur Nachrichtenübermittlung. Berechnen. Physik. Komm. 91, 43–56 (1995).

Artikel Google Scholar

Hess, B., Kutzner, C., van der Spoel, D. & Lindahl, E. GROMACS 4: Algorithmen für hocheffiziente, lastenausgleichende und skalierbare molekulare Simulation. J. Chem. Theorieberechnung. 4, 435–447 (2008).

Artikel Google Scholar

Pronk, S. et al. GROMACS 4.5: ein Open-Source-Toolkit für die molekulare Simulation mit hohem Durchsatz und hoher Parallelität. Bioinformatik 29, 845–854 (2013).

Artikel Google Scholar

Abraham, MJ et al. GROMACS: Hochleistungsmolekularsimulationen durch mehrstufige Parallelität vom Laptop bis zum Supercomputer. SoftwareX 1–2, 19–25 (2015).

Artikel Google Scholar

Bussi, G., Donadio, D. & Parrinello, M. Kanonische Abtastung durch Geschwindigkeitsneuskalierung. J. Chem. Physik. 126, 014101 (2007).

Artikel Google Scholar

Parrinello, M. & Rahman, A. Polymorphe Übergänge in Einkristallen: eine neue Methode der Molekulardynamik. J. Appl. Physik. 52, 7182–7190 (1981).

Artikel Google Scholar

Michaud-Agrawal, N., Denning, EJ, Woolf, TB & Beckstein, O. MDAnalysis: ein Toolkit für die Analyse von Molekulardynamiksimulationen. J. Comput. Chem. 32, 2319–2327 (2011).

Artikel Google Scholar

Gowers, R. et al. MDAnalysis: ein Python-Paket zur schnellen Analyse von Molekulardynamiksimulationen. Proceedings of the Python in Science Conferences (Hrsg. Benthall, S. & Rostrup, S.) 98–105 (2016); https://doi.org/10.25080/issn.2575-9752

Harris, CR et al. Array-Programmierung mit NumPy. Natur 585, 357–362 (2020).

Artikel Google Scholar

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

Virtanen, P. et al. SciPy 1.0: grundlegende Algorithmen für wissenschaftliches Rechnen in Python. Nat. Methoden 17, 261–272 (2020).

Artikel Google Scholar

Müllner, D. Modern hierarchical, agglomerative clustering algorithms. Preprint at arXiv https://doi.org/10.48550/arXiv.1109.2378 (2011).

Jung, H. et al. Maschinengesteuerte Pfadprobenahme zur Entdeckung von Mechanismen der molekularen Selbstorganisation (Softwarekapsel). Code Ocean https://doi.org/10.24433/CO.7949737.v1 (2023).

Jung, H. et al. Maschinengesteuerte Pfadprobenahme zur Entdeckung von Mechanismen der molekularen Selbstorganisation (Trainings- und Validierungsdaten). Zenodo https://doi.org/10.5281/zenodo.7704969 (2023).

Referenzen herunterladen

Wir danken FE Blanc für nützliche Kommentare und der Openpathsampling-Community, insbesondere D. Swenson, für Diskussionen und technische Unterstützung. HJ, RC und GH danken der Max-Planck-Gesellschaft für ihre Unterstützung. RC dankt dem Frankfurt Institute for Advanced Studies für seine Unterstützung. RC und GH bedanken sich für die Unterstützung durch das LOEWE CMMS-Programm des Landes Hessen, den Sonderforschungsbereich 1507 der Deutschen Forschungsgemeinschaft und die International Max Planck Research School on Cellular Biophysics. AA und PGB danken für die Unterstützung des CSER-Programms der niederländischen Organisation für wissenschaftliche Forschung (NWO) und von Shell Global Solutions International. BV, CL und CD bedanken sich für die Unterstützung des Österreichischen Wissenschaftsfonds (FWF) durch den SFB TACO, Fördernummer F 81-N. Die Geldgeber hatten keinen Einfluss auf das Studiendesign, die Datenerfassung und -analyse, die Entscheidung zur Veröffentlichung oder die Erstellung des Manuskripts.

Open-Access-Förderung durch die Max-Planck-Gesellschaft.

Diese Autoren haben gleichermaßen beigetragen: Hendrik Jung, Roberto Covino.

Abteilung für Theoretische Biophysik, Max-Planck-Institut für Biophysik, Frankfurt am Main, Deutschland

Hendrik Jung & Gerhard Hummer

Frankfurt Institute for Advanced Studies, Frankfurt am Main, Deutschland

Roberto Covino

van 't Hoff Institut für Molekularwissenschaften, Universität Amsterdam, Amsterdam, Niederlande

A. Arjun und Peter G. Bolhuis

Fakultät für Physik, Universität Wien, Wien, Österreich

Christian Leitold & Christoph Dellago

Institut für Biophysik, Goethe-Universität Frankfurt, Frankfurt am Main, Deutschland

Gerhard Hummer

Sie können diesen Autor auch in PubMed Google Scholar suchen

Sie können diesen Autor auch in PubMed Google Scholar suchen

Sie können diesen Autor auch in PubMed Google Scholar suchen

Sie können diesen Autor auch in PubMed Google Scholar suchen

Sie können diesen Autor auch in PubMed Google Scholar suchen

Sie können diesen Autor auch in PubMed Google Scholar suchen

Sie können diesen Autor auch in PubMed Google Scholar suchen

HJ entwickelte den Code für maschinelles Lernen, führte Ionen- und Proteinassemblierungssimulationen durch und analysierte die Daten zusammen mit RC und GHAA führte die Gashydratsimulationen durch und analysierte die Daten zusammen mit HJ, RC und PGBCL führten die Polymerfaltungssimulationen durch und analysierten die Daten gemeinsam mit HJ, RC und CDGH konzipierten die Studie. HJ, RC und GH haben die Forschung entworfen. HJ, RC und GH haben das Manuskript unter Mitwirkung aller Autoren verfasst. RC, PGB, CD und GH koordinierten das Projekt.

Korrespondenz mit Gerhard Hummer.

Die Autoren geben an, dass keine Interessenkonflikte bestehen.

Hauptredakteurin für Handhabung: Kaitlin McCardle, in Zusammenarbeit mit dem Nature Computational Science-Team.

Anmerkung des Herausgebers Springer Nature bleibt hinsichtlich der Zuständigkeitsansprüche in veröffentlichten Karten und institutionellen Zugehörigkeiten neutral.

Normalisierte Autokorrelationsfunktion der Übergangszeit tTP als Funktion der Anzahl der Monte-Carlo-Schritte (MC) für TPS-Simulationen an verschiedenen Ionenspezies, bei denen die Schusspunkte durch den Algorithmus (rote und grüne Kurven, unabhängige Läufe) oder durch ausgewählt wurden zufällige einheitliche Auswahl (orange und blaue Kurve, unabhängige Läufe). Jede Autokorrelationsfunktion wurde über MC-Ketten mit einer Länge von jeweils 100.000 Schritten berechnet. Die Dekorrelationszeiten der Markov-Ketten, definiert als die Punkte, an denen die normalisierten Autokorrelationsfunktionen den Wert e−1 erreichen, sind mit gepunkteten vertikalen Linien in der gleichen Farbe wie die jeweilige Autokorrelationsfunktion markiert. Maschinengeführte Probenahme führt in allen Fällen zu Dekorrelationszeiten von vier MC-Schritten, während dieser Wert bei zufälliger gleichmäßiger Auswahl je nach Ionenart zwischen 20 und 60 Schritten liegt.

Quelldaten

Committor-Kreuzkorrelationsdiagramme für Modelle, die durch Transfertraining an verschiedenen Ionenspezies mit unterschiedlicher Anzahl von Trainingsschusspunkten erhalten wurden. Das ursprüngliche Modell wurde auf LiCl trainiert und ist in allen Fällen gleich (siehe Abb. 1 für ein Committor-Kreuzkorrelationsdiagramm des ursprünglichen Modells zur Vorhersage von LiCl). Das Modell hat die Architektur „ResNet I“ und verwendet die Eingaben „SF longranged II“ (siehe Ergänzungstabellen 5 und 6). Das Transfertraining erfolgte durch Randomisierung der letzten Schicht des Originalmodells mit einem einzelnen Neuron und anschließendes Training an der angegebenen Anzahl von Schusspunkten mit Daten für jedes Ionensystem separat. Der Durchschnitt der Stichproben-Committors (blaue Linie) +/- eine SD (orange schattiert) wurde berechnet, indem die Stichproben-Committors in den durch die vertikalen Schritte angegebenen Bereich des erlernten Committors eingeteilt wurden. Als Referenz zeigt die rote Linie die Identität an.

Quelldaten

Dargestellt sind Kreuzkorrelationen von abgetasteten Committoren (x-Achse) und vorhergesagten Committoren (y-Achse) für Modelle, die auf verschiedenen Teilmengen der sechs verschiedenen Ionenpaare trainiert wurden. (Zeile 1) Training für alle sechs Ionenpaare (Li+ Cl−, Li+ I−, Na+ Cl−, Na+ I−, Cs+ Cl− und Cs+ I−) gleichzeitig. (Zeile 2) Modell ohne Daten für Li+ trainiert, d. h. ohne Li+ Cl− und Li+ I−. (Zeile 3) Modell ohne Daten für Na+ trainiert, d. h. ohne Na+ Cl− und Na+ I−. (Zeile 4) Modell ohne Daten für Cl− trainiert, d. h. ohne Li+ Cl−, Na+ Cl− und Cs+ Cl−. (Zeile 5) Modell ohne Daten für I− trainiert, d. h. ohne Li+ I−, Na+ I− und Cs+ I−. Alle Modelle haben die Architektur „ResNet I“ (siehe Ergänzungstabelle 6). Für alle Modelle wurden die verschiedenen Ionenspezies unterschieden, indem die Lennard-Jones-Parameter ϵ und σ des Kraftfelds als zusätzliche Deskriptoren zum Satz „SF shortranged“ (siehe Ergänzungstabelle 5) hinzugefügt wurden, der zur Beschreibung des Lösungsmittels um das Kraftfeld herum verwendet wurde Ionenpaar. Der Durchschnitt der Stichproben-Committors (blaue Linie) +/- eine SD (orange schattiert) wurde berechnet, indem die Stichproben-Committors in den durch die vertikalen Schritte angegebenen Bereich des erlernten Committors eingeteilt wurden. Als Referenz zeigt die rote Linie die Identität an.

Quelldaten

a, Validierung des auf mehrere Ionen reduzierten Modells für weitere vier Ionenpaare. Der Durchschnitt der ausgewählten Committoren (blaue Linie) und ihre Standardabweichung (orange schattiert) werden durch Binning entlang des vorhergesagten Committors (rote Linie: Identität) berechnet. b, Wichtigste Eingabekoordinaten, die den Committor bestimmen, der auf Assoziationssimulationen eines einzelnen Ionenpaars, in diesem Fall LiCl, trainiert wurde (siehe auch Ergänzungstabelle 13 für eine Auflistung der zehn relevantesten Koordinaten). c, Reduzierte Modelle q0 und q1, die die Assoziation von Li+Cl− in Wasser beschreiben, erhalten durch symbolische Regression bei strikter (λ = 10−6) bzw. sanfter Regularisierung (λ = 10−7). Beachten Sie, dass das erste Modell nicht von den Freiheitsgraden des Wassers abhängt (eine Beschreibung der Koordinaten finden Sie noch einmal in der Ergänzungstabelle 13). d, Kreuzkorrelationsdiagramme zwischen untrainierten Committor-Daten und den symbolischen Regressionsvorhersagen als unabhängige Validierungen der Genauigkeit von q0 und q1. Der Durchschnitt der ausgewählten Committoren (blaue Linie) +/- eine SD (orange schattiert) wird durch Binning entlang des vorhergesagten Committors (rote Linie: Identität) berechnet.

Quelldaten

Kreuzkorrelation zwischen erlerntem Committor und dem Committor, die durch wiederholte Stichprobenentnahme an nicht trainierten Konfigurationen für zwei Modelle erhalten wurde, die nur auf drei der vier im Trainingssatz verfügbaren Temperaturen trainiert wurden. (Links) Das Committor-Modell wurde nur auf Daten für T=270 K, 275 K und 285 K trainiert, um die Fähigkeit des Modells zur Interpolation auf T=280 K zu beurteilen. (Rechts) Modell wurde auf Daten für T=270 K, 275 trainiert K und 280 K, um die Fähigkeit des Modells zu beurteilen, auf T=285 K zu extrapolieren. Die rote Linie stellt die Identität dar.

Quelldaten

a, Committor-Kreuzkorrelation bei nicht trainierten Konfigurationen. b, Eingaberelevanz für neuronale Netze, die auf Polymerkeimbildungsdaten unter Verwendung von Merkmalen mit niedriger bzw. hoher Auflösung trainiert wurden. In a wurden der Durchschnitt der ausgewählten Committoren (blaue Linie) und ihre Standardabweichung (orange schattiert) durch Binning entlang des vorhergesagten Committors (rote Linie: Identität) berechnet. In b zeigen die blauen Balken den Mittelwert von n unabhängigen Wiederholungen der Eingabewichtigkeitsanalyse (n = 100 für Features mit hoher Auflösung und n = 250 für Features mit niedriger Auflösung), und die Fehlerbalken zeigen +/- eine SD an.

Quelldaten

Eingabewichtigkeitsanalysen unter Verwendung aller Trainingspunkte (oberes Feld) und einer Teilmenge mit ncontacts < 0,01 (unteres Feld), entsprechend Trainingspunkten nahe dem ungebundenen Zustand. Die Höhe jedes Balkens ist der Mittelwert über 50 unabhängige Analysen (n = 50), während die Balken +/- eine SD anzeigen. Alle Werte werden auf die größte Wichtigkeit in jedem Satz normalisiert.

Quelldaten

Committor-Kreuzkorrelationsdiagramm für den symbolischen Regressionsausdruck \({q}_{B}({x}_{9},{x}_{2}2)=-\exp ({x}_{9}^{2 })\log ({x}_{9}-\frac{{x}_{9}}{\log ({x}_{22})})\) auf untrainierten Validierungs-Committor-Daten der Mga2-Transmembrananordnung. Der Ausdruck ist eine Funktion des interhelikalen Kontakts 9 (x9) oben an den beiden Helices und des Kontakts 22 (x22) unten (siehe Ergänzungstabellen 2 und 3). Der Durchschnitt der ausgewählten Committoren (blaue Linie) und ihre Standardabweichung (orange schattiert) werden durch Binning entlang des vorhergesagten Committors (rote Linie: Identität) berechnet.

Quelldaten

Trainingsiterationen für den Li+Cl−-Zusammenbau. Die blaue Linie zeigt die aus dem Effizienzfaktor berechnete Lernrate bei jedem Schritt, orangefarbene Kreuze zeigen, wann tatsächlich trainiert wurde. Der Einschub zeigt den Trainingsverlust pro Schusspunkt für jedes Training. Es werden nur die ersten 26000 Monte-Carlo-Schritte angezeigt. Das Modell verwendet die Eingaben „SF longranged II“ und hat die Architektur „ResNet I“ (siehe Ergänzungstabellen 5 und 6).

Quelldaten

Ergänzungstabellen 1–23.

Effizienz von Ionen-TPS-Simulationen: Anzahl der akzeptierten, erzeugten und erwarteten Übergänge pro Monte-Carlo-Schritt für alle durch maschinelles Lernen unterstützten TPS-Simulationen für alle Ionenspezies. Der in der letzten Spalte angezeigte Unterschied besteht zwischen erwarteten und generierten Übergängen, also (erwartet − generiert)/Schritten. Die Werte für die Anfangsphase (Endphase) wurden über die ersten (letzten) 1.000 MC-Schritte berechnet, der Wert für die vollständige Simulation beträgt über alle 100.000 MC-Schritte.

Statistische Quelldaten für Abb. 1b–e.

Statistische Quelldaten für Abb. 2a, c, e.

Statistische Quelldaten für Abb. 3b, c.

Statistische Quelldaten für Abb. 4a,b.

Statistische Quelldaten für Abb. 5b, c, e–g.

Statistische Quelldaten für Extended Data Abb. 1.

Statistische Quelldaten für Extended Data Abb. 2.

Statistische Quelldaten für Extended Data Abb. 3.

Statistische Quelldaten für Extended Data Abb. 4a,b,d.

Statistische Quelldaten für Extended Data Abb. 5.

Statistische Quelldaten für Extended Data Abb. 6.

Statistische Quelldaten für Extended Data Abb. 7.

Statistische Quelldaten für Extended Data Abb. 8.

Statistische Quelldaten für Extended Data Abb. 9.

Open Access Dieser Artikel ist unter einer Creative Commons Attribution 4.0 International License lizenziert, die die Nutzung, Weitergabe, Anpassung, Verbreitung und Reproduktion in jedem Medium oder Format erlaubt, sofern Sie den/die Originalautor(en) und die Quelle angemessen angeben. Geben Sie einen Link zur Creative Commons-Lizenz an und geben Sie an, ob Änderungen vorgenommen wurden. Die Bilder oder anderes Material Dritter in diesem Artikel sind in der Creative Commons-Lizenz des Artikels enthalten, sofern in der Quellenangabe für das Material nichts anderes angegeben ist. Wenn Material nicht in der Creative-Commons-Lizenz des Artikels enthalten ist und Ihre beabsichtigte Nutzung nicht gesetzlich zulässig ist oder über die zulässige Nutzung hinausgeht, müssen Sie die Genehmigung direkt vom Urheberrechtsinhaber einholen. Um eine Kopie dieser Lizenz anzuzeigen, besuchen Sie http://creativecommons.org/licenses/by/4.0/.

Nachdrucke und Genehmigungen

Jung, H., Covino, R., Arjun, A. et al. Maschinengesteuerte Pfadprobenahme zur Entdeckung von Mechanismen der molekularen Selbstorganisation. Nat Comput Sci 3, 334–345 (2023). https://doi.org/10.1038/s43588-023-00428-z

Zitat herunterladen

Eingegangen: 15. Februar 2023

Angenommen: 10. März 2023

Veröffentlicht: 24. April 2023

Ausgabedatum: April 2023

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

Jeder, mit dem Sie den folgenden Link teilen, kann diesen Inhalt lesen:

Leider ist für diesen Artikel derzeit kein Link zum Teilen verfügbar.

Bereitgestellt von der Content-Sharing-Initiative Springer Nature SharedIt