Closing Kok’s cycle of nature’s water oxidation catalysis

The Mn4CaO5(6) cluster in photosystem II catalyzes water splitting through the Si state cycle (i = 0–4). Molecular O2 is formed and the natural catalyst is reset during the final S3 → (S4) → S0 transition. Only recently experimental breakthroughs have emerged for this transition but without explicit information on the S0-state reconstitution, thus the progression after O2 release remains elusive. In this report, our molecular dynamics simulations combined with density functional calculations suggest a likely missing link for closing the cycle, i.e., restoring the first catalytic state. Specifically, the formation of closed-cubane intermediates with all hexa-coordinate Mn is observed, which would undergo proton release, water dissociation, and ligand transfer to produce the open-cubane structure of the S0 state. Thereby, we theoretically identify the previously unknown structural isomerism in the S0 state that acts as the origin of the proposed structural flexibility prevailing in the cycle, which may be functionally important for nature’s water oxidation catalysis.

The Mn 4 CaO 5 (6) cluster in photosystem II catalyzes water splitting through the S i state cycle (i = 0-4).Molecular O 2 is formed and the natural catalyst is reset during the final S 3 → (S 4 ) → S 0 transition.Only recently experimental breakthroughs have emerged for this transition but without explicit information on the S 0 -state reconstitution, thus the progression after O 2 release remains elusive.In this report, our molecular dynamics simulations combined with density functional calculations suggest a likely missing link for closing the cycle, i.e., restoring the first catalytic state.Specifically, the formation of closed-cubane intermediates with all hexa-coordinate Mn is observed, which would undergo proton release, water dissociation, and ligand transfer to produce the open-cubane structure of the S 0 state.Thereby, we theoretically identify the previously unknown structural isomerism in the S 0 state that acts as the origin of the proposed structural flexibility prevailing in the cycle, which may be functionally important for nature's water oxidation catalysis.
For a long time, stories taking place during the crucial S 3 → (S 4 ) → S 0 transition had been rarely accessible to experimentalists.Very recently, remarkable experimental breakthroughs on this transition were reported by Bhowmick et al. 8 and Greife et al. 21, in which major clues for crystallographic structures and reaction kinetics have been provided.While "S 3 → S 4 " for generating the Mn(IV)-O• radical by single-electron, multi-proton transfer is identified as the kinetic bottleneck of "S 3 → S 0 " (ca.2.5 ms), the subsequent O 2 formation and release during "S 4 → S 0 " are claimed to be much faster 21 .Although debates existed [22][23][24][25][26] , the coupling between O5 and Ox (shown in Fig. 1) is regarded as the most viable mechanism for O-O bond formation in both reports.However, the cluster evolution after O 2 release remains unclear with respect to the possible intermediates involved and the atomic-level details of the transformations leading to reconstitution of the S 0 state, i.e., closing Kok's cycle.Consequently, the essential issues remaining in the recovery of the natural catalyst to its first catalytic state call for followup research.
Computational chemistry plays a prominent role in proposing realistic models for extremely fast processes that are difficult to be probed by experiments.Born-Oppenheimer ab initio molecular dynamics (BO-AIMD) simulations combined with minimum energy path (MEP) searches at the level of density functional theory (DFT) are employed in this work.These computations are logically based on the putative "oxo-oxyl coupling" mechanism in the S 4 state 11,21,27 , at present the leading and most widely accepted proposal for O 2 formation (see Supplementary Note 1 for more discussion).Actually, the simulation results would be generally applicable in whatever coupling ways as long as O5 and Ox are substrates, which is most favored by the recent substrate water exchange experiments 28,29 .The main purpose is to investigate how the OEC cluster with a structural cavity left by release of O 2 (Im0 ÀO 2 , Fig. 2) would evolve step by step, and thereby explore feasible pathway to the S 0 state.The study aims to shed light on the nature and origin of substrate water and provide insights on the underlying molecular mechanism resulting in the reconstruction of the Mn 4 CaO 5 cluster, for a more complete understanding of nature's water oxidation catalysis.
To characterize the conformational changes of the cluster during water insertion in detail, our BO-AIMD simulations, without application of any steering force, execute a long simulation time based on a large model using full quantum mechanical treatment of DFT (Supplementary Fig. 1).BO-AIMD allows a priori exploration of the potential energy surface of a low-barrier transition through non-biased sampling, and hence identification of not already presumed reaction pathways and intermediates.The dynamic trajectories are strictly confined along the ground state under adiabatic approximation.It should be highlighted that the simulations cover only part of the "S 4 → S 0 " stage after O 2 release, instead of the "S 3 → S 0 " transition in a millisecond timescale that would not be captured in a picosecond trajectory (see more clarification in Supplementary Note 2).For further chemical reactions that may take place at timescales longer than tens of picoseconds of a typical AIMD simulation, MEP searches for specific reactions were carried out by truncated DFT models.The obtained results are beyond the current knowledge and discussed below.

Water insertion dynamics and Mn 4 Ca cluster evolution
Previous theoretical studies involving the S 4 → S 0 transition have reached a general consensus that it is W3(H 2 O), rather than W2 or any other crystal water molecule, that would refill the vacant site formed by O 2 release [30][31][32][33] .Since stepwise occurrence of O 2 release and water insertion is verified more favorable than the concerted mechanism 30,33 , Im0 ÀO 2 is validated as the starting state pending water insertion (see more detailed discussion in Supplemenatry Note 4).Our BO-AIMD simulations for both octet/αααβ and doublet/αβαβ spin states (see Supplementary Note 5 for spin state definition) have observed interesting phenomena regarding the structural evolution of the cluster that have not been reported before.On that basis, water insertion can be divided into three basic stages, represented by Im0 ÀO 2 and the other two intermediates, Im1 and Im2 (Fig. 2).For the initial protonation state of W2, we here follow the hydroxide form W2(OH − ) that was suggested to be more consistent with the experimental magnetic and electron paramagnetic resonance (EPR) spectroscopic data 34,35 and the pKa predictions by electrostatic energy computations 36 .However, it is noted that other studies that attempted to reproduce other types of spectroscopies such as Fourier transform infrared (FTIR) 37 and X-ray absorption spectroscopy (XAS) 38 , and calculations of proton hyperfine coupling constants 39 favor the doubly protonated form W2(H 2 O) in certain S states.Therefore, the nature of W2 remains an open question (see Supplementary Note 3 for more discussion).
Around the initial 5 ps of the simulations, the Ca-bound W3(H 2 O) was found to insert into the structural cavity of Im0 ÀO 2 , however, it would not directly migrate to the bridging position between Mn3 and Mn4, but moves towards Mn1 instead.Thereby, W3 forms a short, strong hydrogen-bonding (HB) interaction with W2(OH − ), while being in the bridging position between Ca and Mn1.This causes the squarepyramidal to trigonal-bipyramidal conversion of the Mn4(III) coordination, similar to the proposed five-coordinate Mn4(IV) in the local geometry 40,41 .Alongside, simultaneous movements of W5(HOH605) and W6(HOH577) were observed, occupying the original locations of W3 and W5 in close proximity, respectively (arrows in Im0 ÀO 2 ), and in this way, Im1 is formed.In the next 5-10 ps, the spontaneous proton  12,21,30,34,64,90 .
transfer from W3(H 2 O) to W2(OH − ) occurs, which was found to be highly correlated (almost synchronous) to W7(HOH523) binding to Mn4, and W8(HOH529) is also gradually pulled along that way (arrows in Im1).These structural changes are accompanied by a "pivot/carousel"-like reorganization 40,42,43 of the Mn4 ligands and expansion of the cluster with elongated Mn1-Mn4 distance (see Supplementary Note 6 for more analysis on the structural changes and energetics); at the same time, the deprotonated W3(OH − ) approaches closer to Mn3 arriving at a bonding distance.As a result, Im2 is formed, which has a typical closed-cubane structure with a saturated octahedral Mn4 coordination and W3 as the Ca/Mn1/Mn3 μ 3 -OH bridge group, and this conformation is dynamically stable thereafter.All these events were spontaneously take place during the simulations, indicating barrierless (or almost) events that can be easily captured in the dynamic sampling within tens of picoseconds (see Supplementary Movie 1).
The pathway observed for water insertion and structural transformation makes sense on the basis of molecular principles.The early formation of a closed-cubane (termed "B") instead of an open-cubane structure (termed "A"), i.e., that W3 initially binds to Mn1 rather than Mn4, can be expected for two reasons.Firstly, compared with Mn4, Mn1 is placed at a shorter spatial distance and with closer bond connectivity to Ca 2+ (an indispensable cofactor for charge compensation in the cluster [44][45][46][47] ), making Mn1 more positively charged than Mn4 (Supplementary Figs. 2 and 3) and thus a better Lewis acid for W3 coordination.Secondly, W2 was found to rotate moderately towards the cavity in the early phase of the simulations, which impedes W3 ligating to Mn4 but favors W3 binding to Mn1 by forming a strong HB interaction with it, as shown in Im1.Next, closing the Mn 3 CaO 4 cubane is fulfilled by Mn3-W3 bonding, which necessitates W3 deprotonation (to W2) due to the stronger Lewis basicity of the deprotonated W3(OH − ) bonding in the μ 3 -position.Meanwhile, W3 movement to Mn3 further promotes the "pivot/carousel"-like rotation of W2 (in HB interaction with W3) together with W1 around the Mn4 axis, creating a vacant site for W7 coordination from the O4 channel, which is a possible water delivery system 8,48-52 .Thus, it can be seen that the ligand rearrangements on Mn4 and the strongly coupled W7 binding with W2 protonation are attributed to the electrostatic affinity between Mn1/Mn3 and W3.Concurrently, these observations illustrate the critical roles of both Ca and Mn4 in water transport to the core position of the cluster.It is surmised that motions of outer crystal waters in HB network would also be involved, which is, however, not possible to be observed in a finite model (see Supplementary Note 7 for more discussion).Variations of the key interatomic distances with time evolution are displayed in Fig. 3.The electronic configuration of the cluster and metal oxidation states Mn1(III) Mn2(IV) Mn3(III) Mn4(III) remain unchanged throughout; see Mulliken spin populations in Supplementary Figs. 4 and 5.The almost indistinguishable phenomena for the two spin states reflect the insignificant role of the ferromagnetic/ antiferromagnetic couplings between Mn2 and the other Mn centers.Furthermore, the commonality underlying the two simulations represents a mutual testification, eliminating the randomness and indicating the reproducibility of the results, and justifying the reliability of the conclusion.
With respect to the energetics of the three intermediates, additional geometry optimizations were applied to the snapshots extracted from the BO-AIMD simulations, showing an obvious downhill process as much as ca.−30 kcal mol −1 in free energy difference from Im0 ÀO 2 to Im2 (Supplementary Table 1).However, it is inadequate to draw direct correlations in terms of the energetics to experimental measurements based on the present model, for which other parts and events during the S 4 → S 0 transition (e.g., O 2 release and H + discharge leading to HB network rearrangement coupled to protein dynamics) beyond the local structural evolution surrounding the Mn cluster are not explicitly covered.Thus it is emphasized that the computed energetics do not fully represent the complete donor-side reactions of PSII 53 , but only the intrinsic/local energetics of the conformational Mn3, and W7 binding to Mn4; the green arrow denotes H + is suggested to be released to the lumen from W1 via Asp61.The routes of the key atomic motions observed during the simulations are marked with blue arrows.The red, yellow, and blue background colors stand for the three basic phases of Im0 ÀO 2 , Im1 and Im2, respectively, which correspond to those in Fig. 3. changes occurring at/around the inorganic cluster, and therefore the validity of the structural intermediates of the cluster from Im0 ÀO 2 to Im2 remain unaffected (see Supplementary Note 8 for more discussion).Anyway, the results demonstrate the thermodynamic rationality and kinetic feasibility for the ultrafast water insertion and cluster evolution after O 2 release.Thus, it is suggested here that W3 approaching Mn1 followed by formation of a closed-cubane structure, which is shown to occur at a picosecond timescale, is identified as the most favored pathway for water insertion under the present study.W1, now trans to μ-O4 in Im2, can easily deprotonate to the lumen via D1-Asp61 (arrows in Im2) seen from the obtained reaction energetics (Supplementary Note 9, Supplementary Table 2, Supplementary Figs. 8, 14 and 15).This is expected to be the second released H + during the S 3 → S 0 transition, producing a state defined as "pre-S 0 " herein, which will be studied further in the subsequent section.

Attainability of the open-cubane S 0 state
The above change of Im0  3 and 4) by MEP calculations based on the evolved structures from the BO-AIMD simulations (Supplementary Figs. 9, 10, 16 and 19).As depicted in Fig. 4, W2 decoupling from Mn4 is estimated to be slightly endothermic by ca.2-3 kcal mol −1 with a small transition state (TS) barrier of ca.4-5 kcal mol −1 , a magnitude that can be easily overcome by thermal vibrations in the protein matrix.According to indications from previous DFT studies in transition-metal chemistry [54][55][56][57] , as well as our sensitivity test on functionals (see below), the calculations in this case can give quite reliable results with possible errors normally within a range of a few kcal mol −1 that would not affect the feasibility of the reaction.Ligand dissociations are ubiquitous in organometallic chemistry, mostly acting as pre-steps in substitution reactions for formation of catalytically active species as necessary intermediates 58   similar to a recent report on a biomimetic polyoxometalate water oxidation catalyst 59 , and the reactivity is enabled by the presence of Jahn-Teller (J-T) effect at Mn(III) which extends the bonding distance of the J-T axial ligand and thereby facilitates its de-coordination.
Besides, the coincidence of W2 protonation and W7 binding synergistically further weakens the W2 coordination to Mn4, because of the further elongated Mn4-W2 bond and the "structural trans effect" 60,61 in octahedral transition-metal complexes.Formation of the pre-S 0 (W2-unbound) state, with W2 moving out to a distance about 4 Å away from Mn4, is analogous to the last necessary step required in the substrate water exchange mechanism where a water molecule dissociates from Mn4(III) in a closed-cubane structure (either exoergic (S 1 ) or endoergic (S 2 )) 62 .Furthermore, the possibility of W2 decoupling from Mn4(III) in a closed-cubane structure is also reflected in a hypothetical mechanism for water exchange in the S 0 state 63 .Thus we consider the dissociation of W2(H 2 O) at this state as a chemically sensible process (see Supplementary Note 10 for more discussion).In this way, W2 is released to the bulk and the S 0 B (closed-cubane) state is formed.
Thereby, Mn4(III) in the resulting S 0 B structure contains an unsaturated peta-coordination sphere which allows for the interaction with μ 3 -W3(OH − ).As shown in Fig. 5, the μ 3 -W3(OH − ) shift from Mn1 to Mn4 can be realized across a low barrier of ca.4-6 kcal mol −1 , reaching the almost isoenergetic S 0 A structure that is perceived as the generally accepted S 0 state structure 4,30,34,64 , and the intermetallic distances (2.77, 2.76, 2.84 and 3.40 Å) are basically consistent with the structural constraints from the extended X-ray absorption fine structure (EXAFS) experiments 65,66 .According to the Eyring-Polanyi equation, the reactions should occur at a timescale of nanoseconds at room temperature, which is within (far below) the experimental limitation of ca.2.5 ms observed for the S 3 → S 0 transition 21 .Judging from thermodynamics and kinetics, the isomerization is readily reversible, switching the orientation of the J-T elongation axis of Mn4(III) along the W7-W3 or Asp170-Glu333 vector, while retaining those of Mn1(III) and Mn3(III) perpendicular to each other.Therefore, the conformational interconversion in S 0 , with regard to mechanism, resembles that of the S 1 state induced by the J-T effect 67 and fundamentally differs from that of the S 2 state where valence isomerism (III) ↔ (IV) between Mn1 and Mn4 is involved 68 .Moreover, the above findings have been exposed to a sensitivity measurement where the outcomes by different DFT functionals have been examined using a series of different dispersion-parameterized hybrid and nonhybrid functionals (Supplementary Note 11, Supplementary Table 5 and Supplementary Fig. 11).In addition, a similar test has been performed to verify the structural isomerism using another structural model originating from the 3 F XFEL data (PDB ID: 6DHP) 4 with the S 0 state as the major population (Supplementary Table 6 and Supplementary Figs. 12 and 13).The obtained thermodynamic parameters are generally in accordance with that of Fig. 3, even though slight variations depending on model/functional are as expected observed (deviations within 1-2 kcal mol −1 ), and such minor deviations would normally be expected given the inherent approximations and limitations of the DFT methodology [54][55][56][57] .Although it is not possible to offer an absolute quantification of the reaction energetics, a reliable qualitative conclusion can be safely drawn, i.e., W2 dissociation is achievable, and S 0 A and S 0 B are quasi-isoenergetic and interconvertible through a low barrier.Thus we emphasize the structural isomerism as a basic function of the OEC that is already manifested in the first state of the catalytic cycle.

Discussion
Bhowmick et al. 8 and Greife et al. 21have unveiled key structural and kinetic data during the millisecond S 3 → S 0 transition.According to Bhowmick et al., from 1200 to 4000 μs several structural changes indicative of O 2 release and/or water insertion occur, and the extended timescale between 2000 and 4000 μs may be due to pronounced variation in water positions (e.g., the slow reappearance of W20 in the O4 channel) and significant rearrangement of the HB network (including water-water, protein-protein, and water-protein interactions) related to the last proton release.This precisely corresponds to the suggested H + transfer from Im2 to the lumen via D1-Asp61 in our scheme, but the further effect on the broad water/protein environment in PSII cannot be embodied in the present study.This may provide clues for capturing more possible intermediates during the extended timescale for crystallographic snapshot data in future.As no structural evidence was observed for an empty O5 site, it is speculated that refilling of the cavity by water binding is ultrafast, which is also reflected in our simulation results.Within this period, 1200 and 2000 μs are the two essential timepoints that are closely related to (but not entirely covered by) our work.The 1200 μs snapshot signifies the onset of O 2 evolution, and the 2000 μs snapshot, without Ox on the electron density omit map, indicates completion of binding of a water molecule that refills the vacant site formed by O 2 release.On this basis, the series of processes we suggest in the present study should in principle transiently reside between the two timepoints, but represent a very short phase seen from the picosecond water insertion and subsequent nanosecond isomerization of the cluster (see Supplementary Note 16 for more discussion).This complements the details en route to the S 0 state, although the proposed temporarily present species may not accumulate to sufficiently large amounts to be detected experimentally by XFEL crystallography with microsecond intervals for snapshots.Since we show above that the S 0 B → S 0 A closedto-open conversion in nanoseconds appears as the kinetically most demanding step after O 2 release, the whole cluster transformation from Im0 ÀO 2 to S 0 A also satisfies the requirement that the rate-limiting step of S 3 → S 0 in millisecond timescale belongs to the S 3 'Y Z + to S 4 transition prior to O-O bond formation (rather than S 4 → S 0 ), as determined by Greife et al.As the specific locations of water molecules and HB interactions play a very important role for facilitating the critical step of Mn(IV)-O• formation during the S 3 → S 4 transition, they are also crucial for the resetting process of the Mn 4 CaO 5 cluster after O 2 release, especially for the W3-W2 interaction partially accounting for W3 binding to Mn1 and the internal proton transfer further promoting the Mn4 ligand reorganization (Supplementary Note 6).Besides, the HB interaction of W3-W5, of W5-W6, and W7-W8 (and similarly others outside the present model) are also indispensable for recovery of the microenvironment of the OEC cluster.Different from the singleelectron multi-proton transfer with a moderate energetic barrier (13.6 kcal mol −1 ) for the S 3 'Y Z + to S 4 transition, the present case does not involve any apparent electron transfer (no electron-hole to be reduced) and involves only one barrierless (or almost) proton transfer from W3(H 2 O) to W2(OH − ) on the picosecond timescale, which appears much less demanding than the deprotonation of Ox-H.This is also consistent with the fast kinetics for the S 4 → S 0 transition post to O-O bond formation (see Supplementary Note 17 for more discussion).
Previous to the present study, Capone et al. 33 have conducted a molecular dynamic study on the mechanism of oxygen evolution and Mn 4 CaO 5 cluster restoration, based on the oxo-oxyl coupling mechanism for O-O bond formation.The major conclusions include confirmation of the Ca-bound W3 as the inserted water molecule and validation of the two-step mechanism for O 2 release and water insertion, which are both in line with this study.For restoration of the Mn 4 CaO 5 cluster, they assume that W3 (along with deprotonation to W2) moves to the bridge position of Mn3/Mn4, directly resulting in the open-cubane structure of the S 0 state, which coincides with the schemes by refs.30-32 in terms of the binding sites of W3 throughout.However, here we propose a different pathway for how W3 should enter the cavity, i.e., it binds first to Mn1 for formation of the closedcubane intermediates followed by W2 protonation, which triggers the ligand reorganization on Mn4 and the coordination of W7, and ultimately the open-cubane S 0 state forms after W2(H 2 O) dissociation and then W3(OH − ) transfer to Mn4.While uncertainty remains regarding the cause of the differences, the two plausible pathways should be noteworthy currently and further comparative studies may be needed.Anyway, since the target structure is the same, the alternative pathway presented here reveals possible intermediates and the important structural isomerism for the first state of the catalytic cycle (see Supplementary Note 18 for more discussion).
It is noted that the XFEL crystallography resolved only opencubane structures 4,8,69 and only multiline EPR signal is exhibited 64,[70][71][72][73][74][75] for the metastable S 0 state.Besides, the 55 Mn hyperfine, nuclear quadrupole interaction, and EXAFS parameters are best matched by those estimated for the S 0 A conformation 65,66 .However, these do not conflict with the presence of the structural isomerism.Assuming that decontamination is sufficiently complete and the resolution is sufficiently high for the targeted S 0 state, one straightforward interpretation would be that S 0 A dominates over S 0 B in the XFEL samples.According to the relationship between ΔG°and equilibrium constant K, a very minor energy advantage of one isomer over the other (that is even within the DFT error range) would lead to its overwhelming proportion; see Supplementary Note 11 for an analysis on relative populations and estimated energetics.Furthermore, Cox et al. pointed out that it is possible that a certain state that can exist in multiple conformers may not show all these forms under the experimental conditions 76 .Ibrahim et al. stated that the implicit form (if present) may be short-lived due to fast formation and decay kinetics, and/or its fraction may be below the detection limit 6 .Actually, the theoretically proposed S 0 structural isomerism is experimentally suggestive in several aspects.The S 0 multiline signal in spinach is only visible in presence of a few percent methanol [70][71][72][73][74][75] , which is a strong indication for an equilibrium between at least two states that are close in energy.The suggestion is also supported by the fact that in thermophilic cyanobacterium T. vestitus the S 0 multiline signal can be observed also in absence of methanol 74 .Besides, the closed-cubane S 0 B conformation is entailed in elucidating the water exchange mechanism 63 , similar to that of the other S states 16,28,62,77 .Thus the presence of the S 0 isomerism makes sense even without a crystallographic structure available for a closed-cubane form.The situation is reminiscent of a recent experimental (EPR) evidence presented by Kosaki and Mino 78 (in line with Saito et al. 79 but in contrast to Barchenko and O'Malley 80 ) identifying the g = 4.1 S = 5/2 high-spin S 2 state as a strong support for closedcubane conformation, which is, though, also still unidentified in the XFEL studies.
Although not yet resolved by XFEL crystallography [3][4][5][6][7][8]81,82 , closedcubane structures could silently play a vital role in the catalytic cycle; see Supplementary Note 13 for additional discussion on the validity of this hypothesis. As ilustrated in Fig. 6, our theoretical simulations suggest a likely missing link during the S 4 → S 0 transition in which closed-cubane intermediates are involved, and furthermore the isomerism in the S 0 state, in addition to those already well established for the S 1 and S 2 states 67,68 .Hereby, the possibly final knowledge gap of the structural flexibility exhibited by the Mn 4 CaO 5 (6) cluster is filled, thus serving as the open-closed structural precursors for the following S-states.Since the transformation is reversible, the isomerism provides a basic reference for a mechanistic proposal of O-O bond formation in either open-or closed-cubane structure, if they are both catalytically relevant.It is inferred that different from most artificial compounds, the exclusive ability of geometric rearrangement should be a unique characteristic of the tetranuclear manganese OEC cluster, and may be associated with its inherent excellent catalytic efficiency (the correlation to be further studied).The theoretical discovery of such potential structural variants can motivate new designs and improvement of bioinspired water-splitting catalysts.
In brief, we have theoretically explored the chemistry in the S 4 → S 0 transition after O 2 release and provided important information for the S 0 -state reconstitution.The identified pathway for the insertion of the Ca-bound water molecule W3 into the cavity of the OEC cluster involves closed-cubane intermediates and is energetically favorable.This is accompanied with spontaneous W3(H 2 O) deprotonation to W2(OH − ) and a series of water molecule displacements.The subsequent W2(H 2 O) dissociation and W3(OH − ) shift to Mn4 are facile, making the open-cubane S 0 state for the next cycle attainable.More importantly, the results from this study encourage us to propose the reversible isomerization in the S 0 state that is concomitant with alteration of the J-T distorted axis of the dangler Mn4(III).The isomerism, already existing in the first state of Kok's cycle, lays the structural foundation for the subsequent S-states and may contribute to water oxidation catalysis in PSII.

BO-AIMD simulations
The BO-AIMD model for water insertion after O 2 release (depletion of O5 and Ox for the O 2 -released state) are based on the roomtemperature serial femtosecond crystallography of the S 3 state, which was taken from the second flash (200 ms) data by Ibrahim et al. 6 (PDB ID: 6W1V, monomer A) after removal of the mixed S 2 state in minor population.The model consists of the inorganic Mn 4 CaO 4 cluster, 20 amino acid residues Asp61, Asn87, Tyr161, Gln165, Ser169, Asp170, Asn181, Val185, Glu189, His190, Asn298, Lys317, His332, Glu333, Ala336, His337, Asp342, Ala344, Glu354, and Arg357, and 24 crystal water molecules HOH515(W1), HOH530(W2), HOH584(W3), HOH550(W4), HOH605(W5), HOH577(W6), HOH523(W7), HOH529(W8), HOH525, HOH596, HOH511, HOH626, HOH627, HOH522, HOH534, HOH600, HOH514, HOH517, HOH505, HOH574, HOH541, HOH508, HOH526, HOH545 and one chloride ion (Cl − 407), resulting in total 369 atoms and a net total charge of +1 (Supplementary Fig. 1).The quantum mechanical (QM) size is by far the largest one among all molecular dynamic studies on the OEC, and the specially customed GPU acceleration enables a long simulation time on such a large QM model.Since the focus of the computational study is within or closely around the OEC cluster, the large model in full QM treatment is capable of representing the conformational changes of the cluster affected by water insertion.Protonation states were chosen according to the most widely accepted scheme 18,21,34,62,83 .The O 2 -depleted model of the S 4 state, which is derived from the XFEL structure of the S 3 state by removing H + and e − from Ox, renders the initial coordinates for the Im0 ÀO 2 intermediate.The models for the energetic comparisons among Im0 ÀO 2 , Im1 and Im2 were extracted from the starting structure of the first phase, the last snapshot for the second phase, and the last snapshot for the third phase shown in Fig. 3, which were then fully optimized in the same size as used in the BO-AIMD simulations.The dispersion-corrected density functional theory, DFT-D3, using the hybrid functional B3LYP (in its standard form) was performed at double precision with the hybrid DIIS/A-DIIS scheme.The double-ζ effective core potential (ECP) basis set LanL2DZ was applied for the metal atoms Mn, Ca and the mixed full electron basis set 6-31 G*/ 3-21 G was applied for the H, C, N, O, Cl atoms (3-21 G is only used for the alkyl groups of the peripheral residues that are non-bonding to the cluster).Energy minimizations using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method were implemented before the BO-AIMD simulations at the same level of theory.Backbone constraints by fixing α-carbons of peptide bonds were applied throughout the computations.The model is placed in a Closed-cubane (B) Open-cubane (A) in the S 3 state denotes the unidirectional isomerization because the open-cubane structure is largely stabilized 12,91,92 ; however, our recent study suggests the isomerization becomes reversible again in the ensuing S 3 Y Z • state after H + release 93 .O-O bond formation in the S 4 state (peroxo in S 4 ') for either open-11,94-97 or closed-cubane 27,31,98-100 structure (with O5 and Ox as substrates) has been theoretically supported, and the Im0 ÀO 2 structure formed after O 2 release would be basically the same; thus the missing link would apply in either case.Note that the protonation states for some certain ligands are still debatable, such as W2 = OH − /H 2 O [34][35][36][37][38]101,102 , O4/O5 = O 2− /OH − (for S 0 ) 34,37,38,64,103   spherical cavity surrounded by the conductor-like screening model (COSMO) of the polarizable continuum with a dielectric constant ε = 6.0, implicitly mimicking the protein matrix. The Busi-Parrinello Langevin dynamics (T damp = 1000 fs) for the thermostat and timereversible integrator with dissipation for self-consistent field (SCF) were adopted for the NVT canonical ensemble (particle Number, Volume, and Temperature) simulations at 298.15 K.The simulations were pursued for at least 30 ps with velocity-Verlet integration using a 1.0 fs time step. All the BO-AIMD simulations were executed by the commercially available GPU-accelerated package TeraChem 84 (version 1.94) on supercomputers equipped with NVIDIA Tesla V100 cards.

MEP calculations
Truncated DFT models were constructed from the last snapshots of the BO-AIMD simulations for the subsequent MEP searches, without missing out possible influence from surrounding necessary groups.This approximation represents a sensible compromise based on consideration on the feasibility and efficiency for massive Hessian calculations.For W1 (the new W2) deprotonation to Asp61, the model includes the inorganic Mn 4 CaO 6 cluster, 10 amino acid residues Tyr161, Glu189, His190, His332, Glu333, His337, Asp342, Ala344, Glu354, and Arg357, Asp61 and 12 crystal water molecules W1-W8, HOH596, HOH545, HOH525, and HOH526, resulting in total 223 atoms and a net total charge of +1 (Supplementary Fig. 8).According to the test by Retegan et al. 85 , a QM model with such a size is adequate to yield accurate reaction energetics and spectroscopic properties for the OEC system.For the subsequent W2(H 2 O) dissociation from Mn4, the starting geometry for MEP was based on the product state of Asp61 protonation by exclusion of the protonated Asp61 from the model (the reason is described in Supplementary Note 14, Supplementary Figs. 9 and 10).For the later μ 3 -W3(OH − ) ligand transfer for the closed-open-cubane transition, the starting geometry for MEP was based on the product state of W2 dissociation by further removal of W2 (assumed to release to the bulk).Test computations based on the original crystal structure 6DHP including Asp61 in the model justify the validity of the conclusion (Supplementary Figs. 12 and 13).The proton-released state is termed "pre-S 0 ", as shown in Fig. 3, specifically pre-S 0 (W2-bound) and pre-S 0 (W2-unbound) according to the bound state of W2 in Fig. 4. In the next step, the S 0 B (closed-cubane) state is formed after W2 is released to the bulk, and then the S 0 A (open-cubane) state is formed by μ 3 -W3(OH − ) ligand transfer.Geometry optimizations with backbone constraints (α-carbons fixed along the peptide chains) were performed by the unrestricted hybrid functional B3LYP* 86 (15% exact exchange, dispersion parameters taken from B3LYP, Supplementary Note 11), supplemented by various different dispersion-parameterized DFT functionals for comparison where the obtained relative free energies are small (within a few kcal mol −1 ).The LanL2DZ and 6-31 G* basis sets were used for Mn/Ca, and the rest H, C, N, O, and Cl atoms, respectively.Analytic frequency calculations on the optimized structures at the same level of theory verified all local minima; zero-point energies (ZPE) and thermal effects (298.15K, 1 atm) were also extracted for thermal correction to Gibbs free energies.Transition states were confirmed through eigenvectors with adequate and expected negative eigenvalues, single imaginary frequency vibrations, and intrinsic reaction coordinate (IRC) analyses (see Supplementary Note 15 for more details), ensuring that the relationship with reactants and products are logical and correct.TSs were located by the Berny algorithm and synchronous transit-guided quasi-Newton (STQN) method.Finally, more accurate single point energies were computed with SDD (for Mn/Ca) and cc-pvtz (-f) (for H, C, N, O) basis sets under SMD continuum solvation model (solvent-accessible surface, ε = 6.0).DFT-D3 with Becke-Johnson (GD3BJ) damping dispersion correction was applied in both geometry optimizations and single-point energy calculations.All these computations were carried out by Gaussian 16 (version C. 01) 87 .Formal oxidation states for Mn were verified by Mulliken spin populations and localized orbital bonding analysis (LOBA) 88 using Multiwfn 89 (version 3.8).

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Fig. 1 |
Fig. 1 | Overview of the final transition of Kok's cycle.a Exhibitions of the XFEL crystal structures of the OEC cluster in the S 3 and S 0 states 4,8 ; the stage included in the S 4 → S 0 transition after O 2 release (highlighted at the top of the cycle) represents the focus in this study.b Molecular diagram for the S 3 → S 0 transition via a putative S 4 state where O-O bond formation occurs through O5-Ox• radical coupling (shown in the dashed box); the partial ligands in the S 0 state are not labeled because of possible ligand rearrangement during water insertion; the oxidation and protonation states are shown based on the most widely accepted forms12,21,30,34,64,90 .

Fig. 2 |
Fig. 2 | Graphic presentations of water insertion after O 2 release.a The S 4 state based on the XFEL structure of the S 3 state 6 upon removal of one proton and one electron on Ox; the yellow arrow represents the formation and release of O 2 through O5-Ox coupling.b The Im0 ÀO 2 state derived from the S 4 state by removal of O5 and Ox for O 2 release.c The Im1 state formed by W3 insertion to Mn1 and W5 binding to Ca. d The Im2 state formed by W3 deprotonation to W2 and binding to ÀO 2 → Im1 → Im2 is recognized as a very rapid conversion during the S 4 → S 0 transition.The pre-S 0 state formed by the deprotonation of Im2 is still structurally different from the final S 0 A (open-cubane) state, because of its closed-cubane conformation and one additional water ligand.However, further progression from pre-S 0 to S 0 A , via stepwise W2(H 2 O) dissociation and μ 3 -W3(OH − ) ligand transfer, are proven feasible (Fig. 4 and Supplementary Tables

a
Fig. 3 | Structural evolution of the Mn 4 Ca cluster.a Variations of the key interatomic distances with time evolution along the simulation trajectories for the octet/αααβ and doublet/αβαβ spin states.b Chemical structural formulas for the three intermediates; they all correspond to their respective local minima (stationary points) on the potential energy surfaces.

Fig. 4 |
Fig. 4 | Water (W2) dissociation in the pre-S 0 state.Relative Gibbs free energy profiles for the octet/αααβ and doublet/αβαβ spin states are shown in blue and green, respectively; key interatomic distances are displayed in Ångström exemplified by the doublet/αβαβ spin state; these also apply for Fig. 5.

Fig. 5 |
Fig. 5 | Structural isomerism by μ 3 -OH − (W3) ligand transfer in the S 0 state.Metal oxidation states are labeled in Roman numerals and the J-T elongation axes are marked in orange.Spin up/down is marked in yellow/burgundy at the atomic centers involved (contour value for spin density set to 0.1).Note that all these parameters show slight model/functional dependence (see below).

Fig. 6 |
Fig. 6 | Suggested missing link and isomerism in Kok's cycle.The missing link during the S 4 → S 0 transition embedded in the full catalytic cycle is highlighted in red, and the isomerism in green.O5 in S 1 B is at a non-bonding distance with respect to both Mn1 and Mn4 67 .in the S 3 state denotes the unidirectional isomerization because the open-cubane structure is largely stabilized12,91,92 ; however, our recent study suggests the isomerization becomes reversible again in the ensuing S 3 Y Z • state after H + release93 .O-O bond formation in the S 4 state (peroxo in , and Ox = OH − /O 2− /O •− (for S 3 /S 3 Y Z •)4,5,12,38,97,104 , etc.