Water Network Dynamics Next to the Oxygen-Evolving Complex of Photosystem II

Abstract

The influence of the environment on the functionality of the oxygen-evolving complex (OEC) of photosystem II has long been a subject of great interest. In particular, various water channels, which could serve as pathways for substrate water diffusion, or proton translocation, are thought to be critical to catalytic performance of the OEC. Here, we address the dynamical nature of hydrogen bonding along the water channels by performing molecular dynamics (MD) simulations of the OEC and its surrounding protein environment in the S₁ and S₂ states. Through the eigenvector centrality (EC) analysis, we are able to determine the characteristics of the water network and assign potential functions to the major channels, namely that the narrow and broad channels are likely candidates for proton/water transport, while the large channel may serve as a path for larger ions such as chloride and manganese thought to be essential during PSII assembly.

1. Introduction

Photosystem II (PSII) is a 650 kDa complex composed of 20 proteins embedded in the thylakoid membrane of green plant chloroplast and internal membranes of cyanobacteria. PSII sustains aerobic life on Earth by harvesting solar light, initiating the process of photosynthesis, and producing oxygen from water. The water oxidation reaction is catalyzed by the oxygen-evolving complex (OEC), which is a CaMn₄O₅ cuboidal cluster embedded in the D1 protein subunit. Oxygen is produced through the catalytic Kok cycle. In each turn of the cycle, the OEC evolves through five storage states (S-states) of oxidizing equivalents, with S₀ and S₄ being the most reduced and oxidized states, respectively [1,2]. The stable dark-adapted S₁ state (structure and oxidation states shown in Figure 1) is advanced to the S₂ state by a single flash of light and generates a mixture of two redox isomers. Isomer A is an “open” cubane structure where Mn1 and Mn4 have formal oxidation states III and IV, respectively [3]. Isomer B is “closed” cuboid where the oxidation states of Mn1 and Mn4 are switched to IV and III, respectively. S_2A exhibits a characteristic multiline EPR signal at g = 2.0 [4], while S_2B has a broad EPR band at g = 4.1 [5,6]. A second flash moves the OEC into the S₃ state, where all four Mn centers have formal 4+ charges and an additional water is bound between Mn1 and Mn4, setting the complex up for oxygen evolution [7]. A final third flash brings the complex into the S₄ state, however this state is highly unstable and almost immediately decays to the S₀ state, releasing oxygen in the process. As the S₄ state decays faster than it is formed, its structure has not been experimentally characterized, although several theoretical models have been proposed. From the S₀ state, the OEC can return to S₁, either through a fourth light flash or gradual decay in darkness, although some fraction of OEC centers will remain in the S₀ state, even after extensive dark adaptation [8].

The varying shape and oxidation states of the OEC make it a challenging system to model. Typically, quantum mechanics (QM) and two-layer QM/MM models are used to study the OEC and its surrounding protein [9,10]. These methods are computationally expensive and limit the size of the system that can be studied. Classical methods, such as molecular mechanics, provide practical approaches for computational modeling, although with significant shortcomings. Force fields require significant parameterization, especially for nonstandard residues and metal-containing systems, such as the OEC. Classical molecular dynamics simulations have been performed on PSII and get around the parameter problem by making the OEC and other metal cofactors static [11,12]. Valence charges are assigned to a non-bonded model of the OEC and the complex and its ligating residues are locked into their crystal structure or QM-calculated coordinates. While the majority of PSII is indeed dynamic in such simulations, the OEC itself is not. One of our ongoing goals has been to derive classical parameters for the OEC and its ligands. Several methods are available to generate parameters for metal complexes, however none have been tested on a system containing as many metal centers and high oxidation states as the OEC. In an effort to make the OEC more dynamic during these simulations, we use a flexible model with ESP-derived partial charges and strong force constants for all metal containing bonds and angles, 500 kcal/mol/Ų and 100 kcal/mol/rad², respectively.

Dynamic simulations of the OEC are of major interest because they allow for the study of water passages and the effects of changing S-states on these passages. Previous studies have examined the pathways water takes to and from the OEC [13,14,15]. Molecular dynamics simulations have indicated that the two sites for water binding to the OEC has completely separate sources, with waters bound to the calcium (W₃ and W₄) originating from the large channel and the water bound to Mn4 (W₁ and W₂) being supplied by the narrow or broad channels [13]. While the broad and large channels connect near D1-E189, there is little, if any, crossover between the two [14]. While the large channel has been widely accepted as the source for W₃ and W₄, it is unclear whether W₁ and W₂ flow from the narrow or broad channel, or even the large channel. Though the broad channel is significantly wider than the narrow channel, as indicated by their names, the narrow actually has fewer and lower permeation barriers than the broad channel [14]. Part of the interest in these channels arises from the implications their purposes have on the mechanism of forming the S-states of the OEC, which requires two new waters to bind to the OEC for every completed turn of the Kok cycle. One of these waters binds during the S₂-to-S₃ transition, as described above [7]. W_X, sharing a hydrogen bond with O4, is the primary candidate for a newly-bound water during this transition, although the mechanism and where it binds is unclear. One mechanism proposes that W_X originates from the narrow channel and binds to Mn4 in the S_2B isomer as part of a “carousel” mechanism [16]. A second possible mechanism proposes that W_X flows up from the large channel and replace W₃ on the calcium after that water binds to Mn1 [17]. Studies of the binding of methanol [18] and ammonia [19,20] near or directly to the OEC support the carousel mechanism. Ammonia binds to Mn4 as an additional ligand in the S₂ state, and neither ammonia nor methanol shows any interaction with Ca²⁺. These are strong indicators as to where water binds in the S₃, although it does not necessarily show which channel substrate water molecules come from before entering the W_X pocket. In addition to supplying waters, channels are also needed for proton shuttling, oxygen evacuation, and OEC reassembly. In the present work, we utilize a flexible model of OEC and network theory to characterize the water dynamics in the surroundings of the complex in the S₁ and S₂ states. Our analysis clearly identifies three channels of highly diffusive and interconnected waters in direct contact with the OEC complex. Each of these channels display different water mobility and interconnectivity, as well as different sensitivity to the advance from S₁ to S₂ state. It is shown that restricted mobility might lead to a stronger hydrogen bonding network. In light of these results we postulate different roles played by the identified channels in the OEC activity.

2. Results

2.1. Water Network Dynamics

The mobility of the interstitial water molecules in PSII provides insights on potential diffusion pathways for H₂O, H⁺, or O₂. Therefore, we computed the Root Mean Square Fluctuation (RMSF) per water molecule to quantify the mobility of waters in the system. Figure 2 (left column) shows that there are two main channels of highly diffusive waters that connect the Mn₄CaO₅ complex with the lumen. The narrow channel shows a straight pathway of hydrogen-bonding water molecules with moderate mobility. The second large channel shows a branching point and much higher mobility towards the protein exterior. A third channel bridges the narrow and large channels with intermediate water mobility values. This bridge channel matches the entry of the broad channel described in Reference [18], which splits into two more channels (the broad and a proposed proton-exit channel) between D1–D61 and D1–D65 (not shown here). Interestingly, the structure and mobility of these water channels show sensitivity with respect to the S₁ and S₂ states studied here, suggesting a potential regulatory role of the water dynamics on OEC activity.

In order to further analyze the influence of the surrounding waters on the Mn₄CaO₅ complex, we computed the Eigenvector Centrality (EC) in terms of the Mutual Information (MI) for all the waters in the system [21]. The measure of centrality gives an idea of how much each atom influences the dynamics of neighboring atoms, i.e., how connected each water molecule is with its neighbors. It is expected that higher connectivities can be related with high proton mobility, and therefore, an increased effectiveness in proton transport between the OEC and the exterior of the protein.

This analysis shows that the same three channels that exhibit the highest mobility also show the highest EC within the water network (Figure 2, right column). Again, it is shown that the EC distribution is very sensitive to both the state transition from S₁ to S₂ and the isomerization S_2A to S_2B. On the other hand, it is interesting to note that the waters in these two channels manage to solvate the Mn₄CaO₅ complex, and therefore, long chains of water with great centrality reaching this complex may indicate water-assisted long-range modulation of the reaction taking place in the OEC.

2.2. Diffusion vs. Connectivity

To provide a comprehensive analysis on how the properties of the previously identified water channels change in the different states of the OEC analyzed here, we computed the number of hydrogen bonds per water molecule in each channel (Figure 3). The observed number of hydrogen bonds per water is considerably smaller than the corresponding value in bulk conditions (~4). This is reasonable, considering that the water coordination is restricted inside the protein cavities in comparison with the bulk. Furthermore, several of the hydrogen bond interactions between water molecules are replaced by specific interactions with the amino acids near the OEC. On the other hand, it is clear that the number of hydrogen bonds between water molecules is larger in the narrow channel than in the large or broad cavities, except in the S₁ state. This can be explained by the smaller mobility of the water molecules in the narrow channel (Figure 2). The more rigid cavity constrains the waters in positions where the hydrogen bonding is favored in comparison with the large or broad channels. Therefore, while the narrow channel could be more appropriate for proton transport, the large channel would be the most favorable for the diffusion of larger ions through the water network. With respect to changing S-states, the large channel is largely unaffected.

The narrow channel is the most sensitive, with the distribution of hydrogen bonds narrowing and increasing from the S₁ to S_2A states. In S_2B, the mean similarly shifts upward and distribution noticeably broadens compared to the S_2A state. Considering the placement of the narrow channel close to Mn4, it is reasonable that the cavity would be so affected by the state transitions since Mn4 is oxidized or changes its coordination in the S_2A and S_2B states, respectively. Additionally, the bridge channel has its most dramatic change in the S_2A state, with its distribution of hydrogen bonds slightly narrowing and shifting downward, although this is nowhere near as dramatic as in the narrow channel. Again, this modification can be attributed to the changes to the OEC in this state, with the oxidation of Mn4, which the bridge channel runs alongside. It is important to note that the great mobility and connectivity between the waters in the intersection of the narrow and bridge channels (Figure 2) indicates that these channels are not isolated from one another and water or protons can flow between them.

Figure 4 (left panel) shows the total root mean square deviation (RMSD) of the atomic positions associated to the water molecules in the three channels. Similar to the RMSF per water, the RMSD accounts for the overall standard deviation of atomic displacements, reflecting the total mobility of the waters in the system. It is clear that water diffusion is more favorable in the large channel. This supports the hypothesis that the large channel is more suitable for the migration of larger ions. On the other hand, the normalized first eigenvalue of the correlation matrix (Figure 4, right panel) is a measure of the overall connectivity between waters in each channel. This value indicates that the bridge channel has the highest connectivity, with the exception of the S_2A state where the narrow channel is more connected. This observation supports the idea that the broad channel is part of a proton exit pathway.

3. Discussion

We have analyzed the dynamic networks established by interstitial waters adjacent to the OEC of photosystem II by using molecular dynamics simulations. We have identified three channels of highly diffusive and interconnected waters. One of these channels, named the narrow channel, is located in a rigid cavity where water movement is restricted and hydrogen bonding between waters is enhanced. Conversely, in the large channel, waters are much more mobile and hydrogen bonding between waters is more fluctuating. The bridge channel between the two, whose location matches well with the entry of the broad channel, has intermediate character between the other two channels. Therefore, while the structural and dynamical properties of the narrow channel, and the broad channel to a lesser extent, should be ideal to enhance Grotthuss-like proton diffusion, the large channel might be more appropriate for diffusion of larger ions.

Previously, both the narrow and large channels have been identified as possible candidates for supplying substrate waters to the OEC, specifically the W_X molecule near O4 [16,17]. Generally, the narrow channel is favored by experiments [22], while the large channel has computational support [17]. Meanwhile, the broad channel has been proposed as a possible exit channel for proton release to the lumen during the Kok cycle [23]. Based on our analysis, the large channel should be ideal for transport of both water and larger ions. This could be key for the periodic reassembly of the OEC and for the delivery of chloride ions. The distant chloride ion, Cl-2, near D1–N184, is located very close to this channel, further indicating that this channel is very likely to be its source. However, the other chloride, Cl-1, near D2–K317, is closer to the broad channel and has no clear path to the large channel. Therefore, although the large channel is the best for delivering ions, it is more likely that the broad channel is the preferred pathway for Cl-1. The channel we have identified as the bridge is part of the broad channel and should be capable of facilitating such a transport.

The narrow channel has been shown by our analysis to be the most rigid of the three near-OEC channels. Though part of the broad channel has been put forward as the proton exit channel, our results suggest that the well-ordered narrow channel is the best pathway for protons. The increased rigidity of the narrow channel in the S₂ state could explain the reduced occupancy of the water molecule next to W_X in the narrow channel in the 1F (S₂) structure observed by Jan Kern et al. [24] and Jimin Wang et al. [25]. If the water molecules are stabilized and unable to move in the S₂, water may temporarily be unable to fill this small pocket in the narrow channel until its reappearance in the 3F (S₀) structure. This missing water in the S₂ state would create a gap in the chain of hydrogen-bonded water molecules in the narrow channel. If this channel is the proton exit path, this gap would be a barrier for exiting protons which would explain why a proton is not released in the S₁→S₂ transition. However, it was also shown in one structure (3ARC/5V2C) [26,27] that the narrow channel can be blocked by glycerol without apparent inhibition of oxygen evolution. Assuming the presence of this molecule is not an artifact of the experiment, then this channel cannot participate directly in water or proton transport. Instead, the narrow channel could simply be stabilizing the hydrophilic pocket near O4 that W_X occupies before binding to the OEC in the S₂→S₃ transition.

Experimentally, it has been shown that depletion of chloride from PSII prevents the OEC from progressing past the S₂ state, indicating that the proton exit pathway is blocked [28,29]. Mutation of D2–K317 to alanine eliminates the chloride dependency, showing that it is the formation of a salt bridge between D1–D61 and K317 that prevents proton release [30]. D61 is located at the junction between the narrow and broad channels, so the aspartate-lysine salt bridge could block either or both channels for a proton released from substrate water on Mn4. As far as the delivery of water molecules to the OEC, our analysis indicates that the broad and large channels are the best candidates, with the narrow channel displaying low water mobility. A recent computational study looked into how water rehydrates a dry PSII and found that the water molecules occupying the W_X position near O4 originated from the broad channel, not the large [13]. If W_X does enter through the broad channel, then the rigidity of the narrow channel would not rule out the carousel mechanism as the method by which new water molecules bind to the OEC. To summarize, our results indicate that the large channel is the most likely supplier of ions to the OEC cavity, with the exception of Cl-1, which is likely supplied by the broad channel. The narrow channel should be suitable for proton removal, with the broad channel being another possibility. The water mobility of the large and broad channels makes them almost definite suppliers of substrate waters to the OEC. Future analysis will be done using a larger model, such as a complete PSII, in order to determine how far the effects of these channels reach.

4. Materials and Methods

4.1. MD Simulations

Simulations were run using a sphere of residues and cofactors within 25 Å of the OEC. Chain cuts were capped with acetyl (ACE) and N-methylamide (NME) residues and these caps were frozen during simulations. Additionally, heavy atoms within water molecules and cofactors more than 15 Å away from the OEC were also frozen to their crystallographic coordinates. All hydrogen and heavy atoms within 15 Å were free to move. The calcium and manganese centers of the OEC were fit with the ESP charges using the Merz-Singh-Kollman (MK) method [31] in Gaussian09 [32]. ESP charges were fit to a sidechain and backbone model of the OEC and all its directly-bound ligands. This model was calculated with DFT/B3LYP [33,34] using the LanL2DZ [35] basis set and pseudo-potential for calcium and manganese atoms and 6-31(d) [36,37] for all other atoms. Ligand and µ-oxo bridge partial charges were calculated using the RESP [38] method with metal centers present but restrained to their ESP values. A RESP fit of the entire OEC was attempted, however, this resulted in near-zero metal charges and negatively-charged hydrogens. Holding the metal centers to ESP charges solved both of these issues. All metal-containing bonds and angles were held to their QM/MM optimized structure [7,8] with 500 kcal/mol/Ų and 100 kcal/mol/rad², respectively. These force constants were used in lieu of those calculated by the Seminario method [39], as the values estimated by the latter method were too weak (<100 kcal/mol/Ų) to hold the system together. The stronger force constants retained the complex structure without completely restraining it in space. The TIP3P force field [40] was used for all water molecules in the system. Bonds and partial charges were calculated for the S₁ (III IV IV III), S_2A (III IV IV IV), and S_2B (IV IV IV III) states. Parameters for metal-containing cofactors were created using the MCPB.py [41,42] script and other cofactors using Antechamber, both available through the AmberTools package [43]. As per the MCPB protocol, metal centers and the atoms bound to them were assigned arbitrary atom types (e.g., M2 or Y6) for which parameters were created. The remaining ligand atoms were assigned their usual Amber atom type, as determined by their amino acid and atom names (e.g., C_γ in an aspartate has Amber atom type CO). Simulations were run in Sander for 100 ns under the Amber 14 [44] force field at 300 K after heating from 0 K for 20 ps, after which the variance of the temperature was less than 3 K, indicating achievement of a stable simulation temperature. Example files containing the charge and bonding parameters for the S₁ state, as well as a sample simulation, can be found in the Supplementary Materials.

4.2. Eigenvector Centrality Distribution for Water Networks

To shed light on the collective organization of interstitial waters in PSII, we utilized the eigenvector centrality (EC) methodology recently developed in our group [45]. This technique has successfully identified allosteric pathways in protein residue networks. The EC of a node, $c_i$, is defined as the sum of the centralities of all nodes that are connected to it by an edge, $A_{ij}$,

$$c_i=\frac{1}{\lambda }{{\displaystyle \sum }}_{j=1}^n{A}_{ij}{c}_j$$

(1)

where the edges $A_{ij}$ are elements of the so-called adjacency matrix $A$, and $\lambda$ is the eigenvalue associated to the eigenvector composed by $c_i$ elements. In particular, our approach relies on assigning the functional dynamics to the major collective mode of the system (i.e., the first eigenvector of $A$), and hence in this context, $\lambda$ is the highest eigenvalue of $A$. There are several ways of defining an adjacency matrix, and for insight into the interconnectivity of a water network, the focus was put on the correlated motions of the system. Therefore, we utilized an adjacency matrix weighted by the dynamical correlations between protein residues, the elements of which can be quantified by the generalized correlation coefficient:

$$G_{MI}\left[r_i,r_j\right]={\left[1-\left({\mathrm{e}}^{-\frac{2}{3}I\left[r_i,r_j\right]}\right)\right]}^{-\frac{1}{2}}$$

(2)

where $I\left[r_i,r_j\right]$ are expressed in terms of the mutual information [21],

$$I\left[r_i,r_j\right]=S\left[r_i\right]+S\left[r_j\right]-S\left[r_i,r_j\right]$$

(3)

being $S\left[r_i\right]$, $S\left[r_j\right]$, and $S\left[r_i,r_j\right]$ the marginal and joint Shannon entropy [46], for atomic displacements $\left(r_i,r_j\right)$, obtained as ensemble averages over the MD simulation. The EC is a measure of how well connected a node is to other well-connected nodes in the network. In addition, the EC serves as a measure of the connectivity against a fixed scale when normalized, and hence, at variance with other centrality metrics, it can be used to reliably compare different networks.

5. Conclusions

Using a flexible model of the OEC, we conducted an EC analysis on the water channels surrounding the OEC in the S₁ and S₂ states. Due to its low mobility and high connectivity, we concluded that the narrow channel near O4, which had previously been proposed as a substrate water channel, is more suited to transporting protons rather than molecules such as water. The narrow channel is also the only channel to show significant change between S-states, becoming more rigid in the S₂ state, possibly in preparation to shuttle away a proton in the S₂→S₃ transition. The large and broad channels are mobile enough to easily transport molecules and ions. Because of this and the lack of accessibility by the large channel to the Mn4 binding site, we propose the broad channel as a supplier of substrate water molecules.