Wyniki i dyskusja
Zdarzenie autojonizacji zostało zbadane przy użyciu symulacji ab initio RETIS, jak opisano w Materiałach i Metodach. W symulacjach RETIS zastosowaliśmy stosunkowo prosty geometryczny parametr porządku odległości, λ, jak pokazano na Rys. 1: Gdy układ składa się tylko z H2O, λ jest największą odległością wiązania kowalencyjnego O-H, a gdy układ zawiera OH- i H3O+, λ jest najkrótszą odległością pomiędzy tlenem w OH- i atomami wodoru w H3O+. W dalszej części pracy atom tlenu używany do określenia parametru porządku nazywamy Oλ. Rodzaj gatunku (OH-, H2O, lub H3O+) został określony poprzez przypisanie każdemu wodorowi pojedynczego wiązania łączącego go z najbliższym tlenem. Zauważmy, że definicja parametru porządku nie wymaga progu dla zdefiniowania wiązania chemicznego, ani nie ogranicza parametru porządku do konkretnych cząsteczek wody na czas trwania symulacji. Oznacza to, że obliczamy szybkość dysocjacji dowolnej cząsteczki wody w układzie zamiast pojedynczego ukierunkowanego wiązania O-H lub cząsteczki wody.
Warunki inicjacji i lokalne zmienne kolektywne. (A) Rozkłady reaktywne (rλc,λr(ξ)) i niereaktywne (uλc,λr(ξ)) dla ξ={w4,na} oraz λc=1.16 Å i λr=2.0 Å. Dla celów wizualizacji przedstawione rozkłady zostały znormalizowane . Wstawki górna i prawa przedstawiają jednowymiarowe projekcje rozkładów. Widać wyraźne rozdzielenie obu rozkładów wzdłuż współrzędnej w4, co wskazuje, że trajektorie reaktywne są bardziej skompresowane w porównaniu z trajektoriami niereaktywnymi. Dodatkowo, atom tlenu użyty do obliczenia parametru porządku (Oλ) przyjmuje średnio większą liczbę wiązań wodorowych w trajektoriach reaktywnych, w porównaniu z trajektoriami niereaktywnymi. (B) Ilustracyjna migawka z trajektorii reaktywnej, gdzie Oλ jest pokazany na niebiesko. Cztery otaczające atomy tlenu, które są używane do obliczenia tetraedrycznego parametru q, są pokazane na pomarańczowo. Przewód wodny jest zaznaczony żółtą linią (i szarymi przezroczystymi kulami), a parametr kątowy qcos jest wskazany. W tym ujęciu przewód wodny jest ściśnięty, q wykazuje odchylenie od struktury tetraedrycznej, qcos wskazuje, że trzy atomy tlenu ustawiają się w przewodzie, a Oλ przyjmuje trzy wiązania wodorowe i jedno oddaje (pokazane zielonymi liniami).
Jeżeli weźmiemy pod uwagę współrzędną q, zauważymy, że rλc,λr jest przesunięta w kierunku niższych wartości q w porównaniu z uλc,λr, co wskazuje, że zniekształcenie z układu tetraedrycznego wokół dysocjującej formy wody może również inicjować zdarzenie. Odkrycie to jest nieco zaskakujące, gdyż w niektórych innych reakcjach chemicznych w fazie wodnej stwierdzono odwrotny efekt (31). Podobne wnioski można wyciągnąć dla rozkładu ξ=(w4,qcos). W tym przypadku występuje pik wzdłuż współrzędnej qcos dla rozkładu reaktywnego bliższego liniowemu ułożeniu cząsteczek wody. Na Rys. 4B przedstawiamy reprezentatywną migawkę, otrzymaną we wczesnym stadium (po 3 fs) trajektorii reaktywnej. Ogólnie rzecz biorąc, wyniki przedstawione na Rys. 3 wskazują, że kompresja drutu wodnego (mierzona przez w4) i hiperkoordynacja (mierzona przez na) lub dystorsja (mierzona przez q i qcos) są koniecznymi warunkami inicjacji dla autojonizacji. Jednakże, nie są to warunki wystarczające, jak pokazują wartości TAλc,λr na Rys. 3B: Nadal 60% trajektorii rozpoczynających się w idealnym zakresie parametrów ξ nie ustanawia zgodnego skoku protonu.
Uczenie maszynowe (ML) zastosowane do danych próbkowania ścieżek (33, 34) jest obiecującym podejściem do znalezienia ważnych zmiennych zbiorczych, które mogą być łatwo przeoczone przez ludzką intuicję. Aby zbadać tę możliwość, zbudowaliśmy modele ML do przewidywania wyniku trajektorii, biorąc pod uwagę stan systemu wodnego na początku trajektorii. Skupiamy się na tym samym zakresie, co w analizie mocy predykcyjnej i używamy stanu systemu, kiedy λ>1.15 Å jest po raz pierwszy osiągnięta, aby przewidzieć wynik. Zastosowaliśmy kilka technik ML, w których każdy zespół ścieżek nieparzystych został włączony do kalibracji, a zespoły ścieżek parzystych zostały użyte do zbioru testowego. Alternatywny podział, w którym dane w ramach każdego zespołu ścieżek zostały równo podzielone na dwie części, dał podobne wyniki. Ponadto, ponieważ silnie skośne rozkłady są trudne do traktowania za pomocą ML, pominęliśmy ponowne ważenie zbiorów danych z wagami statystycznymi odpowiednich zespołów ścieżek. Jednakże, zastosowaliśmy techniki ML jako podejście jakościowe w celu znalezienia nowych parametrów, które mogłyby być testowane ilościowo w ramach metody mocy predykcyjnej (19).
Dodatkowo, aby uniknąć potencjalnego ryzyka nadinterpretacji, zdecydowaliśmy się ograniczyć złożoność procesu decyzyjnego ML i narzuciliśmy maksymalnie cztery parametry rzędu podczas obliczania TAλc,λr. Na przykład, doskonałe wyniki predykcyjne (>90%) uzyskano przy użyciu maszyn gradientowo-boostingowych opartych na zespole (35, 36). Interpretacja modelu jest jednak problematyczna, ponieważ wykorzystuje się zespół 100-150 głębokich drzew decyzyjnych (dodawanych w sekwencji). Mimo poprawy wydajności, wzrasta szansa na przepasowanie z przypadkowymi korelacjami. W związku z tym ograniczyliśmy się do modeli decyzyjnych opartych na pojedynczych drzewach decyzyjnych klasyfikacji i regresji (CART) (20). Ograniczenie do czterech parametrów porządkowych dla funkcji TAλc,λr wynika z podobnych przesłanek. Dodanie większej liczby parametrów daje bardziej rzadkie macierze reprezentujące rozkłady reaktywny/niereaktywny, a w rezultacie całkowanie numeryczne do obliczania nakładania się tych rozkładów staje się bardzo wrażliwe na rozmiar bloku i może zaniżać nakładanie się z powodu pustych bloków z powodu niewystarczającej statystyki.
Rozważaliśmy 138 zmiennych zbiorczych składających się z odległości tlen-tlen; odległości tlen-wodór dla początkowo związanych cząsteczek wody; wszystkie kąty utworzone przez Oλ i jego czterech najbliższych sąsiadów tlenowych; oraz parametry porządku Steinhardta dla zamówień 3, 4, i 6 (32) (patrz Materiały i Metody po więcej szczegółów). Dodatkowo, dodano parametry porządkowe już uwzględnione. Rys. 5A przedstawia wynikowe drzewo decyzyjne. Co godne uwagi, spośród wszystkich parametrów wejściowych, parametr w4 znajduje się zarówno na szczycie drzewa decyzyjnego, jak i jest najważniejszą zmienną mierzoną redukcją błędu klasyfikacji przypisywanego każdej zmiennej przy każdym podziale drzewa decyzyjnego (20) (Dodatek SI, Rys. S9). W drzewie decyzyjnym pojawia się również uporządkowanie tetraedryczne i liczba akceptowanych wiązań wodorowych. Aby opisać pierwszy efekt, podejście ML nadało priorytet parametrowi uporządkowania Steinhardta q4 powyżej podobnego parametru q używanego przez nas wcześniej. Niektóre odległości, które również pojawiają się w drzewie decyzyjnym, takie jak d25, odległość pomiędzy Oλ i jego 25. najbliższym tlenem, są najprawdopodobniej wynikiem przypadkowych korelacji spowodowanych ograniczonym rozmiarem zbioru danych. Można to zweryfikować sprawdzając ważność tej zmiennej: d25 nie pojawia się wśród 20 najważniejszych zmiennych (Dodatek SI, Rys. S9), a w rzeczywistości inne podobne zmienne (np. d24) znajdują się wyżej w rankingu, choć z niską ważnością. Ważniejszym i intuicyjnie uzasadnionym parametrem sugerowanym przez podejście ML jest λ2, odległość OH pomiędzy tlenem najbliższym Oλ a jego wodorem o największym wiązaniu wewnątrzcząsteczkowym. Ponowne obliczenie zdolności predykcyjnej przy użyciu parametrów z drzewa ML (Rys. 5B) nie dało lepszych wyników niż kombinacja w4, q, na, i qcos, ale powinno być postrzegane jako równie dobre, biorąc pod uwagę niepewności statystyczne.
Wyniki analizy uczenia maszynowego. (A) Drzewo klasyfikacyjne i regresyjne dla przewidywania wyniku zainicjowanych trajektorii. Uwzględniliśmy tu kilka dodatkowych zmiennych zbiorczych (opis w Materials and Methods), ale tylko niewielki ich podzbiór jest ostatecznie potrzebny do skonstruowania drzewa: w4, q4 , λ2 (długość rozciągniętego wiązania wodorowego w cząsteczce wody najbliższej gatunkowi Oλ), di (odległość od Oλ do i-tego najbliższego tlenu) oraz d¯i (średnia odległość uwzględniająca i najbliższe oksygeny). Notacja dla węzłów jest wyjaśniona za pomocą samodzielnego węzła w lewym górnym rogu. Drzewo to przewiduje trajektorie jako reaktywne, tj. osiągające λ≥2, lub niereaktywne w oparciu o zmienne zbiorcze uzyskane w ramce w trajektorii, gdy λ jest po raz pierwszy ≥1.15. Węzły przewidujące trajektorie reaktywne są oznaczone kolorem niebieskim (klasa 1), natomiast węzły przewidujące trajektorie niereaktywne są oznaczone kolorem zielonym (klasa 0). Zwróć uwagę, że procenty na dole kwadratów nie odzwierciedlają fizycznie poprawnych frakcji, ponieważ zespoły ścieżek nie były ponownie ważone przy użyciu ich wag statystycznych. Reguły są tekstowymi reprezentacjami przemierzania drzewa; na przykład, reguła 5 (która przewiduje trajektorie reaktywne) może być wyrażona jako w4≥7.6 i λ2≥1.1. Reguły te dają różne warunki inicjacji i są one wymienione w Dodatku SI, Tabela S1, dla dolnego rzędu węzłów. (B) Moc predykcyjna i prawdopodobieństwo przecięcia jako funkcja λr dla λc=1.16 Å i różnych kombinacji zmiennych kolektywnych. Tutaj porównujemy moc predykcyjną przy użyciu zidentyfikowanych przez nas zmiennych kolektywnych ze zmiennymi oznaczonymi jako ważne przez analizę uczenia maszynowego. (C) Reaktywne (rλc,λr(ξ)) i niereaktywne (uλc,λr(ξ)) rozkłady dla ξ={λ2,d¯2} oraz λc=1.16 Å i λr=2.0 Å. Dla celów wizualizacji, przedstawione rozkłady zostały znormalizowane. Górna i prawa wstawka pokazują jednowymiarowe projekcje rozkładów.