Dynamics of Long-Distance Hydrogen-Bond Networks in Photosystem II

of local, lower-occupancy H-bonds that were not present in the crystal structure. All-atom molecular dynamics simulations allow us to derive a detaile...
0 downloads 9 Views 3MB Size
Subscriber access provided by UNIVERSITY OF TOLEDO LIBRARIES

B: Biophysical Chemistry and Biomolecules

Dynamics of Long-Distance Hydrogen-Bond Networks in Photosystem II Federico Guerra, Malte Siemers, Christopher Mielack, and Ana-Nicoleta Bondar J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b00649 • Publication Date (Web): 28 Mar 2018 Downloaded from http://pubs.acs.org on March 28, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Dynamics of Long-Distance Hydrogen-Bond Networks in Photosystem II

Federico Guerra, Malte Siemers, Christopher Mielack, and Ana-Nicoleta Bondar*

1

Freie Universität Berlin, Department of Physics, Theoretical Molecular Biophysics, Arnimallee 14, D-14195 Berlin, Germany

*Corresponding Author E-mail: [email protected], Phone number: +49-30-838-53583

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 55

Abstract Photosystem II uses the energy of absorbed light to split water molecules, generating molecular oxygen, electrons and protons. The four protons generated during each reaction cycle are released to the lumen via mechanisms that are poorly understood. Given the complexity of photosystem II, which consists of multiple protein subunits and cofactor molecules and hosts numerous waters, a fundamental issue is finding transient networks of hydrogen bonds that bridge potential proton donor and acceptor groups. Here, we address this issue by performing all-atom molecular dynamics simulations of wild type and mutant photosystem II monomers, which we analyze using a new protocol designed to facilitate efficient analysis of hydrogen-bond networks. Our computations reveal that local

protein/water

hydrogen-bond

networks

can

assemble

transiently

in

photosystem II, such that the reaction center connects to the lumen. The dynamics of the hydrogen-bond networks couple to the protonation state of specific carboxylate groups, and are altered in a mutant with defective proton transfer. Simulations on photosystem II without its extrinsic PsbO subunit provide a molecular interpretation of the elusive functional role of this subunit.

ACS Paragon Plus Environment

2

Page 3 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Introduction Photosystem II is a nanometer-size biomolecular complex found in plants and other oxygenic photosynthetic organisms that use the energy of sun light to split waters, generating molecular oxygen that sustains biological life. The chemical reaction of water splitting generates protons and electrons, which need to be removed from the vicinity of the reaction site (denoted as the oxygen-evolving complex, OEC, in Fig. 1A, Fig. 1B). The mechanism that ensures transfer of the protons to the lumen is poorly understood, though it is thought to involve hydrogen (H)-bond networks of protein groups and water molecules. Identifying the H-bond networks that could serve as proton-transfer pathways is a key step towards deciphering the atomic mechanism by which photosystem II works. To identify such networks and characterize their dynamics, here we have carried out all-atom molecular dynamics simulations of photosystem II and implemented an algorithm for efficient analyses of dynamic H-bond networks. Crystal structures of photosystem II

1-6

have provided invaluable information on

the architecture of photosystem II, including H-bond networks of interest as potential proton transfer pathways

2, 7-12

. The analysis of hydrophilic channels found in crystal

structures of photosystem II allowed potential exit channels to be identified. One such channel includes carboxylate groups of subunits D1 (D61, E65 and E333) and D2 (E312), and other protein groups that can hydrogen bond

8, 13

. Water molecules,

which could play important roles in proton transfer reactions of photosystem II

11, 13

,

were identified in a hydrophilic channel that could be involved in proton transfer 7. Numerous internal water molecules could be identified in the high-resolution crystal structure from ref. 2, including water molecules that participate in a long-distance

ACS Paragon Plus Environment

3

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 55

hydrogen-bond network that connects the OEC to the lumen via groups of subunits D1, CP43 and PsbV 2. A particularly important network that includes protein groups and water molecules extends from D61 of subunit D1 to the lumen via the carboxylate pair D1-E65 and D2-E312

10

, and D224 of the PsbO subunit (Fig. 1C, Fig. 1D, Fig. 1E). The potential

involvement of D1-E65 has been recognized based on earlier structural analyses D1-E65 is also part of a potential proton channel found in ref.

7

14

;

based on channel

calculations. Crystal structures indicate that the carboxylate groups of D1-E65 and D2-E312 are within H-bond distance (2.5Å in ref. 2, 3.3Å in ref. 1), and pKa computations suggest that D2-E312 might be protonated

9, 15

. At the lumen end of

the extended H-bond network (Fig. 1E), D224 is part of a larger cluster of PsbO carboxylate groups thought to function as a proton antenna 10, 12, 16 - that is, a cluster of carboxylate groups located within close distance, which can transiently bind a proton on a protein’s surface

17

. Thus, the molecular picture arising from the crystal

structures is that an extended H-bonded network that includes carboxylate sidechains and water molecules (Fig. 1D), might help release protons from the vicinity of the OEC to the lumen. Based on the crystal structure alone, however, it remains unclear whether this extended H-bonded network indeed serves as proton transfer pathway: Complex Hbond networks can be rather dynamic, and the dynamics of H-bond clusters in protein environments cannot be anticipated based on static crystal structures Indeed, recent Fourier-Transform Infrared (FTIR) spectroscopy

21-22

18-20

.

data on

photosystem II were interpreted to suggest that near the OEC there is an extended network of H-bonds that includes carboxylate groups whose pKa values or environment change during specific steps of the reaction cycle of photosystem II

ACS Paragon Plus Environment

22

.

4

Page 5 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The H-bond network proposed based on FTIR data includes protein groups identified based on inspection of the crystal structure – e.g., D1 groups E65 and E329, and D2-E312, and water molecules 22. The dynamics associated with altered pKa values or altered environment of specific carboxylates from the extended H-bonded network 22 remains unclear at the molecular level of detail. Knowledge of these dynamics is essential for the transfer mechanism, because motions of potential proton carriers allow proton-transfer paths to be established, and can lower the energy barrier for proton transfer 23-24. A related important question is how water molecules, whose hydrogen atoms are not solved in the crystal structure 2, contribute to dynamic H-bond networks and how they interact with carboxylate groups that might carry a proton. A plausible scenario here is that dynamics of protein groups and waters associate with transient breaking and forming of H-bonds within the extended network (Fig. 1E), and/or with the formation of local, lower-occupancy H-bonds that were not present in the crystal structure. All-atom molecular dynamics simulations allow us to derive a detailed view of the motions of a protein and bound waters at room temperature, in fluid environments. Importantly, such simulations allow us to probe the dynamics of H-bonded protein/water networks

18, 25

, and to explore the response of the H-bond networks to

changes in the protonation state of specific protein groups 19, 26. Here, we relied on all-atom molecular dynamics simulations to explore the dynamics of photosystem II in a hydrated lipid membrane environment (Fig. 1C), and the response of protein and bound waters to changes in the protonation state of specific carboxylate groups or to mutation. The simulation system includes the 19 protein subunits of the photosystem II monomer from Thermosynechococcus vulcanus 2, with more than 70 cofactor molecules and special lipids, embedded in a

ACS Paragon Plus Environment

5

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 55

lipid membrane patch composed of phosphatidylcholine (POPC) lipids. Since the identity of protonated carboxylate groups that could be involved in proton transfer remains unclear, we studied the dynamics of photosystem II for three different protonation states chosen based on considerations of carboxylate sidechain interactions in the starting crystal structure. We augmented the simulations on protonation-dependent dynamics of photosystem II with independent simulations where we probed the dynamics of the D1-E65A and ∆PsbO mutants. We chose the D1-E65A mutant by considering interactions observed in simulations, and experiments that indicate altered reaction cycles of this mutant photosystem II complex 22. The ∆PsbO mutant allows us to probe the functional role of PsbO. A specific challenge when studying H-bond dynamics of photosystem II is the size of the H-bond data sets that need to be analyzed: The lumen-exposed domain of photosystem II contains >500 sidechain–sidechain, sidechain–backbone Hbonds. To identify and characterize efficiently the dynamics of H-bond networks in photosystem II, we designed and implemented a data analysis tool inspired by graph theory. The all-atom molecular dynamics simulations allow us to identify two extended Hbond networks that connect the region of the OEC to the bulk, with exit groups contributed by the PsbO and CP43 subunits, respectively. The D1-E65A mutation impacts significantly the dynamics of the H-bond networks, with extensive changes in the dynamics of specific H-bonds, including at the PsbO interface. In simulations where PsbO is absent, we observe rearrangements of local H bonds and rapid loss of an internal chloride ion thought to be important for function. Overall, our simulations underlie the importance of protein and water dynamics for the transient formation of H-bonded networks that could participate in proton transfers. Indeed, H-

ACS Paragon Plus Environment

6

Page 7 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

bond networks that have been discussed before based on static crystal structures appear to have complex dynamics, whereby some protein groups bridge persistently via H-bonding waters, whereas for other groups water bridging can be rather dynamic.

Methods Protein structure. We used as starting coordinates monomer 1 of the crystal structure of T. vulcanus photosystem II from ref.

2

(PDB ID: 3WU2 monomer 1, i.e.,

all protein subunits whose chains are labeled with capital letters, and the cofactors, special lipids, waters and ions associated with these protein chains). This monomer contains coordinates for 19 protein subunits, all of which were included in our computations. In the crystal structure, chains A (subunit D1), chain K (PsbK), and chain M (PsbM) contain mutations R279P, F33L/F39W and F8L, respectively. To facilitate direct comparison of simulation results with the starting crystal structure, these mutations were preserved. For simplicity, in what follows we refer to the photosystem II complex that retains the crystal structure mutations as ‘wild-type photosystem II’; we specify explicitly a mutant photosystem II only when we introduced an additional mutation. The total number of amino acid residues included in our photosystem II model is 2791. The crystal structure of photosystem II includes coordinates for 54 cofactor molecules, and 20 special lipid molecules. We included in our model all cofactor and special lipid molecules that were solved in the crystal structure for monomer 1, except for those molecules located at the external membrane interface; for simplicity, these surface molecules were removed. Several cofactor and special lipid molecules are incomplete or missing from the crystal structure of ref. 2. To construct

ACS Paragon Plus Environment

7

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

coordinates for these molecules, we used structure information from ref.

Page 8 of 55

27

. Our

resulting model of photosystem II includes 34 chlorophyll-a molecules, 2 pheophytin, 2 heme and 2 plastoquinone molecules, one non-heme iron/bicarbonate complex, and 8 β-carotene molecules. In addition to the protein, cofactor and special lipid molecules, we included all ions and water molecules belonging to monomer 1 from ref. 2. Hydrogen atoms coordinates were built using the CHARMM Hbuild command. To prepare the lipid membrane simulations, we first used the OPM server

28

to

orient the photosystem II monomer relative to the membrane normal. In the next step, we used CHARMM-GUI

29

to embed photosystem II in a 200x200 Å hydrated

membrane patch composed of 910 phosphatidylcholine (POPC) lipid molecules, ~135,000 water molecules, and potassium ions for charge neutrality. The simulation system of photosystem II in hydrated lipid membrane contains a total of ~580,000 atoms. Our choice of POPC for the bulk lipid molecules, which are not present in thylakoid

membranes,

is

reasonable,

because

compatible with the functioning of photosystem II

phosphatidylcholine

30

appear

, and the thickness of POPC

membranes is within 1Å of the cyanobacterial thylakoid membrane 31. The photosystem II mutants E65A and ∆PsbO were prepared with the CHARMM software 32 using as starting coordinates the simulation system described above.

Force field description. We used the CHARMM amino acid residues

33

32

with the CMAP correction

force field parameters for protein 34

, CHARMM lipid parameters

35

and ions 36, and the TIP3P water model 37.

ACS Paragon Plus Environment

8

Page 9 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

For chlorophyll-a, pheophytin-a and plastoquinone-9, we employed parameters previously derived by us

38

. Of the 34 chlorophyll molecules included in our

simulation system, 28 have their magnesium ion coordinated by histidine or asparagine sidechains, and 6 chlorophyll molecules that are coordinated by a water molecule. To maintain the structural integrity of the chlorophyll/histidine and chlorophyll/asparagine

pairs,

we

introduced

coordination

bonds

between

magnesium and the coordinating nitrogen atom of the histidine/asparagine moieties. Such an approach is used in standard CHARMM to describe the coordination of a heme molecule by a histidine sidechain. For simplicity, all bonded parameters required to describe the coordination of chlorophyll by histidine or asparagine moieties were taken from the corresponding CHARMM values for heme coordinated by a histidine sidechain. For the 6 chlorophyll molecules that in the crystal structure are coordinated by water molecules, no additional parameters were introduced, i.e., there are no coordination bonds between chlorophyll and water molecules. We used the ParamChem web software special

lipids

of

photosystem

II

39

to assign force-field parameters for

(digalactosyldiacylglycerol,

DGDG;

phosphatidylglycerol, PG; monogalactosyl-diglyceride, MGDG; and sulfoquinovosyl diacylglycerol, SQDG), and for the bicarbonate and β-carotene molecules. The complex composed of bicarbonate, iron, and four histidine moieties is located at the stroma side of photosystem II, i.e., far away from the region of interest for proton transfer (Fig. 1C); to maintain the structural integrity of this complex, we introduced coordination bonds with associated force constants of 100 kcal/mol · Å2 between iron, the bicarbonate carboxylate oxygen atoms, and the histidine Nε atoms. Likewise, to preserve the internal geometry of the manganese cluster, we introduced bonds with a force constant of 500 kcal/mol · Å2 between all members of

ACS Paragon Plus Environment

9

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 55

the manganese cluster, and valence angles with force constant of 200 kcal/mol · Å2. As reference values for the bond lengths and valence angles we used the crystal structure from PDB ID: 4UB6, ref. 3. We initially assigned partial charges to the atoms of the manganese cluster according to the oxidation state proposed in ref.

3

for the S1 oxidation state, with

Mn1(III), Mn2(IV), Mn3(IV) and Mn4(III). It had been suggested that one of the oxygens of the OEC (O5) is protonated 3; for simplicity, all five oxygen atoms of the OEC were described as dianion oxygens, each with a charge of -2. The calcium ion was assigned a charge of -2. On the basis of test molecular dynamics simulations of photosystem II, whereby we assessed interactions between the manganese cluster and its immediate protein environment, we reduced by one unit the absolute value of all charges within the manganese cluster. These reduced charges were then used in all simulations presented here. As discussed in the Results section, the description of the manganese cluster used here, albeit somewhat simplistic, allows the environment of the manganese cluster to be stable.

Protonation states. We use molecular dynamics simulations to probe how local H bonds of photosystem II respond to changes in the protonation state at the D1E65/D2-E312 site (Fig. 1B, Fig. 1C, Fig. 1D, Fig. 1E): In the high-resolution crystal structure the carboxylate groups of D1-E65 and D2-E312 are within H-bond distance, suggesting that this carboxylate pair might be protonated computations have suggested that D2-E312 might be protonated

10

, and pKa

9, 15

. For simplicity,

we assigned standard protonation states to all titrable groups in the protein complex, that is, Asp and Glu sidechains are negatively charged, Arg and Lys sidechains are

ACS Paragon Plus Environment

10

Page 11 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

positively charged, and His sidechains are singly protonated on the Nδ atom. Exceptions here are the histidine sidechains CP43-H157 and CP47-H201, which coordinate chlorophyll-a molecules at the Nδ atom (and are therefore protonated on the Nε atom), and the glutamate pair D1-E65/D2-E312, for which we considered three different protonation states as follows. We tested the likely protonation state of the D1-E65/D2-E312 pair by performing two independent simulations in which either D1-E65 or D2-E312 was protonated, and a third simulation in which both carboxylate groups were negatively charged (Table 1).

Molecular dynamics simulations. We performed all-atom molecular dynamics simulations using NAMD

40-41

with the smooth particle mesh Ewald summation for

the Coulomb interactions

42 43

, a switch function with cutoff at 12 Å for the van der

Waals interactions. Following geometry optimization, we used velocity rescaling with a frequency of 500 fs during an equilibration phase in the NPT ensemble (temperature T = 300K, pressure P = 1 bar) with an integration step of 1 fs and slow release of weak harmonic constraints as follows. During the first 50 ps we placed harmonic constraints of 3 kcal/mol · Å2 on atoms of the protein backbone and of cofactor molecules, and constraints of 2 kcal/mol · Å2 on protein sidechains. During the next 100 ps we decreased the constraints to 1.5 kcal/mol · Å2 on the protein backbone and cofactors molecules, and switched off all constraints from sidechains atoms. Finally, we equilibrated the system for 1ns without harmonic constraints. The simulations then proceeded with the production run, which was performed in the NPT ensemble using a Langevin dynamics scheme piston

45

and a multiple time step algorithm

46-47

44

and a Nose-Hover Langevin

with bonded interactions computed

each 1 fs, van der Waals every 2 fs, and PME every 4 fs.

ACS Paragon Plus Environment

11

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 55

Molecular dynamics to probe fast water motions. H-bonding between water molecules can be associated with short, sub-picosecond timescales

12, 48-50

. In order

to probe fast motions of water molecules that could participate in H-bonding networks, we augmented the NPT simulations described above by simulations in the NVE ensemble (constant volume V and constant energy E). Starting from five different coordinate snapshots taken from Sim1 at 128 ns, 132 ns, 136 ns, 140 ns and 144 ns respectively we performed, for each snapshot, 5x100 ps independent NVE simulations by saving coordinate sets every 10 fs.

Analysis of protein/water H-bond networks and definition of water Wires. Analyses of dynamic H-bonds in photosystem II pose the challenge of identifying the mobile water molecules that could transiently bridge groups of photosystem II whose sidechains can H-bond. In order to make these analyses computationally efficient, we developed an algorithm to detect water wires connecting pairs of amino acid residues. As summarized in Scheme 1, the basic idea of the algorithm is to map, onto the protein H-bond sidechains, chains of H-bonding waters that have been precomputed, and then subject this ensemble of protein sidechains and H-bonded water chains to a search criterion which ensures that each pair of protein H-bond sidechains is bridged by the shortest water chain. The maximum number of water molecules in such a chain is 5, which corresponds to essentially zero probability of proton transfer between acid/base pairs in water 51. As H-bonding criterion we consider a distance of 3.5 Å between the heavy atoms, and allowed for a deviation of up to 60º of the H-bonding angle. For

ACS Paragon Plus Environment

12

Page 13 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

simplicity, as amino acid residues whose sidechains may be connected by a hydrogen-bonded water chain (i.e., the end points of a Wire) we consider Asp, Glu, and His groups; a Wire can also bridge directly to the OEC. Hydroxyl groups of Ser, Thr and Tyr sidechains may be an intermediate of the Wire if they are directly hydrogen bonded to one of the end points of the Wire. A water Wire is thus defined as two protein groups being connected via the shortest-distance hydrogen-bonded water chain of 1-5 water molecules; a Wire can also connect a protein group to the OEC. The analysis of H-bonds proceeds in three steps (Scheme 1). In the first step, we compute all water–water, water–sidechain, sidechain–sidechain H-bonds. The output of this calculation is used to build up a graph of H-bonding water molecules (Scheme 1A). In the second step, for each sidechain si, the algorithm queries the graph for shortest water Wires that H-bond to si (Scheme 1B). In the third step, the algorithm checks whether waters belonging to any of the shortest water Wires detected in step 2 as being connected to si, fulfill the criterion for H-bonding to a second sidechain, sj. In case this criterion is fulfilled, we consider that si and sj are connected via a H-bonding water Wire (Scheme 1C). The protein/water H-bond networks identified with the algorithm described above were then subjected to analyses of their occupancy. We considered that the network is present (i.e., occupied) for a certain time interval (i.e, number of consecutive coordinate snapshots) only if the water Wire of the network is exactly the same: that is, the water wire consists of the same water molecules in the same order. The size of this time interval used for analysis was different in the NPT vs. NVE simulations. In the case of the NPT simulations (Sim1-5, see Table 1), we considered a minimum interval of two consecutive frames. In Sim1-5 coordinates were saved each 10 ps,

ACS Paragon Plus Environment

13

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 55

this implies that we select, as occupied, only water Wires that remain the same for at least 10 ps. Thus, the protocol used cannot describe the dynamics of the network within the 10 ps interval, but it allows us to find out whether, at the beginning and the end of the 10 ps time interval, the water Wire connecting a specific pair of residues is identical or not. If the water wire is identical, we consider that the Wire either did not change within the 10 ps time interval, or it might have undergone small, transient rearrangements. We used the protocol as summarized above to analyze the last 48 ns of each of the NPT trajectories, that is, for each NPT trajectory we used 4,800 coordinate snapshots to analyze the H-bond networks. In the case of NVE simulations, we analyzed the entire trajectories, i.e., 10,000 coordinate snapshots for each NVE simulation, and computed occupancy rates by using a minimum interval of 100 consecutive frames (i.e., a time interval of 1ps). Water molecules that are part of a protein/water H-bond network can be rather dynamic on the picosecond time scale

12

, and fast motions of water molecules may

introduce noise in the analysis of the H-bond bridges. To circumvent this issue, when analyzing protein/water H-bonds based on the NPT trajectories, we introduced an error function whereby, if a shortest path between two residues is broken or increased in length the algorithm checks if the starting shortest path is present in a second graph that is build using larger cut-offs typical of weak H-bonds: H-bonding distances (3.8 Å) and angles (70º). This graph is queried for maximum 10 frames in a row. If the starting shortest path does not meet the H-bond criteria, without the error function being applied, within this frames interval, a new shortest path is selected. Unless specified otherwise, in our analysis we include only water-mediated bridges with an occupancy rate of at least 5%.

ACS Paragon Plus Environment

14

Page 15 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

We denote the number of water molecules in a Wire as the length of that Wire. To calculate the number of water molecules in the first and second hydration shell of the water molecules of a water Wire, we first queried the graph of H-bonded water molecules for the first neighbors of the water molecules involved in a specific water Wire. This gives us all water molecules that are H-bonded (in the first hydration shell) to water molecules constituting the Wire. The same procedure is then performed to identify water molecules in the second hydration shell of water molecules constituting the Wire, this time starting from the water molecules belonging to the first hydration shell. The endurance time of a particular water Wire is defined as the length of time that that water Wire remains composed of the same discreet water molecules during the last 48ns of the simulation. Processing of the raw H-bond data was done using the python MDAnalysis toolkit 52-53

. Molecular graphics were prepared using VMD 54 and PyMol 55.

Results and Discussion During the ~150ns of each of our main simulations, the protein structure remains largely stable. Overall, the average Cα root-mean-square distances (rmsd) calculated relative to the starting crystal structure coordinates remain within ~1.5 Å and ~2.5 Å (Table 1, Fig. S1). These relatively small values of average Cα rmsd likely originate from the fact that a significant part of photosystem II consists of hydrophobic

transmembrane helices

packed against

chlorophyll

molecules

coordinated by histidine or asparagine amino acid sidechains (Fig. 1A, Fig. 1C). For chlorophyll molecules whose magnesium ion interacts with water in the starting crystal structure, we find that the average distance of water–magnesium interactions

ACS Paragon Plus Environment

15

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 55

is well described, such that the average value computed from simulations, 2.06 Å, is very close to the 2.07 Å average distance computed from the crystal structure (Fig. S2). The structure of the environment of the manganese cluster remains stable (Fig. S3).

Local dynamics of photosystem II depends on the protonation state A particularly important site for proton transfer in photosystem II is that of the carboxylate pair D1-E65/D2-E312: In the starting crystal structure

2

the carboxylate

oxygen atoms of D1-E65 and D2-E312 are within 2.5 Å distance (Fig. 1E), raising the possibility that this carboxylate pair could be protonated

10

. Moreover, FTIR

spectroscopy data suggest that D1-E65 and D2-E312 are part of an extended Hbond network that gives rise to a protonated carboxylate vibrational mode 21-22. In the crystal structure, the D1-E65/D2-E312 pair is part of a H-bonded network that includes groups from subunits D1, D2 and PsbO

10

(Fig. 1E). D1-E65 is at the

center of this H-bonded network: its carboxylate group H-bonds to D2-E312 and D2N335, and its backbone H-bonds to D1-D59 (Fig. 1E). D2-R334 bridges via H-bonds to D2-E312 and to the backbone carbonyl groups of PsbO-L157 and PsbO-D158 (Fig. 1E). In total, the H-bonded network of the D1-E65/D2-E312 pair depicted in Fig. 1E includes 13 H-bonds. For simplicity, in what follows we will refer to this Hbond network as the central H-bonded network of photosystem II. As first step towards exploring the coupling between H-bond dynamics and protonation state at the central H-bonded network (Fig. 1E), we performed three independent simulations of photosystem II, in which we consider different protonation states of D1-E65 and D2-E312. In Sim1, D1-E65 is protonated and D2-

ACS Paragon Plus Environment

16

Page 17 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

E312 is negatively charged (Fig. 2A); in Sim2, the protonation is reversed (Fig. 2B, Table 1). Both D1-E65 and D2-E312 are negatively charged in Sim3 (Fig. 2C). We find that, when D1-E65 is protonated (Sim1), no fewer than 11 of the H-bonds observed in the crystal structure remain present. The carboxylate group of D1-E65 makes water-mediated interactions with D1-D61, whereas D2-E312 interacts with D1-R334 (Fig. 2A). The dynamics of H-bonds involving D1-E65, D2-E312 and D1R334 is overall well converged (Fig. 1F), suggesting that these three groups involve in persistent H-bonds within the central H-bonded network (Fig. 2A, Fig. 2D). Sim2, in which D2-E312 is protonated, preserves the H-bond between D1-E65 and D2-E312 (Fig. 2B, Fig. 2D). Nevertheless, only 9 of the H-bonds of the central H-bond network are preserved, and H-bonding of D1-R334 and D1-N335 is perturbed (Fig. 2B). That is, protonated D2-E312 (Sim2) appears to associate with the central H-bond network (Fig. 1E) being overall less preserved than when D1E65 is protonated (Sim1). Sim3, in which both D1-E65 and D2-E312 are negatively charged (Table 1), preserves only 8 of the H-bonds of the central H-bond cluster; D1-E65 rotates towards and interacts with PsbO-R152, whereas D2-E312 interacts with D1-R334 (Fig. 2C). That is, in simulations where the D1-E65/D2-E312 pair carries no proton, important features of the central H-bonded network are poorly described. To further assess the impact of the protonation state on interactions within photosystem II, we analyzed the inter-subunit H-bonds between subunits that either directly participate to the central H-bond cluster, or are located close to this cluster; the subunits included in this analysis are D1, D2, CP43, CP47, PsbO, and PsbU. As summarized in Supplementary Table ST1, the three simulations on wild-type photosystem II (Sim1-Sim3) provide a reasonable description of the inter-subunit H-

ACS Paragon Plus Environment

17

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 55

bonds included in the analysis: Overall, Sim1 preserves 60.2% of inter-subunit Hbonds present in the starting crystal structure, as compared to 53.6% in Sim2, and 56.1% in Sim3. Simulations probing the response of the protein structure and dynamics to changes in the protonation state could, in principle, be used to think about the likely protonation state in the crystal structure, as the protonation state that best preserves structural features of the protein might be considered to be the more likely protonation state. A caveat of such an approach is that details of H bonding in a protein crystal structure might be different when protein motions are studied at room temperature, with the protein embedded in a fluid lipid membrane environment 18. In the analyses summarized above, Sim1 (D1-E65 protonated) associates with overall best description of the H-bonds, but it does not preserve a direct H bond between D1-E65 and D2-E312 (Fig. 2A). Sim2 (D2-E312 protonated), although it preserves the direct H-bond between D1-E65 and D2-E312 (Fig. 2B), preserves overall fewer of the H-bonds observed in the crystal structure. In what follows, we consider Sim1 as a reference for discussing the H-bond networks of wild-type photosystem II. We note, however, that our current simulations and data analyses do not exclude the possibility that D2-E312 is protonated, as suggested in refs.

9, 15

based on crystal

structure analyses.

Dynamic H-bonded networks and water Wires of wild-type photosystem II A key open question that is left open by analyzing H-bonds from static crystal structures alone (Fig. 1B, Fig. 1C, Fig. 1D, Fig. 1E) is that of how dynamic the Hbonds are. This question is particularly important for H-bond clusters because, in a cluster, H-bond sidechains can engage in dynamic interactions whereby H-bonds

ACS Paragon Plus Environment

18

Page 19 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

break and reform, and H-bond partners may change

18, 20

. An illustration of the

dynamics of H-bonds for selected amino acid residues of the central H-bonded network is given in Fig. 1F. Such H-bond dynamics are important, because they can result in proton-transfer pathways being established or interrupted. Analyses of H-bonding in Sim1 allowed us to characterize two H-bond networks that extend from the vicinity of the manganese cluster to the lumen bulk (Fig. 3A, Fig. 3B, Fig. 3C, Fig. 3D). The first of these networks, which we denote as the PsbO network, extends from the vicinity of the OEC to PsbO via four water-mediated Wires that include D1-D61, D1-E65, D2-E310, and D2-E312 (Fig. 3A, Fig. 3B, Table 2). For each of these Wires, we observe more than one possible connection characterized by a different end group and by its average occupancy (Table 2, Table 3). For example, Wire 1 connects D1-D61 to D1-E333, which is part of the coordination shell of the manganese cluster; in Sim1, Wire 1 has an overall occupancy of 67.9% (Table 2, Table 3). But D1-D61 can also connect to D1-D170, or directly to the OEC, via Wires with lower occupancies, of ~20-27% (Table 2, Table 3). For simplicity, these two other Wires that connect D1-D61 to the OEC or its immediate vicinity are labeled Wire1b and Wire1c, respectively (Table 3). The higher occupancy Wire1 is considered as the main Wire1 and used for illustrating the long-distance PsbO H-bond network that reaches from the vicinity of the OEC to PsbO (Fig. 3A). Likewise, D2-E312 has a main connection - Wire3 - with D2-E310 (50.1% occupancy, Fig. 3A, Fig. 3B, Table 2, Table 3) further bridges via H-bonded water to PsbO-D222 (Wire3a, ~33% occupancy) and to PsbO-D224 (Wire3b, ~45% occupancy, Table 3).

ACS Paragon Plus Environment

19

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 55

Water-mediated bridges between D1-D61 and D1-E65, and between D2-E310 and PsbO-D224 (Wire2 and Wire4 in Fig. 3B, Table 2), have an occupancy rate higher than 90%. The two other wires, Wire1 and Wire3, are more dynamic, being present for ~50-68% of Sim1 (Fig. 3B, Table 2). The second H-bonded network, denoted here as the CP43 network (Fig. 3B), starts with D1-E329 and it extends, via five intermediate wire segments, to CP43E104 (see Wire5, Wire6, Wire7, Wire8 and Wire9 in Table 2). Compared to the PsbO network (Fig. 3A, Fig. 3B), the CP43 network (Fig. 3C, Fig. 3D) is overall more dynamic, the highest occupancy of an intermediate wire being ~74% (Wire6, see E329-E413 in Fig. 3D). D1-E329 also establishes a watermediated H-bond with D1-H190 (12.6%), which is constantly H-bonded to D1-Y161. Both D1-H190 and D1-Y161 are linked to the manganese cluster via water Wires (Fig. 3C). The fact that D1-Y161 (also denoted as Yz) and D1-H190 are part of the CP43 H-bond network is compatible with experiments suggesting that these two amino acid residues could participate in a proton exit pathway 56-58. The analysis above suggests a complex molecular picture of the internal Hbonded networks of photosystem II, whereby long-distance H-bond networks between the vicinity of the OEC and the lumen consist of several segments that can have largely different occupancies (Fig. 3, Table 3). Moreover, alternative, lowoccupancy connections may be sampled (Table 3). This dynamic picture of the water-mediated H-bonds between protein groups of the long-distance H-bond networks would have been difficult to predict from the crystal structure.

ACS Paragon Plus Environment

20

Page 21 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Dynamics of the water Wires To characterize the dynamics of water Wires we used Sim1 to calculate, for each of the main Wires from the PsbO and CP43 H-bonded networks (see Table 2, Table3), the average length of the Wire, the number of waters in the first and second hydration shell of the Wire waters, and estimates of the endurance time of the Wires. Results of these computations are illustrated in Table 2, Fig. 4, Fig. S4, Fig. S5, Fig. S6, and summarized below. Overall, we observe that shorter Wires tend to have higher average occupancies (Table 2). Wire2 and Wire4, which have average occupancies of ~94-96%, are short Wires in which carboxylate groups bridge via one water (Fig. 4B, Fig. 4D, Fig. S4B, Table 2). A somewhat lower, though still significant (~68%), occupancy is observed for Wire1, which is also a one-water Wire (Table 2), but there is infrequent sampling of conformations in which Wire1 includes two H-bonded waters (Fig. S4A). By contrast, Wire3 (Fig. 4A), which has an average occupancy of ~50% (Table 2), consists most of the time of three H-bonded waters, but there is noticeable sampling of conformations in which Wire3 consists of two or four waters (Fig. 4C). The three lowest-occupancy Wires, Wire5, Wire7 and Wire9 (with average occupancies