Results and Discussion
材料と方法にあるように、第一原理RETISシミュレーションを用いて、自己イオン化現象を調べた。 RETISシミュレーションでは、図1に示すように比較的単純な幾何学的距離秩序パラメータλを用いた。H2O種のみからなる系ではλは最大の共有結合O-H結合距離であり、OH-とH3O+を含む系ではλはOH-中の酸素とH3O+中の水素原子間の最短距離とされた。 以下では、秩序変数に用いた酸素原子をOλと呼ぶことにする。 水素の種類（OH-, H2O, H3O+）は、各水素に最も近い酸素と結合する単結合を割り当てることで識別した。 なお、秩序変数の定義には化学結合を定義するための閾値は必要なく、シミュレーションの期間中、秩序変数を特定の水分子に拘束するものでもない。 つまり、単一のターゲットとなるO-H結合や水分子ではなく、系内の任意の水分子の解離速度を計算することになる。
今回のケースでは、最後の界面を我々のシステムで得られる最大距離の先に置いた。 そのため、すべての軌道は、システムが再びH2O種のみを含むようになるまで伝播された。 分離されたイオンは、分離が大きくても数フェムト秒以内に再結合する可能性があり(16)、この観測は我々の解析でも確認された(SI Appendix, Fig. S1)。 準安定なイオン化状態をよりよく識別するために、我々は経路再重み付け（21）を用いて、軌道の長さ（フェムト秒）に等しい別の秩序変数λ′に交差確率を投影した
図1では、我々のシミュレーションから計算した交差確率を秩序変数の関数として示している。 原理的には、最初のプロトンジャンプの後、反応座標λを増加させる2つのメカニズムが考えられる。 イオン種がさらに別のプロトンジャンプによって分離する、いわゆるGrotthussメカニズムで、H3O+またはOH-イオンが別の酸素に再割り当てされ、反応座標が突然不連続に増加する。 もう一つの可能なメカニズムは、最初のイオン種をそのままにして、拡散によって互いに遠ざかるようにし、より緩やかな反応座標の増加をもたらすものである。 1.5 Åから3.2 Åの間の完全に平坦な中間プラトー領域に基づいて、我々は第一の機構のみが有効であると結論付けることができる。 λ>3 Åについては、λを秩序変数と考え、λ′≧1 psを安定解離事象を識別する基準として使用した。 この選択は、逆組換え反応の時間スケールが明確に分かれていないため、交差確率の平坦なプラトー領域が発生するため、かなり恣意的である。 閾値を1psとすると、交差確率はPA=4.0×10-15となる。 シミュレーションで計算した初期フラックスfA=2.9×10-3 fs-1と合わせると、解離定数はkD=fA×PA=1.1×10-2 s-1となった。 時間閾値を必要としない別の解離速度定数は、水素の入れ替わりを起こした軌道をカウントすることで定義できる。 この定義の根拠は、プロトンの交換は水素結合ネットワークの重要な再編成を意味し、逆反応を独立した再結合反応とみなすことができるためである。 逆に、順方向の反応は、相関のある逆方向の反応に追随しないので、準安定状態を確立していることになる。 この定義から、kD=0.16 s-1の速度が得られる。
25℃で実験的に決定された解離定数と比較すると、速度定数を500倍過大評価している（ただし、大きな閾値を選択するとシミュレーション速度は低下して実験速度に近くなる）。 精度に影響を与えるすべての要因（統計誤差、関数、小さなシステムサイズ、プロトンの純粋な古典的処理、時間のしきい値）を考慮すると、実験との偏差は満足できるもので、他の密度汎関数理論研究と同程度です。 第一原理計算で考慮した関数によっては、エネルギー障壁は10-20 kJ/molの誤差を持つことがあり(22)、室温ではすでに実験と理論の速度定数の55-3000倍の差に相当することになる。 それでも、密度汎関数理論は一般に実験と合理的に一致する傾向や機構情報を再現することができる(23)。
我々はまた、生成軌道の平均エネルギーを秩序変数の関数として計算した（図1）。 このエネルギーは速度定数の温度微分から導かれる活性化エネルギーに収束すると考えられる(24, 25)。 この活性化エネルギーは、秩序変数の選択に依存する自由エネルギー障壁よりも、より直接的な実験との比較を可能にすることに注目したい。 受容経路の平均エネルギーから得られる活性化エネルギーは約17.8 kcal/molである。 比較のために、Natzle and Moore (8) の実験データのアレニウスプロットでは、活性化エネルギーは約17.3 kcal/mol となり、Eigen and Maeyer (7) では活性化エネルギーは 15.5 kcal/mol と 16.5 kcal/mol の間になったと報告されています。 我々の結果との乖離は、上記の典型的な誤差よりも小さく、また、速度定数が低いにもかかわらず、実験の活性化障壁が我々のシミュレーション結果よりも低いという事実は、むしろ驚くべきことである。 4264>
パスサンプリング法は反応性（および非反応性）軌道を生成し、可能なメカニズムや開始条件を発見するために使用することができます。 これらの条件を特徴付けるために、我々はξ=(ξ1,ξ2,…)とラベル付けした追加の集団変数を考えた。 原理的には、これらのξは系内のすべての位置と運動量の関数であり、必ずしも単純な物理的解釈を持つとは限らない。 水素結合を形成する能力は水の特徴の一つであり(26)、以前の計算機研究によってイオン化した種をつなぐ水素結合線の関連性が示されているので(16、17)、我々は水素結合ネットワークと四面体形状からの歪みを定量化する一連の比較的単純な集団変数に注目した(
最初の集団変数は新生イオン種に架けられた水素結合線の長さである。 我々の目的は、開始された軌道の結果を予測することであり、特に反応的な事象の開始条件を予測することである。 したがって、水素結合ワイヤーをイオン種をつなぐものとして定義することはできない。これは、我々が予測したい結果の一つであるからだ。 一つの軌道に対して、λが与えられた閾値、λt=1.15Åより大きい最初の時点で、Oλ種とi-1個の他の水種を含む最短のワイヤーを水素結合ワイヤーと定義しました。 これは、長さwiが連続するメンバーのO-O距離の合計として得られるi個の水種を含むワイヤーを定義する。
さらに、Oλ種を取り巻く局所構造を記述する以下の4つの集団変数を考慮した。 (i) O-λを含む水種が受け入れた水素結合の数na、(ii)供与した数nd; (iii) O-λとその4つの最近接酸素原子によって定義される角度を用いて得られた四面体秩序パラメータq (27, 28)(by the definition q=1 for a perfect tetrahedral structure and otherwise q≠1); (iv) ワイヤの二つの内部角度の余弦の最小値と定義されている角度順序パラメータ、qcos; とした。 4264>
余分なξを定義した後、我々は予測力法（19）を用いて軌道を分析した。 この方法は、まず、λr>λc>λ0となるように定義した2つの閾値λrとλcに基づいて軌道を反応性か非反応性に分類するものである。 軌道は指定されたλrに到達すればreactive、そうでなければnon-reactiveとみなされる。 rλc,λr(ξ) は、λc を通過する軌道のうち、点 ξ で λc を通過し λr に到達する割合、 uλc,λr(ξ) は、λc を通過する軌道のうち、点 ξ で λc を通過するが λr には到達しない割合である。 この2つの分布は、追加された秩序変数と反応性の関係についての情報を与えます。 例えば、uλc,λr(ξ)=0の場合、ξはアクセスできないが、ξでλcを横切ることができれば、その軌道は反応性を持つということがあり得ます。 異なるξの重要性を定量化するために、予測能力TAλc,λrを計算すると、(19)TAλc.と定義される。λr=1-1PA(λr|λc)∫rλc,λr(ξ)uλc,λr(ξ)+uλc,λr(ξ)dξ,such that 1≧TAλc,λr≧PA(λr|λc). 集合変数が反応性に相関しない場合は下限に達しますが、ξが反応に関係する場合はTAλlc,λr>PA(λr|λlc)となります。 交差確率のみを用いた場合と比較して、余分なξを考慮した場合にどの程度予測能力が高まるかを、比率TAλλc,λr/PA(λr|λc)≧1としています。 なお、式2の定義から、2つの分布の重なりが小さければ予測能力は高くなることがわかる。
まず、3、4、5個の水分子を含む水素結合線の長さを調査した。 これらの集合変数（それぞれw3、w4、w5）の予測能力を比較すると、w4とw5は反応性と相関が高く、w4はより大きなλrに関連することが分かった（SI Appendix、図S2）。 そこで、以下では、4つの水分子を含むワイヤに着目する。 水のワイヤでは、イオン種が少なくとも2つの水分子によって分離されている場合、1つの水分子によって分離されている場合と比較して、イオン状態がより長い時間存続することが観察された。 これは、少なくとも3回のプロトン移動が起こったことを意味する。 そこで、最初に共有結合していたO-H結合の距離をモニターし、1回目（｜OH｜1）、2回目（｜OH｜2）、3回目（｜OH｜3）のプロトン移動について図2に示しました。 グロッタス機構(29, 30)から予想されるように、最初の自己イオン化の後、数回のプロトン移動が起こり、イオン種がワイヤーに沿って分離していくことがわかる。 図2は、この現象が、協調的にも段階的にも起こりうることを示している。 1個目と2個目のプロトンの移動はほとんど協奏的に起こるが、3個目のプロトンの移動は（起こるとしても）段階的に起こることも協奏的に起こることもある。 このことは、これらのイベント間の待ち時間にも反映されています（SI Appendix, Fig.S3）。 2回目のプロトン移動と3回目のプロトン移動の間の待ち時間分布は、1回目と2回目の移動と比較して広くなっています。 また、ワイヤーの安定性を調べるために、時間を反転させた軌道における水素結合ワイヤーの計算も行いました（SI Appendix, 図S4）。 その結果、Hassanaliら（17）が報告したように、軌道は確かに収縮したワイヤー（w4<7.6Å）で始まり、終わるが、最後にはこれらのワイヤーは必ずしも同じ酸素原子を含んでいないことが分かった。 これは、水素結合ワイヤーが実際に切断されたためか、より少ない破壊（例えば、5員環のワイヤー内で連続する4つの酸素の選択がずれるなど）により発生する可能性がある。 長い軌道の大部分は別のワイヤーを経由して再結合しますが、それでもかなりの数の長い軌道(>1 ps)があり、それらは解離経路と全く同じになります。 これは、ワイヤの切断が準安定状態に到達するための必要条件であるという仮説(16)と矛盾しています。 また、目視では、水素結合のごく短いオンオフの揺らぎを除いて、水素結合ワイヤーがそのまま残っている比較的長い軌道が存在することがわかる。 したがって、上記の仮説を擁護することは困難であると判断した。 逆に、実際に切断された場合、常に長寿命の準安定状態になるかどうかも調べることができる。 そのために、水素結合が交換された軌道はすべて、必ず水素結合の切断を意味するという仮定を再び採用しました。 SI Appendix, Fig. S5 は、プロトン交換を伴う軌道は平均して長くなるが、それでも比較的短い（35fs）ことを示している。
自己イオン化イベントの協調的な振る舞い。 (A)4つの軌道における1回目(i=1)、2回目(i=2)、3回目(i=3)のプロトン移動の際の初期共有結合O-H結合の距離(|OH |i). 矢印は時間の方向を示し、異なる軌道は異なるタイプの水素移動を例示しています：失敗したステップワイズ（濃い灰色、｜OH｜1-｜OH｜2に対してのみ表示）、協奏的（薄い灰色）、｜OH｜1-｜OH｜2に対してのみ協奏的（青色）、協奏的ステップワイズ（オレンジ色）です。 (B）最終パスアンサンブルのすべての軌道を使用して、｜OH｜1-｜OH｜2と｜OH｜1-｜OH｜3の密度を求めた（左の列）。 右の列は、長さtpath>60fsの軌道を考慮した場合の密度である。 4264>
追加の集団変数（SI Appendix, Figs S6 and S7）を比較すると、ndは他の変数よりも関連性が低いことがわかり、これ以上考慮しないことにした。 他の集合変数は反応性との相関が高く、Fig.3Aではそのいくつかの組み合わせで予測能力を示している。 Fig.3Bでは、TAλc,λrをλc=1.16Åにおけるλr≤2Åの関数として、集合変数のいくつかの組み合わせによる交差確率と比較している。 Fig. 3のA,Bは、交差確率と比較して、予測能力を107倍高めることができることを示しています。 この場合、交差確率はTAλc,λr∼0.4と小さいので、完全に予測することはできないことに注意してください。 このことは、Geisslerら(16)が示唆したように、記述に重要な他の集団変数が存在することを示している。 図1Cから明らかなように、λ=2.0 Åに達する軌道の大部分は、長寿命の準安定状態に至らないことがわかります。 最初のスナップショットからこれを予測するのは、最初に伸びたOH結合から遠く離れた多くのMDステップの後の水分子間の衝突に依存するので、まだ一歩が足りないように思われます
付加的な集団変数を考慮することによって、水の自己イオン化の予測力を高めることができる。 (A)水素結合長（ξ=w4）、配向秩序パラメータ（ξ=q）、角度秩序パラメータ（ξ=qcos）、Oλ種が受け入れる水素結合数（ξ=na）という追加の集団変数を用いた交差確率に対する予測力（TAλc、λr）。 (B) λc=1.16 Åと集団変数の異なる組み合わせに対する予測力と交差確率をλrの関数として表したもの。 4264>
開始条件をより詳細に調べるため、図中の反応・非反応分布rλc,λr（ξ）、uλc,λr（ξ）を調査したところ、λcは1.15ÅでPA=1になっている。 図4はλc=1.16Å, λr=2.0Åのときのもので、ここでは早く再結合するものも含めてすべての解離事象を調べ、ξ=(w4,na)の分布を図4Aに示している。 これは、λc=1.16Åを横切る軌道が、短いワイヤー（w4が小さい）ほど反応する確率が高いことを示しています。 これは、Hassanaliら(17)が最初に示唆したように、自己イオン化の重要な条件として「圧縮された」ワイヤーを仮定していることを裏付けている。 このことは、Hassanaliら(17)が提唱した、Oλが過配置されたワイヤーが反応性を持つ確率が高いという仮説を支持するものである。 それでも、Fig.4A のどの点（w4,na）でも反応しない確率が大きくなっています。 例えば、(i) 7.15<w4<7.6 で、同時に na=3 であれば、反応性がある確率は 3.6⋅10-6 であり、小さいがそれでも λc のランダム点から反応性がある確率より 58 倍も大きいことがわかる。 より極端な場合、(ii)w4<7.3と同時にna=4とすると、チャンスは0.15に増加する。 予測能力TAλλc,λrは、これらのチャンスの加重平均を提供し、その重みは関連性に比例する(19)。すべての反応軌道のうち、領域iで45％がλcを横切り、領域iiでは0.6％だけなので、後者は75倍低い重みとなる
開始条件と局所的な集団変数。 (A) ξ={w4,na}, λc=1.16Å, λr=2.0Å における反応性（rλc,λr(ξ))、非反応性（uλc,λr(ξ))の分布。 上と右の挿入図は、分布の一次元投影図を示しています。 2つの分布はw4座標に沿って明確に分離しており、反応性の軌道は非反応性の軌道に比べてより圧縮されていることがわかる。 また、秩序変数の計算に用いた酸素原子（Oλ）は、反応性軌道では非反応性軌道に比べて平均してより多くの水素結合を受け入れていることがわかる。 (B）反応性軌道のスナップショット（Oλは青色で表示）。 四面体秩序パラメータqの計算に使用される4つの酸素原子はオレンジ色で示されています。 水線は黄色の線（と灰色の透明な球）でハイライトされ、角度パラメータqcosが示されています。 このスナップショットでは、水ワイヤーは圧縮され、qは四面体構造からのずれを示し、qcosは3つの酸素原子がワイヤーに並んでいることを示し、Oλは3つの水素結合を受け入れ、1つを提供する（緑の線で示す）ことがわかります。
q座標を考えると、rλc,λrはuλc,λrに比べてq値が低い方にシフトしており、解離する水種の周りの四面体配置からの変形もイベントの起点になる可能性があることがわかります。 この発見は、他のいくつかの水相化学反応では逆の効果が見いだされたので、やや意外である(31)。 同様の結論は、ξ=(w4,qcos)の分布についても導き出せる。 ここでは、水分子の線形配列に近い反応性分布に対して、qcos座標に沿ったピークが存在する。 図4Bは、反応軌跡の初期（3fs後）に得られた代表的なスナップショットである。 図3の結果から、水ワイヤの圧縮（w4で測定）、過配置（naで測定）または歪み（qとqcosで測定）が自己イオン化に必要な開始条件であることが分かりました。 しかし、図3BのTAλc,λrの値が示すように、これらは十分な条件ではない。理想的なξパラメータ範囲内で出発した軌道の60%は、まだプロトンジャンプを協調的に確立できていない。
パスサンプリングデータに機械学習（ML）を適用すると、人間の直感で見落としやすい重要な集団変数を見出せる可能性がある (33,34) 。 この可能性を探るため、軌跡の初期に水系の状態を与えて軌跡の結果を予測するMLモデルを構築した。 予測力の解析と同じ範囲に着目し、λ>1.15Åに初めて到達したときの系の状態を用いて結果を予測する。 我々は、すべての奇数パスアンサンブルをキャリブレーションに含め、偶数パスアンサンブルをテストセットに使用するいくつかのML技術を使用した。 各パスアンサンブル内のデータを均等に2分割する代替分割でも、同様の結果が得られた。 また、大きく歪んだ分布はMLで扱いにくいため、対応するパスアンサンブルの統計的重みによるデータセットの再重み付けはさらに省略した。 しかし、予測力法(19)の中で定量的にテストできる新しいパラメータを見つけるために、定性的アプローチとしてML技術を適用した。
さらに、過剰解釈の潜在的リスクを避けるために、ML決定プロセスの複雑さを制限することを選び、TAλc、λrの計算時に最大4次までのパラメータを課した。 例えば、アンサンブルベースのgradient-boostingマシンを用いて優れた予測性能(>90%)が得られた(35, 36)。 しかし、100～150本の深い決定木のアンサンブル（順番に追加される）を用いるため、モデルの解釈に問題がある。 性能は向上するが、偶発的な相関を持つオーバーフィッティングの可能性が高まる。 そこで、分類回帰決定木（CART）（20）に基づく単一木ベースの決定モデルに限定した。 TAλc,λr関数の次数パラメータを4つに制限したのも同様の理由によるものである。 パラメータを増やすと、反応性/非反応性分布を表す行列がより疎になり、その結果、これらの分布間のオーバーラップを計算する数値積分は、ビンサイズに非常に敏感になり、不十分な統計でビンが空になることによりオーバーラップを過小評価する可能性があります。
我々は、酸素-酸素距離、最初に結合した水分子の酸素-水素距離、Oλとその4近傍酸素が形成するすべての角度、オーダー3、4、6のスタインハート秩序パラメータ（32）からなる138の集団変数を検討した（詳細は材料と方法を参照）。 さらに、すでに考慮されている次数パラメータを追加した。 Fig. 5Aは、その結果得られた決定木である。 驚くべきことに、すべての入力パラメータのうち、w4パラメータは決定木の頂点にあり、決定木の各分岐において各変数に起因する分類誤差の減少によって測定される最も重要な変数である(20) (SI Appendix, Fig. S9)。 また、四面体秩序と受容された水素結合の数も決定木に登場する。 最初の効果を説明するために、MLアプローチはSteinhardt q4順序パラメータを、我々が以前使用した同様のqパラメータよりも優先させた。 Oλとその25番目に近い酸素との距離であるd25のように、決定木にも現れるいくつかの距離は、データセットのサイズが限られているため、偶然の相関関係によるものである可能性が最も高い。 このことは、この変数の重要度を調べることによって検証される：d25は、最も重要な20の変数の中に現れず（SI Appendix、図S9）、実際、他の同様の変数（例えば、d24）は、重要度は低いものの、より上位にランクされている。 ML 法で示唆されたより重要で直感的なパラメータは、Oλ に最も近い酸素とその分子内結合が最も大きい水素との間の OH 距離である λ2 である。 MLツリーのパラメータを用いた再計算（図5B）では、w4、q、na、qcosの組み合わせより高い性能は得られなかったが、統計的不確実性を考慮すると同等と考えられる
機械学習による解析結果。 (A）開始された軌道の結果を予測するための分類と回帰の木。 ここでは、いくつかの追加の集団変数（材料と方法に記述）を考慮したが、木を構築するために最終的に必要とされるのは小さなサブセットだけである：w4、q4、λ2（Oλ種に最も近い水分子の伸長した水素結合の長さ）、di（Oλからi番目に近い酸素への距離）およびd¯i（i最も近い酸素を考慮した平均距離）。 ノードの表記は、左上の単体ノードを用いて説明する。 この木は、λが最初に≧1.15となる軌道のフレームで得られた集団変数に基づいて、軌道が反応性、すなわち、λ≧2に達するか、非反応性であるかを予測するものである。 反応性の軌道を予測するノードは青色（クラス1）、非反応性の軌道を予測するノードは緑色（クラス0）に着色されています。 なお、パスアンサンブルは統計的な重みで再重量化されていないため、四角の底のパーセントは物理的に正しい割合を反映していない。 例えば、ルール 5 (反応性の軌道を予測) は w4≥7.6 と λ2≥1.1 のように表現できる。 これらのルールは異なる開始条件を与え、SI Appendix, Table S1 に最下段のノードについて記載されている。 (B) λc=1.16 Åと集団変数の異なる組み合わせにおけるλrの関数としての予測力および交差確率。 ここでは、我々が同定した集合変数と機械学習解析で重要とされた変数を用いて予測力を比較している。 (C) ξ={λ2,d¯2}, λ=1.16 Å, λr=2.0 Åにおける反応性（rλc,λr(ξ))と非反応性（uλc,λr(ξ))の分布です。 上と右の挿入図は、分布の一次元投影図である。