Investigating the Many Roles of Internal Water in Cytochrome c

Publication Date (Web): July 16, 2018. Copyright ... Moreover, Dowser++ predicts many more internal water molecules than is typically seen in the expe...
0 downloads 0 Views 1MB Size
Subscriber access provided by - Access paid by the | UCSB Libraries

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

Investigating the Many Roles of Internal Water in Cytochrome C Oxidase Ardavan Farahvash, and Alexei Stuchebrukhov J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b11920 • Publication Date (Web): 16 Jul 2018 Downloaded from http://pubs.acs.org on July 18, 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 37 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

Investigating the Many Roles of Internal Water in Cytochrome c Oxidase Ardavan Farahvash and Alexei Stuchebrukhov* Department of Chemistry, University of California – Davis, One Shields Avenue, Davis, CA, 95616 *Email: [email protected]

ABSTRACT Cytochrome c oxidase (CcO) is the terminal enzyme in the respiratory electron transport chain. As part of its catalytic cycle, CcO transfers protons to its Fe-Cu binuclear center (BNC) to reduce oxygen and, in addition, it pumps protons across the mitochondrial inner, or bacterial, membrane where it is located. It is believed that this proton transport is facilitated by a network of water chains inside the enzyme. Here we present an analysis of the hydration of CcO, including the BNC region, using a semi-empirical hydration program, Dowser++, recently developed in our group. Using high resolution X-ray data, we show that Dowser++ predictions match very accurately the water molecules seen in the D- and K- channels of CcO, as well as in the vicinity of its BNC. Moreover, Dowser++ predicts much more internal water molecules than is typically seen in the experiment. However, no significant hydration of the catalytic cavity in CcO described recently in the literature is observed. As Dowser++ itself does not account for structural changes of the protein, this result supports the earlier assessment that the proposed wetting transition in the catalytic cavity can only be either due to structural rearrangements of BNC, possibly induced by the charges during the catalytic cycle, or occurs transiently, in concert with the proton transfer. Molecular dynamics simulations were performed to investigate the global dynamic nature of Dowser++ waters in CcO, and the results suggest a consistent explanation as to why some predicted water molecules would be missing in the experimental 1 ACS Paragon Plus Environment

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 37

structures. Furthermore, in the light of the significant protein hydration predicted by Dowser++, the dielectric constant of the hydrated cavities in CcO was also investigated using the FröhlichKirkwood model; the results indicate that in the cavities where water is packed sufficiently densely the dielectric constant can approach values comparable even to that of bulk water.

Introduction The mechanism of oxygen reduction by cytochrome c oxidase and its connection to proton pumping has been a topic of intense research since the discovery of the enzyme’s proton pumping in 19771. During its catalytic cycle CcO consumes 8 total protons from the (negative) N-side of the mitochondrial inner, or bacterial, membrane where it is located; four of these protons are taken to the binuclear Fe-Cu center (BNC) where they are used for the reduction of oxygen, and another four are pumped across the membrane2-3 to the (positive) P-side, + . 4e- + 8Hin+ + O2  → 2H2O + 4Hout

(1)

Two hydrophilic input channels, named the K and D channels, are responsible for the conduction of protons to the BNC4-5. The D-channel is believed to conduct all pumped protons through a roughly half of the membrane depth, from the N-side up to the lower boundary of BNC region. However, the further path of the pumped protons in and around BNC, as well as their exit pathways (the output channels) are less clearly defined, and were primarily addresses in computational studies6-7. Several mutational studies have shown that the alterations of the residues in the input D- and K- channels lead to either a significantly lower enzyme activity, a complete inhibition, or a selective loss of the pumping4, 8-14, depending on the type of mutation. The conceptual understanding of the pumping mechanism, as well as its catalytic chemistry at BNC, has emerged in the past two decades2-3,

13, 15-23

, however, the detailed molecular

mechanism still remains elusive, primarily due to the difficulty in monitoring detailed proton 2 ACS Paragon Plus Environment

Page 3 of 37 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

motions in the enzyme. The measurements of the potential generated by the enzyme upon single electron injections, in combination with site-directed mutations, has been most fruitful technique so far, generating most of data on the performance of the pump24-26. One of the key theoretical difficulties is the ill-defined internal water in the enzyme that facilitates proton transport20-21,

27-32

, and more specifically its metastable dynamic nature and

apparently transient behavior along the catalytic cycle32. of the enzyme. In addition, two new water molecules are produced in BNC in every cycle of the enzyme, and constant water transport occurs inside the enzyme, which is still poorly understood. Even the simplest equilibrium property, i.e. the mere question of how much water is actually present in the enzyme is still unsettled33. Recent computational studies have advanced our understanding of CcO, due to both more accurate simulation of water in the enzyme, and the accurate calculations of proton energetics 2022, 31-32, 34-39

. Two new theoretical ideas involving internal water were particularly useful in

advancing the mechanism of CcO, namely the demonstrated fluctuational nature of the water chains in the BNC region31, 40, and a possible reversible wetting transition that is proposed to occur during the catalytic cycle of the enzyme21,

32

. The unusual properties of water in the

constrained hydrophobic environment, clarified only recently41, plays the crucial role here. While monitoring proton motions presents the major experimental challenge, experimental studies of internal water in CcO are also non-trivial, and the predicted new dynamic roles of internal water remain experimentally unconfirmed. It is known that not all water molecules are resolved in the crystal structures even at the highest resolutions42-43. It is also not clear how much the (usually low-temperature) crystal structures correspond to the enzyme natural conditions; this is particular true for equilibrated internal water, as the environment of the

3 ACS Paragon Plus Environment

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

enzyme in the crystal is different from that of the active state enzyme in its natural conditions in the membrane, and the partial “drying out” of some of the internal water cannot be excluded15. Characteristically, no stable water chains have been observed in the hydrophobic catalytic cavity of BNC (roughly above Glu242, see below Figure 2) in all experimental x-ray crystal structures taken to date. Many of such structures are of sufficient resolution to expect the detection of such water molecules44. Overall, however, the scarcity of water molecules in the core of the enzyme seems to agree with our understanding of the structure of cytochrome c oxidase, as the enzyme’s residues rapidly turn more hydrophobic as we approach the catalytic BNC cavity, Figure 1. In this sense the sparsity of water around BNC is natural and justified.

Figure 1. (A) The structure of cytochrome c oxidase (PDB code: 1V54). The catalytic BNC is located near the center of the high spin heme. The z-axis shown is the same coordinate axis used to construct B and C. (B) The hydrophobicity of the standard protein residues calculated as a function of their depth in the membrane. The Whimley-White hydrophobicity scale45 was used to calculate individual amino-acid hydrophobicities and these hydrophobicities where averaged as a function of their corresponding z-coordinate. (C) The density of buried water molecules in cytochrome c oxidase plotted as a function of the depth in the membrane. In both B and C, 20 equally

4 ACS Paragon Plus Environment

Page 4 of 37

Page 5 of 37 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

spaced rectangular cross sections corresponding a roughly 5Å slice were used to bin amino-acids/water molecules.

On the other hand, water is absolutely vital to the chemical function of the BNC, and one would expect there should be a plethora of water molecules in that part of the protein. Yet, such water molecules are distinctly sparse in experimental structures, as seen in Figure 1. Predicting internal water in proteins has been a long standing challenge for computational research33, 46. A number of techniques for free energy evaluation of the inserted water molecules, that includes proper configurational sampling of the protein, have been advanced recently32-33, 38, 46

. However, simple protein hydration programs, such as Dowser, remain to be very popular in

the research community for their efficiency and simplicity of use. Recently, we have developed an advanced version of Dowser, which is called Dowser++47. We have verified (see Methods section) the performance of Dowser++ by testing it on a set of sub-1Å resolution structures in the protein databank, and such tests have shown that Dowser++’s water predictions indeed correspond to waters seen in super high resolution crystal structures. Here we describe an application of Dowser++ to study the internal water in CcO. The objective is three-fold. First, we explore the global solvation of CcO – i.e. to address the question of how much (equilibrated) water is inside the enzyme, and to characterize its dynamic properties. Of particular interest are water molecules that are predicted but not seen in the x-ray structures. Dowser++ predicts a surprising plethora of such water molecules. Molecular dynamics simulations show that these additional Dowser++ predicted water molecules are characteristically more dynamic than the experimentally reported ones, which provides a possible explanation as to why these functionally important water molecules are missing in the xray structures. Second, we specifically analyze the hydration predictions of Dowser++ in the regions of D- and K-channels, and in the catalytic cavity of BNC, where many previous studies have 5 ACS Paragon Plus Environment

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

focused. As in the global hydration picture, we observe some additional water molecules in both D- and K- channels, compared to those reported in the experiment, and in some previous computation studies33; at the same time, we demonstrate a remarkable accuracy with which Dowser++ reproduces those water molecules that are reported in the high-resolution X-ray structures. Despite the prediction of many new water molecules in the protein overall, our calculations failed to reproduce the recently described wetting transition of the hydrophobic catalytic cavity either upon stable protonation of propionate D of heme a3, as in Son et al.21, or upon its transient protonation during the proton transfer through this group32. As Dowser++ does not include structural changes or configurational sampling of the protein (it works with the original static crystal structure of the protein), we interpret this result as supporting earlier analysis21 that the proposed transition is due to the structural rearrangements of the catalytic cavity of CcO, induced by the introduced charge, rather than the due to only additional electrostatic interactions between the introduced charge and water molecules in the cavity. Alternatively, the wetting transition may occur transiently, together with the proton transfer, as suggested in Liang et al.32 Of course, our water placement program is unable to capture either of these effects, as it works with only equilibrium hydration. Third, given the large number of internal water molecules predicted by our program, we calculated their contribution to the dielectric constant of the highly hydrated cavities in CcO. In principle, the core of CcO should be a relatively low dielectric region, with ε ≈ 4 − 8 , however in regions where water molecules are closely packed, the dielectric constant may differ drastically from the typical dry protein values48-49. To test this assumption we used K-means clustering50 to identify groups of tightly packed water molecules inside CcO, and then extracted the dielectric constant from the fluctuations in the dipole moment of individual water molecule using Fröhlich-

6 ACS Paragon Plus Environment

Page 6 of 37

Page 7 of 37 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

Kirkwood theory51.

Methods 1. Dowser++ Predictions Dowser++ is a protein cavity hydration predictor47 based on a combination of two other hydration programs Dowser and WaterDock. Dowser was released by Zhang and Hermans in 199652. The program requires a PDB file with atomic coordinates of a protein as input, and as output it produces another PDB file with a list of coordinates of placed water molecules together with their binding energies. The binding energy of the waters is evaluated by using a specific model of water-protein interaction which includes electrostatics as well as 6-12 Lennard Jones interactions, and only molecules with binding energy greater than a threshold value of 10 kcal/mol are retained. To improve the accuracy of Dowser, two significant modifications have been made. First, a new energy-evaluation scheme was incorporated. This new scheme follows from a physical model of the expected change in the solvation energy upon a transition of a water molecule from the bulk phase to a protein medium53. Following the model, the maximum binding energy was also decreased from 10 to 4 kcal/mol. It should be noted that no explicit entropic contributions are included in Dowser/Dowser++ energy criterion, however we expect this contribution to be rather small for a static protein structures (~1 kcal/mol)53. Dowser++ works with a static protein structure taken from an x-ray data, and therefore cannot capture possible hydration effects due to variations of the structure that may transiently occur during the catalytic cycle of the enzyme, or occur transiently, concerted with the proton transfer 21, 32. Its predictions also only correspond to the equilibrium hydration state of a protein. However, the advantage Dowser++ comes from the 7 ACS Paragon Plus Environment

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

fact that it can hydrate the entire protein without making any reference to the existing crystal structure waters, unlike many free-energy hydration methods. The accuracy of such predictions will be discussed below in this paper. The second change made to Dowser was an improvement on its search algorithm used to find internal hydration sites in the protein. Dowser uses a simple grid based search to determine potential hydration sites. This algorithm limits Dowser’s predictive capability as the grid can only be made so fine without risking placing water molecules on top of each other. In order to improve Dowser’s internal search algorithm Dowser++ uses WaterDock, a tool originally designed to predict the locations of water molecules in protein-ligand binding sites54. WaterDock alone is insufficient to hydrate a whole protein as it is a machine learning based program that has not been calibrated for such a task. However, WaterDock’s use of the modern Monte-Carlo based docking techniques of Autodock Vina55 makes it very efficient at searching the interior of a protein for potential water binding sites. Dowser++ utilizes WaterDock’s search algorithm, and refines its binding energy calculation with the improved 4 kcal/mol energy criterion. Additionally, parameters for the heme, heme-ligands, and metal centers of CcO were added to Dowser++’s standard parameter set to improve prediction accuracy on CcO. In order to evaluate the quality of Dowser++ predictions against experiment, a set of sub1Å crystal structures with sizeable internal cavities were selected from the PDB. The accuracy of Dowser++ on this set was determined by comparing the positions of predicted water molecules to experimental ones. A prediction was considered a successful match if the positions matched within 2.8Å (the van der Waals diameter of water). Results are shown in Table 1.

8 ACS Paragon Plus Environment

Page 8 of 37

Page 9 of 37 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

Table 1. Performance of Dowser++ (D++) on a set of sub-1Å resolution PDB structures. % Exp matched refers to the percent of experimental water molecules within 2.8Å of a water placed by Dowser++. PDB Code Res (Å) Exp waters D++ waters % Exp matched 1us0

0.66

22

19

72.7

4rek

0.74

80

84

86.2

1gci

0.78

12

16

100

4ua6

0.79

45

42

77.8

3d43

0.8

37

47

94.6

2j8t

0.82

21

33

71.4

4y9w

0.83

16

25

100

4ayp

0.85

39

46

92.3

4igs

0.85

15

15

86.7

4xxg

0.85

97

76

77.3

4uaa

0.86

44

40

61.4

2p74

0.88

48

43

75

As seen in Table 1, the accuracy of Dowser++ based on the comparison between the number of water molecules predicted and the number seen in the experiment is reasonably high. Through a widespread analysis of the entire PDB, Carugo and Bordo43 have shown that the number of water molecules seen per residue in a crystal structure has a high correlation (dramatically increases) with the resolution the structure was taken at. One of the best resolution structures of the entire CcO molecule taken to date is 1.8Å (PDB Code: 1V54). This is significantly lower quality than the sub-1Å proteins used in Table 1 to verify the performance of Dowser++. Thus we expect much more water molecules would be predicted than seen in such a crystal structure. It is also worth noting that the highest resolution crystal structures of bovine 9 ACS Paragon Plus Environment

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

cytochrome c oxidase (including 1V54) are all taken at temperatures at or below 100K, therefore the water occupancy seen in these structures may not be perfectly indicative of the occupancy of the enzyme at room temperature, at which simulations were performed in our study.

2. Molecular Dynamics MD simulations were used to explore the dynamic nature of the buried water molecules predicted by Dowser++ and compare them to those found in the native crystal structure. The entire thirteen subunit bovine enzyme (PDB: 1V54) was used to construct the simulation system, and additionally the enzyme was placed in an explicit membrane environment of POPC lipids. CcO was properly positioned in the membrane following data from the OPM (Orientation of Proteins in Membranes) database56. Sodium chloride ions were added to neutralize the system and raise the ionic concentration to 0.100M. All simulations were performed with NAMD using the CHARM 27 force field and the TIP3P model for water molecules44, 57-58. Simulations were done in a NPT ensemble using a Berendsen barostat and Langeven dynamics (pressure = 1atm, temperature = 312K), and periodic boundary conditions were employed with full periodic electrostatics. Hydrogen atoms were held rigid to allow for a 2 femtosecond timestep. A cutoff distance of 12 Å and a switching distance of 10 Å was applied for nonbonded interactions. Electrostatics were evaluated with the particle mesh Ewald (PME) method implemented in NAMD. An initial energy minimization of 5000 steps was performed before each simulation. Simulations were performed for a total 5 nanoseconds, with atomic coordinates output every 2 picoseconds. Heme, heme-ligand, and metal center forcefield parameters, including charge parameters, were adopted from Johansson et al59. These parameters have been widely used and tested in the literature60-62. In order to probe the sensitivity of our predictions to various charge states of the

10 ACS Paragon Plus Environment

Page 10 of 37

Page 11 of 37 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

enzyme, two different states – fully reduced and fully oxidized – were selected. It should be noted that these states are likely not corresponding to any specific intermediate (and metastable) redox/protonation states along the enzyme cycle2. These latter states are not clearly defined, and still debated in the literature; in any case, their study is beyond the scope of this paper. The fully reduced state corresponds to the following formal charges: CuB(I), Fea3(II), Fea(II), and CuA(II); the fully oxidized state corresponds these formal charges: CuB(II), Fea3(III), Fea(III), and CuA(III). These states were also chosen as they correspond to the ones most commonly found in the x-ray crystal structures. A single water molecule ligand was placed on the CuB center of both states, as Dowser++ was found to predict water molecule in this position in both the fully reduced and fully oxidized states. Glu 242, a critical residue in the proton transport and proton gating of CcO, was modeled as protonated and initially placed in “up” configuration without constraint in all simulations. Asp 364, Asp 369, Lys 319, His 207, and Tyr 244 were all kept protonated in both states in correspondence with available experimental and theoretical studies59, 63-64

. While some recent reports suggest23,

64

that in the oxidized enzyme the tyrosine may be

deprotonated, we believed it was best to keep it protonated in order to avoid introducing an unnatural imbalance of charge, and to avoid the issues of specific metastable states of the enzyme along its cycle65. All other residues were kept in their standard protonation states. Despite some prior studies indicating that molecular charges should be scaled in accordance with missing electronic polarizability of the protein medium66-67, we have chosen not to scale charges here in order to keep our MD results consistent with previous literature on CcO, and focus only on the differences arising from the presence of water molecules.

3. Dielectric Constant Calculations The Fröhlich-Kirkwood, or Fröhlich-Kirkwood-Onsager (FKO) theory of dielectric

11 ACS Paragon Plus Environment

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 37

polarization considers in a self-consistent manner a dipole moment surrounded by a continuous dielectric medium51. From this treatment we then have: 2 4π δ p (ε − 1)(2ε + 1) . = 3 v0 k BT 3ε

(2)

Equation 2 relates the dielectric constant of the medium ε , to the fluctuations of a dipole moment

p ; v 0 represents the volume of the dipole moment within the dielectric medium. The left-hand side of this equation is often referred to as the Kirkwood G-factor. Equation 2 was used to calculate the dielectric environment of individual water molecules in cytochrome c oxidase. In order to calculate the dielectric constant of water clusters, the result for each individual water molecule belonging to the cluster was averaged. In order to assign water molecules to clusters, a standard K-means clustering algorithm50 was employed with 18 clusters. The number of clusters was selected using the Bayesian Information Criterion68, but was capped so that no clusters were formed containing fewer than 3 water molecules. Additionally, any cluster with fewer than 7 water molecules was excluded from further analysis for being too small. Each cluster was assigned a molecular volume v0 by taking the average distance between each water molecule and its four closest neighbor, averaging that result over all water molecules belonging to a cluster, and cubing it. Dielectric constant calculations were carried out on Dowser++ water molecules simulated in the fully oxidized state of the enzyme, the simulation details of which are elaborated on in the previous section. In order to correlate these dielectric properties of internal water with those of bulk water, calculations were also done in the same manner as described on bulk water, using a 60Å box with periodic boundary conditions containing 6845 water molecules and 5 nanosecond simulation. 12 ACS Paragon Plus Environment

Page 13 of 37 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

Results and Discussion 1. Dowser++ Prediction Analysis on CcO Dowser++ reproduces most of the molecules seen in x-ray structures as well as placing many new ones. In total Dowser++ predicts 216 buried water molecules, while only 122 molecules are seen in experiment (a difference of roughly 80%), and these predictions match 85% of the experimental waters when subject to the 2.8 Å distance criterion. This prediction accuracy is roughly equivalent to the average accuracy on the sub-1Å resolution set, however the total amount of new water molecules predicted is now significantly higher on the account of the lower resolution CcO structure (PDB code 1V54, resolution 1.8Å). Figure 2 (see also Figures S1-S4 in the SI) shows a comparison between the water predicted by Dowser++ and the crystal structure water molecules in the fully oxidized state of the enzyme both near the BNC and in the K and D channels. Neither the positions of the crystal structure waters molecules nor the positions of Dowser++ predictions change significantly upon the transition to the fully reduced state. In general, Dowser++ finds a few additional unseen water molecules providing increased connectivity throughout both the D and K channels. However, only two water molecules were predicted between Glu 242 and the BNC (in addition to one bound between CuB and Fea3 in BNC); these water molecules were spaced too far apart to support nominal hydrogen bonding ( >4Å). This result shows that thermodynamically stable, tightly bound water chains, expected to be there in order to conduct protons, are not formed even with an improved quality of prediction of Dowser++. This is in line with some earlier conclusions 30. We thus further support the forming literature consensus that the water chains in this region are of metastable (fluctuational) character19-20,

30, 40

. Free energy calculations show

that the formation of such water chains may be the result Grotthuss shuttling of water molecules 13 ACS Paragon Plus Environment

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 37

by an excess proton and/or a hydration transition that occurs due to the protonation of propionate D21, 32, 69. Alternatively, newly formed water molecules at BNC could also transiently contribute to these chains 28. It is also of interest to note that the hydration of the K channel predicted by Dowser++, although has additional water molecules compared to those seen in the X-ray structures, does not change in any significant way upon the transition from the fully oxidized to fully reduced state. In both states we kept the Lys that defines the K channel protonated. Thus the possible difference between the resting oxidized state “O” and the functional so-called high-energy “Oh”

70

is not

due to water structure in the K channel, which is quite stable as shown here, but likely has a different origin. To characterize the dynamic behavior of the predicted water molecules - admittedly too few to form stable chains in BNC - we still sought to quantify the persistence of their connectivity using a weaker hydrogen bonding criterion, i.e. putting nearest distance around 4Å, for water chains connecting Glu 242 to the BNC (Chem), Glu 242 to propionate D (on the way to Pump site PLS), and Thr 365 to the BNC (K Channel) in our simulations. In the K-channel, where comparison possible, the Dowser++ waters form far more persistent chains than those formed with only the crystal structure water molecules. Overall, however, in the region between Glu and BNC/PLS predicted waters are too few to form stable well-connected chains. Also, it is difficult to make any definitive conclusions about the persistence of these water chains by only observing 5 nanosecond trajectories used in the simulations. At longer times we expect to see fluctuations in the number of water molecules in BNC; however here we were concerned only with the dynamics of the predicted water molecule, thus short trajectories were sufficient and justified here.

14 ACS Paragon Plus Environment

Page 15 of 37 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

A more detailed analysis of internal water predicted by Dowser++ in D- and K-channels and in the hydrophobic cavity of BNC, along with its dynamic properties and comparison with other studies, are given in the Supplementary Information (SI). Overall, we observe quite remarkable performance (given the simplicity of the method) of Dowser++ in reproducing all 9 experimental water molecules in the D-channel, and those seen in the K-channel, while predicting a few more molecules that are not reported in the experimental structures. In SI we go one-by-one with each predicted water molecules in the D-channel, and discuss comparison with experiment and other studies. Despite many additional water molecules predicted overall, see below, the catalytic cavity remains poorly hydrated by Dowser++ - predicting no more than three water molecules in the equilibrium state of either reduced or oxidized enzyme. An additional proton charge on propionate D, as in Son et al21 and Liang et al32 does not change the situation, and we do not observe the wetting transition of the BNC cavity. The failure of Dowser++ to reproduce such a transition is obviously related to its neglect of possible changes of protein structure that can occur in response to introduced charge. Indeed such a change of structure and related increased volume of the cavity was observed in Son et al21. Such subtle effects obviously require more sophisticated treatment than Dowser++ can offer in its present form. Some additional details and comments related to this issue can be found in the SI. In the following we will focus on the analysis of the global hydration of CcO predicted by Dowser++, which appear to be more suitable for this program. In global hydration, in addition to predicting almost 80% additional water molecules overall, the density of water predicted by Dowser++ in the D-channel is almost 40% greater (3 additional waters) than what is seen in the experiment, and in the K-channel it is twice (4 additional waters) what is seen in experiment. This difference is made even more apparent by the

15 ACS Paragon Plus Environment

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

z-coordinate water density plot (Figure 2B) which shows that, particularly near center of the enzyme, Dowser++ predicts a significantly higher water density (but probably still not enough to conduct protons) than what is observed in the experimental crystal structure. Such substantial differences is address in the next section.

Figure 2. (A) Comparison between predictions of Dowser++ and experimental crystal structure water molecules (PDB Code: 1V54) in the D and K channels and near the Heme-a3. This PDB structure corresponds to fully oxidized enzyme. Dowser++ predictions are shown as red and white water molecules while crystal structure waters as shown as blue spheres. (B) Comparison between the z-coordinate water 16 ACS Paragon Plus Environment

Page 16 of 37

Page 17 of 37 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

density of Dowser++ predictions and experimental water molecules. The same z-coordinate was chosen as is illustrated in Figure 1A. A detailed analysis of water in D- and K-channels, as well as in the catalytic cavity of BNC, is given in SI.

2. Dynamics of Predicted Water Molecules The ability for crystallography to resolve any section of the protein depends on the dynamic stability of that section42, 71-72. If part of a protein exhibits a high degree of dynamic disorder, the electron density for that part will be fragmented or smeared, making it impossible to fit it to an atomic model. The same principles apply to waters inside a protein as well. Highly dynamic water molecules without unique energy minima are very difficult to resolve as their electron density is usually quite smeared42. Thus if the water molecules predicted by Dowser++ were characteristically more dynamic than the water molecules found in experiment, it would explain why they are not seen in experiment. And if not it would at least falsify Dowser++’s prediction, indicating a better protein-water interaction model was needed. We have used molecular dynamics simulations to investigate the dynamic nature of the Dowser++ predicted water molecules in comparison to the experimental molecules in both the fully oxidized and fully reduced states of the enzyme. Figure 3 outlines the mean squaredisplacement (MSD) of the water molecules from these simulations. These MSD were linearly fit in order to calculate the associated diffusion coefficients (Table 3). It is crucial to note that only water molecules that remained within the surface of CcO during the entire simulation where included in these calculations. In both the fully oxidized and the fully reduced states the water molecules predicted by Dowser++ are significantly more dynamic than the experimental molecules. Across the entire protein the Dowser++ predictions had a 120% larger diffusion in the reduced state and a 25% larger diffusion coefficient in the oxidized state. Even more importantly, the set of Dowser++ 17 ACS Paragon Plus Environment

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

predictions that did not agree with experimental molecules, notated by the green D++ - Exp curves, had a diffusion coefficient which was 185% larger in the reduced state and 75% larger in the oxidized state.

Figure 3.

Comparison of MSD between Dowser++ and Experiment. Comparisons were drawn both

amongst all buried water molecules and those only within 10Å of the two heme-groups (i.e. near the catalytic cavity). The green D++ - Exp curves represents the set of Dowser++ waters which do not match with experimental molecules within a distance of 2.8Å. The black dotted lines represent the lines of best fit for each respective colored curve.

If instead of analyzing waters globally across the entire protein, comparisons were only drawn between Dowser++ and crystal structure waters within 10Å of the catalytic heme groups3, the differences become even more pronounced. Table 2 illustrates how, in this hydrophobic core of CcO, the new, D++ - Exp, waters have 4 times the diffusion coefficient of the experimental 18 ACS Paragon Plus Environment

Page 18 of 37

Page 19 of 37 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

molecules in the reduced state and 2.8 times the diffusion coefficient in the oxidized state. It should be noted that this set of Dowser++ molecules include the previously unseen water chains connecting to the BNC, and that such a significant difference in the dynamic nature of the unseen water molecules, as shown by the diffusion coefficients, indeed may explain why they are not seen in experiment.

Table 2. Comparisons of diffusion coefficient (D) given in (cm2s-1·107) between Dowser++ predicted water molecules and experimental crystal structure water molecules. The diffusion coefficient was extracted from the linear regression of the MSD curves shown in Figure 2. D / (cm2s-1·107) Red – All H2O

Red – H2O near Heme-a3

Ox – All H2O

Ox – H2O near Heme-a3

Exp

1.04

0.73

2.04

1.60

D++

2.31

1.93

2.54

2.26

D++ - Exp

2.97

2.94

3.64

4.51

Besides the marked difference in the MSD/diffusion coefficients between Dowser++ and experimental waters, another noteworthy trend is that, across the board, water molecules in the fully oxidized state exhibit a higher diffusion coefficient than waters in the fully reduced state. This result is well justified if one considers crystal structures of CcO taken in different redox states, but with the same resolution. Only one such pair of structures exists in the PDB (PDB codes 2EIK and 2EIL) and in that case over 50 additional buried water molecules are observed in the reduced state as compared to the oxidized state. A possible physical justification for this result is that, near the oxidized state of the enzyme, water molecules are thought to undergo a reorientation that allows protons to be selectively delivered to the BNC or PLS. The higher degree of dynamicity of the fully oxidized state water molecules, seems to support the notion that 19 ACS Paragon Plus Environment

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

water chains in cytochrome c oxidase are sensitive to the redox state of the enzyme.

Figure 4. Comparison of binding energies between Dowser++ and Experiment. Binding energies were calculated as sum of total of non-bonded potential energies between each water molecule and the surrounding protein environment. The numbers displayed to the right of each colored curve represent the time averaged binding energy values taken over the last nanosecond of simulation. Once again the green D++ - Exp curves represent the set of Dowser++ waters which do not match with experimental molecules within a distance of 2.8Å.

To further elucidate on the trends in the diffusion coefficients and why the unseen water molecules were so dynamic, the binding energies during simulations were calculated (Figure 4). Smaller binding energies correspond to lower retention times for water molecules in their respective energy minima, which correlate to faster diffusion coefficients and more dynamic disorder. As expected the average binding energies were significantly smaller for D++ - Exp 20 ACS Paragon Plus Environment

Page 20 of 37

Page 21 of 37 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

waters than they were for the experimental molecules, with a difference of about 4 kcal/mol in both the oxidized and reduced states. This degree of difference in energy is itself significant, however once again the differences become exacerbated if we zoom-in on the water near the heme groups in the catalytic core of the enzyme. In this region the difference in energy becomes in excess of 12 kcal/mol in both the oxidized and reduced states. This result only further justifies our theory as it illustrates that only the water molecules which are significantly bound by the protein can be seen in crystal structures, while the more loosely bound Dowser++ molecules remain invisible73. It is worth noticing that despite the large difference in energy, the binding energies for Dowser++ are still very feasible for buried water molecules especially in hydrophobic protein regions. Bulk water has a binding energy of -20 kcal/mol, so that the energies in the range of -30 kcal/mol of water molecules predicted by Dowser++ still have a significant energy bias towards the interior of the protein. However, the absolute binding energies must be taken with caution, as they only reflect the forcefield used in MD simulations – in our case TIP3P water model, and CHARM 27 force field for protein and water interactions. Specifically, they do not address the potential polarizability of water molecules in the protein environment. We address these issues in the following section.

3. Dielectric Environment of Water Molecules In total 8 water clusters were identified in Dowser++’s CcO predictions that contained 8 or more water molecules, these clusters are shown in Figure 5. Table 3 illustrates that in clusters where the density/molecular volume of water molecules was an order of magnitude larger than that of bulk water, the Kirkwood G-factor and overall dielectric constant of the water molecules are negligible. However, in regions where the water molecules are packed in density comparable to that of bulk, the G-factor and dielectric constants are large and significant. In particular two

21 ACS Paragon Plus Environment

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

water clusters, one near the top of heme-a3 and the other in the D-channel near Ser 157, have particularly large dielectric constants, especially when compared to what is commonly thought to be the dielectric constant of the dry protein matrix ( ε ≈ 4 − 8 ). Most surprisingly, the D-channel water cluster has a dielectric constant nearly equivalent to that of bulk due to how tightly packed the water molecules in that region are. Formally, the water molecules in that region of the Dchannel (near Ser 157) are packed even a bit more densely than they would be in the bulk; but this seems to be due to a formal algorithm used and the result should be taken as “density comparable to bulk water”. The surrounding protein residues squeezing those water molecules closer together certainly contribute to the effect.

Figure 5.

Snapshot of Dowser++ water clusters identified in CcO taken from fully oxidized state

simulation at timestep 1350. The red spheres represent the center of a particular cluster. The cluster numbers shown correspond to the cluster numbers in Table 3.

22 ACS Paragon Plus Environment

Page 22 of 37

Page 23 of 37 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

Table 3.

Dielectric constant and related parameters of water clusters shown in Figure 5. Quantities

shown are defined in subsection 3 of the Methods section. The error values in the dielectric constant σε w are the standard deviation of the dielectric constant of the water molecules belonging to a particular cluster, while the dielectric constants ε w are the averages. Convergences of G-factors/dielectric constants shown here are presented in the supporting information. Cluster #

Cluster Location

 (Å )



〈  〉 (D2) G-factor

 ± 

1

Heme-A3-Upper

25

52.1

3.4

6.3

10.0 ± 5.9

2

Heme-A3 Lower

11

334.3

4.1

1.2

2.5 ± 0.6

3

Heme-A

10

314.3

3.7

1.1

2.5 ± 0.8

4

D-Channel Near Ser 157

8

18.5

3.8

20.2

30.8 ± 11.6

5

Proton Exit Channel 1

9

227.8

2.9

1.2

2.6 ± 0.9

6

Proton Exit Channel 2

16

105.4

3.3

3.1

5.2 ± 2.2

7

Near His 78

8

299.9

4.2

1.4

2.8 ± 0.8

8

Near Thr 196

8

342.5

4.0

1.1

2.4 ± 0.6

Bulk

-----------------------------

6845

24.0

6.4

26.1

39.7± 0.1

It should be noted these estimates for the dielectric constants are likely to be low, by a factor of roughly 2, due to the fact that calculating the dielectric constant from the TIP3P water model leaves out the electron polarizability of the water molecule67, 74. Including this electronic polarization raises our calculations for the dielectric constant of bulk water to the accepted value of 80. Additionally, and of key importance, to obtain the total dielectric of the protein ( ε ) one needs to add the dielectric of internal water ( ε w ) to the dielectric a dry protein matrix itself ( ε p ), which should be at least as high as purely electronic dielectric constant ε el ~2, but most likely

23 ACS Paragon Plus Environment

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

has a value of between 4 to 848-49.

With such a picture it seems that the notion of a low

dielectric environment of ε ≈ 4 − 8 seems is well supported by our calculations, as half the water clusters analyzed had low dielectric constants of roughly 2.5 (scaled to ~5). Once again, it is only in specific cavities where water in densely packed that the water dielectric constant can be very high, perhaps as great as 60. This result highlights the importance of correctly accounting for all water molecules in a proteins structure when conducting electrostatic calculations. There is currently much interest in developing ways to properly account for and calculate the heterogeneous dielectric environment of the protein75-76; our work suggests that only extremely high resolution structures should be used in order to account for all (or most) interior water molecules, which is of course not always possible. In such a case, Dowser++ appears to provide a useful alternative. Finally it is also worth noting that this method of calculation of the dielectric constant within a protein, while effective and simple, does have its drawbacks. As previously stated it does not take into account the electronic polarizability of water molecules, but that can be treated by introducing a scaling factor that corrects for the electronic polarizability67, 77. The continuum model also assumes each water molecule to be surrounded by a homogenous isotropic dielectric medium, which a protein like cytochrome c oxidase is most certainly not. Clearly, additional studies are needed to elaborate on how best to calculate the dielectric constant in proteins given these limitations. However for the purposes of this study, our results clearly demonstrate the importance of internal water for the dielectric constant of the protein.

Conclusions The predictions of Dowser++ on CcO show many additional water molecules compared with those reported in the best resolution structure available today; this includes some additional 24 ACS Paragon Plus Environment

Page 24 of 37

Page 25 of 37 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

water molecules in the K- and D- channels. These results indicate that the actual hydration of CcO is likely significantly higher than it follows from the current highest resolution crystal structures. However, Dowser++ does not predict significant hydration of the hydrophobic cavity of BNC to form stable water wires that would connect Glu242 to the BNC or PLS, in agreement with previous equilibrium studies. As Dowser++ does not account for protein structural changes, or for a transient hydration, it fails to reproduce the proposed wetting transition21,

32

upon

protonation of propionate D of heme a3. Further hydration of the hydrophobic cavity may indeed require structural changes of BNC and/or rare fluctuations of internal water that will carry proton transfer in the BNC cavity, as shown by several studies in the literature

21, 32, 69

; or entirely

different mechanism is involved. The reason why many of Dowser++ predicted water molecules are not observed in x-ray structures is likely due to the fact that many such molecules are characteristically more dynamic than the crystal structure waters; this could be seen by the difference in the average diffusion coefficient between the Dowser++ molecules and the experimental molecules. In particular, the set of new water molecules predicted by Dowser++ have roughly two to three times the diffusion coefficient of the crystal structure water molecules. Near heme-a3 this difference becomes even more pronounced as the set of new Dowser++ waters have roughly three to four times the diffusion coefficient of the crystal structure waters. Finally, the binding energies of the Dowser++ predictions are also significantly weaker than the binding energies of the crystal structure waters, further supporting the idea that the reason many of Dowser++’s predictions are not seen in experiment is that they are too dynamic and disordered to be resolved. The many unresolved water molecules in CcO produce previously unseen densely packed water cavities. In these densely packed water cavities, water molecules can have very high

25 ACS Paragon Plus Environment

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

dielectric constants, up to ε w

30 (or ~60 when electronic polarizability is taken into account)

for the most densely packed water cluster identified by our calculations. These results highlight the importance of internal water molecules for not only providing the medium and means for proton transfer in the enzyme, but also for the all-important electrostatics of the protein.

Acknowledgements We would like to acknowledge Dr. Alex Morozenko for the development of Dowser++ and for providing the instructions on its use. We would also like to acknowledge Dr. Igor Leontyev for the initial probing the predictions of Dowser++ with MD simulations. This work has been supported in part by NIH GM054052.

Supporting Information SI.docx – Document containing supporting figures that are referenced in the main text.

References 1. Wikstrom, M. K. F., Proton pump coupled to cytochrome c oxidase in mitochondria. Nature 1977, 266 (5599), 271-273. 2. Michel, H., Cytochrome c oxidase: Catalytic cycle and mechanisms of proton pumping-a discussion. Biochemistry 1999, 38 (46), 15129-15140. 3. Konstantinov, A. A., Cytochrome c oxidase: Intermediates of the catalytic cycle and their energycoupled interconversion. FEBS Lett. 2012, 586 (5), 630-639. 4. Zaslavsky, D.; Gennis, R. B., Proton pumping by cytochrome oxidase: Progress, problems and postulates. Biochim. Biophys. Acta, Bioenerg. 2000, 1458 (1), 164-179. 5. Rich, P. R.; Maréchal, A., Functions of the hydrophilic channels in protonmotive cytochrome c oxidase. J. R. Soc., Interface 2013, 10 (86), 20130183. 6. Popović, D. M.; Stuchebrukhov, A. A., Proton exit channels in bovine cytochrome c oxidase. J. Phys. Chem. B 2005, 109 (5), 1999-2006. 7. Sugitani, R.; Stuchebrukhov, A. A., Molecular dynamics simulation of water in cytochrome c oxidase reveals two water exit pathways and the mechanism of transport. Biochim. Biophys. Acta 2009, 1787 (9), 1140-1150. 8. Jünemann, S.; Meunier, B.; Fisher, N.; Rich, P. R., Effects of mutation of the conserved glutamic acid-286 in subunit i of cytochrome c oxidase from rhodobacter sphaeroides. Biochemistry 1999, 38 (16), 5248-5255. 26 ACS Paragon Plus Environment

Page 26 of 37

Page 27 of 37 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

9. Jünemann, S.; Meunier, B.; Gennis, R. B.; Rich, P. R., Effects of mutation of the conserved lysine362 in cytochrome c oxidase from rhodobacter sphaeroides. Biochemistry 1997, 36 (47), 14456-14464. 10. Zhu, J.; Han, H.; Pawate, A.; Gennis, R. B., Decoupling mutations in the d-channel of the aa(3)type cytochrome c oxidase from rhodobacter sphaeroides suggest that a continuous hydrogen-bonded chain of waters is essential for proton pumping. Biochemistry 2010, 49 (21), 4476-4482. 11. Svahn, E.; Faxen, K.; Gennis, R. B.; Brzezinski, P., Proton pumping by an inactive structural variant of cytochrome c oxidase. J. Inorg. Biochem. 2014, 140, 6-11. 12. Lepp, H.; Salomonsson, L.; Zhu, J. P.; Gennis, R. B.; Brzezinski, P., Impaired proton pumping in cytochrome c oxidase upon structural alteration of the d pathway. Biochim. Biophys. Acta, Bioenerg. 2008, 1777 (7-8), 897-903. 13. Brzezinski, P.; Gennis, R. B., Cytochrome c oxidase: Exciting progress and remaining mysteries. J. Bioenerg. Biomembr. 2008, 40 (5), 521-531. 14. Brändén, G.; Gennis, R. B.; Brzezinski, P., Transmembrane proton translocation by cytochrome c oxidase. Biochim. Biophys. Acta, Bioenerg. 2006, 1757 (8), 1052-1063. 15. Wikstrom, M., Cytochrome c oxidase: 25 years of the elusive proton pump. Biochim. Biophys. Acta, Bioenerg. 2004, 1655 (1-3), 241-247. 16. Popovic, D. M.; Stuchebrukhov, A. A., Proton pumping mechanism and catalytic cycle of cytochrome c oxidase: Coulomb pump model with kinetic gating. FEBS Lett. 2004, 566 (1-3), 126-130. 17. Siegbahn, P. E. M.; Blomberg, M. R. A.; Blomberg, M. L., Theoretical study of the energetics of proton pumping and oxygen reduction in cytochrome oxidase. J. Phys. Chem. B 2003, 107 (39), 1094610955. 18. Kim, Y. C.; Hummer, G., Proton-pumping mechanism of cytochrome c oxidase: A kinetic masterequation approach. Biochim. Biophys. Acta 2012, 1817 (4), 526-536. 19. Wikström, M.; Sharma, V.; Kaila, V. R. I.; Hosler, J. P.; Hummer, G., New perspectives on proton pumping in cellular respiration. Chem. Rev. 2015, 115 (5), 2196-2221. 20. Goyal, P.; Lu, J.; Yang, S.; Gunner, M. R.; Cui, Q., Changing hydration level in an internal cavity modulates the proton affinity of a key glutamate in cytochrome c oxidase. Proc. Natl. Acad. Sci. U. S. A. 2013, 110 (47), 18886-18891. 21. Son, C. Y.; Yethiraj, A.; Cui, Q., Cavity hydration dynamics in cytochrome c oxidase and functional implications. Proc. Natl. Acad. Sci. U. S. A. 2017, 114 (42), E8830-E8836. 22. Liang, R.; Swanson, J. M. J.; Wikström, M.; Voth, G. A., Understanding the essential protonpumping kinetic gates and decoupling mutations in cytochrome c oxidase. Proc. Natl. Acad. Sci. U. S. A. 2017, 114 (23), 5924-5929. 23. Blomberg, M. R. A., Mechanism of oxygen reduction in cytochrome c oxidase and the role of the active site tyrosine. Biochemistry 2016, 55 (3), 489-500. 24. Siletsky, S. A.; Konstantinov, A. A., Cytochrome c oxidase: Charge translocation coupled to single-electron partial steps of the catalytic cycle. Biochim. Biophys. Acta, Bioenerg. 2012, 1817 (4), 476488. 25. Belevich, I.; Bloch, D. A.; Belevich, N.; Wikstrom, M.; Verkhovsky, M. I., Exploring the proton pump mechanism of cytochrome c oxidase in real time. Proc. Natl. Acad. Sci. U. S. A. 2007, 104 (8), 26852690. 26. Sugitani, R.; Medvedev, E. S.; Stuchebrukhov, A. A., Theoretical and computational analysis of the membrane potential generated by cytochrome c oxidase upon single electron injection into the enzyme. Biochim. Biophys. Acta, Bioenerg. 2008, 1777 (9), 1129-1139. 27. Riistama, S.; Hummer, G.; Puustinen, A.; Brian Dyer, R.; Woodruff, W. H.; Wikström, M., Bound water in the proton translocation mechanism of the haem-copper oxidases. FEBS Lett. 1997, 414 (2), 275-280.

27 ACS Paragon Plus Environment

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

28. Zheng, X.; Medvedev, D. M.; Swanson, J.; Stuchebrukhov, A. A., Computer simulation of water in cytochrome c oxidase. Biochim. Biophys. Acta, Bioenerg. 2003, 1557 (Supplement C), 99-107. 29. Wikstrom, M.; Verkhovsky, M. I.; Hummer, G., Water-gated mechanism of proton translocation by cytochrome c oxidase. Biochim. Biophys. Acta, Bioenerg. 2003, 1604 (2), 61-65. 30. Tashiro, M.; Stuchebrukhov, A. A., Thermodynamic properties of internal water molecules in the hydrophobic cavity around the catalytic center of cytochrome c oxidase. J. Phys. Chem. B 2005, 109 (2), 1015-1022. 31. Sharma, V.; Enkavi, G.; Vattulainen, I.; Róg, T.; Wikström, M., Proton-coupled electron transfer and the role of water molecules in proton pumping by cytochrome c oxidase. Proc. Natl. Acad. Sci. U. S. A. 2015, 112 (7), 2040-2045. 32. Liang, R.; Swanson, J. M. J.; Peng, Y.; Wikström, M.; Voth, G. A., Multiscale simulations reveal key features of the proton-pumping mechanism in cytochrome c oxidase. Proc. Natl. Acad. Sci. U. S. A. 2016, 113 (27), 7420-7425. 33. Henry, R. M.; Yu, C. H.; Rodinger, T.; Pomes, R., Functional hydration and conformational gating of proton uptake in cytochrome c oxidase. J. Mol. Biol. 2009, 387 (5), 1165-1185. 34. Yamashita, T.; Voth, G. A., Insights into the mechanism of proton transport in cytochrome c oxidase. J. Am. Chem. Soc. 2012, 134 (2), 1147-1152. 35. Popović, D. M.; Stuchebrukhov, A. A., Coupled electron and proton transfer reactions during the o→e transiMon in bovine cytochrome c oxidase. Biochim. Biophys. Acta 2012, 1817 (4), 506-517. 36. Lu, J.; Gunner, M. R., Characterizing the proton loading site in cytochrome c oxidase. Proc. Natl. Acad. Sci. U. S. A. 2014, 111 (34), 12414-12419. 37. Johansson, A.-L.; Chakrabarty, S.; Siöberg, C. B.; Högbom, M.; Warshel, A.; Brzezinski, P., Protontransport mechanisms in cytochrome c oxidase revealed by studies of kinetic isotope effects. Biochim. Biophys. Acta 2011, 1807 (9), 1083-1094. 38. Chakrabarty, S.; Warshel, A., Capturing the energetics of water insertion in biological systems: The water flooding approach. Proteins 2013, 81 (1), 93-106. 39. Supekar, S.; Gamiz-Hernandez, A. P.; Kaila, V. R. I., A protonated water cluster as a transient proton-loading site in cytochrome c oxidase. Angew. Chem., Int. Ed. 2016, 55 (39), 11940-11944. 40. Wikstrom, M.; Ribacka, C.; Molin, M.; Laakkonen, L.; Verkhovsky, M.; Puustinen, A., Gating of proton and water transfer in the respiratory enzyme cytochrome c oxidase. Proc. Natl. Acad. Sci. U. S. A. 2005, 102 (30), 10478-10481. 41. Rasaiah, J. C.; Garde, S.; Hummer, G., Water in nonpolar confinement: From nanotubes to proteins and beyond. Annu. Rev. Phys. Chem. 2008, 59, 713-740. 42. Levitt, M.; Park, B. H., Water: Now you see it, now you don't. Structure 1993, 1 (4), 223-226. 43. Carugo, O.; Bordo, D., How many water molecules can be detected by protein crystallography? Acta Crystallogr., Sect. D 1999, 55 (2), 479-483. 44. Tsukihara, T.; Shimokata, K.; Katayama, Y.; Shimada, H.; Muramoto, K.; Aoyama, H.; Mochizuki, M.; Shinzawa-Itoh, K.; Yamashita, E.; Yao, M., et al., The low-spin heme of cytochrome c oxidase as the driving element of the proton-pumping process. Proc. Natl. Acad. Sci. U. S. A. 2003, 100 (26), 1530415309. 45. Wimley, W. C.; White, S. H., Experimentally determined hydrophobicity scale for proteins at membrane interfaces. Nat. Struct. Mol. Biol. 1996, 3 (10), 842-848. 46. Woo, H. J.; Dinner, A. R.; Roux, B., Grand canonical monte carlo simulations of water in protein environments. J. Chem. Phys. 2004, 121 (13), 6392-6400. 47. Morozenko, A.; Stuchebrukhov, A. A., Dowser++, a new method of hydrating protein structures. Proteins 2016, 84 (10), 1347-1357. 48. Simonson, T.; Brooks, C. L., Charge screening and the dielectric constant of proteins:  Insights from molecular dynamics. J. Am. Chem. Soc. 1996, 118 (35), 8452-8458. 28 ACS Paragon Plus Environment

Page 28 of 37

Page 29 of 37 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

49. Pitera, J. W.; Falta, M.; van Gunsteren, W. F., Dielectric properties of proteins from simulation: The effects of solvent, ligands, ph, and temperature. Biophys. J. 2001, 80 (6), 2546-2555. 50. Lloyd, S., Least squares quantization in pcm. IEEE Trans. Inf. Theory 1982, 28 (2), 129-137. 51. Frohlich, H., General theory of the static dielectric constant. Trans. Faraday Soc. 1948, 44 (0), 238-243. 52. Zhang, L.; Hermans, J., Hydrophilicity of cavities in proteins. Proteins: Struct., Funct., Bioinf. 1996, 24 (4), 433-438. 53. Morozenko, A.; Leontyev, I. V.; Stuchebrukhov, A. A., Dipole moment and binding energy of water in proteins from crystallographic analysis. J. Chem. Theory Comput. 2014, 10 (10), 4618-4623. 54. Ross, G. A.; Morris, G. M.; Biggin, P. C., Rapid and accurate prediction and scoring of water molecules in protein binding sites. PLoS One 2012, 7 (3), e32036. 55. Trott, O.; Olson, A. J., Autodock vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization and multithreading. J. Comput. Chem. 2010, 31 (2), 455-461. 56. Lomize, M. A.; Lomize, A. L.; Pogozheva, I. D.; Mosberg, H. I., Opm: Orientations of proteins in membranes database. Bioinformatics 2006, 22 (5), 623-625. 57. Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K., Scalable molecular dynamics with namd. J. Comput. Chem. 2005, 26 (16), 17811802. 58. Brooks, B. R.; Brooks, C. L.; MacKerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S., et al., Charmm: The biomolecular simulation program. J. Comput. Chem. 2009, 30 (10), 1545-1614. 59. Johansson, M. P.; Kaila, V. R. I.; Laakkonen, L., Charge parameterization of the metal centers in cytochrome c oxidase. J. Comput. Chem. 2008, 29 (5), 753-767. 60. Fee, J. A.; Case, D. A.; Noodleman, L., Toward a chemical mechanism of proton pumping by the b-type cytochrome c oxidases: Application of density functional theory to cytochrome ba3 of thermus thermophilus. J. Am. Chem. Soc. 2008, 130 (45), 15002-15021. 61. Ghosh, N.; Prat-Resina, X.; Gunner, M. R.; Cui, Q., Microscopic pka analysis of glu286 in cytochrome c oxidase (rhodobacter sphaeroides): Toward a calibrated molecular model. Biochemistry 2009, 48 (11), 2468-2485. 62. Samudio, B. M.; Couch, V.; Stuchebrukhov, A. A., Monte carlo simulations of glu-242 in cytochrome c oxidase. J. Phys. Chem. B 2016, 120 (9), 2095-2105. 63. Tuukkanen, A.; Verkhovsky, M. I.; Laakkonen, L.; Wikström, M., The k-pathway revisited: A computational study on cytochrome c oxidase. Biochim. Biophys. Acta, Bioenerg. 2006, 1757 (9), 11171121. 64. Gorbikova, E. A.; Wikström, M.; Verkhovsky, M. I., The protonation state of the cross-linked tyrosine during the catalytic cycle of cytochrome c oxidase. J. Biol. Chem. 2008, 283 (50), 34907-34912. 65. Sharma, V.; Wikström, M., The role of the k-channel and the active-site tyrosine in the catalytic mechanism of cytochrome c oxidase. Biochim. Biophys. Acta, Bioenerg. 2016, 1857 (8), 1111-1115. 66. Leontyev, I. V.; Stuchebrukhov, A. A., Electronic continuum model for molecular dynamics simulations of biological molecules. J. Chem. Theory Comput. 2010, 6 (5), 1498-1508. 67. Leontyev, I.; Stuchebrukhov, A., Accounting for electronic polarization in non-polarizable force fields. Phys. Chem. Chem. Phys. 2011, 13 (7), 2613-2626. 68. Pelleg, D.; Moore, A. In X-means: Extending k-means with efficient estimation of the number of clusters, Proceedings of the Seventeenth International Conference on Machine Learning, Morgan Kaufmann Publishers Inc.: 2000; pp 727-734. 69. Peng, Y.; Swanson, J. M. J.; Kang, S.-g.; Zhou, R.; Voth, G. A., Hydrated excess protons can create their own water wires. J. Phys. Chem. B 2015, 119 (29), 9212-9218.

29 ACS Paragon Plus Environment

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

70. Wikström, M.; Verkhovsky, M. I., Mechanism and energetics of proton translocation by the respiratory heme-copper oxidases. Biochim. Biophys. Acta, Bioenerg. 2007, 1767 (10), 1200-1214. 71. Mackay, J. P.; Landsberg, M. J.; Whitten, A. E.; Bond, C. S., Whaddaya know: A guide to uncertainty and subjectivity in structural biology. Trends Biochem. Sci. 2017, 42 (2), 155-167. 72. Wlodawer, A.; Minor, W.; Dauter, Z.; Jaskolski, M., Protein crystallography for noncrystallographers, or how to get the best (but not more) from published macromolecular structures. FEBS J. 2008, 275 (1), 1-21. 73. Sharpe, M. A.; Ferguson-Miller, S., A chemically explicit model for the mechanism of proton pumping in heme–copper oxidases. J. Bioenerg. Biomembr. 2008, 40 (5), 541-549. 74. Leontyev, I. V.; Stuchebrukhov, A. A., Electronic polarizability and the effective pair potentials of water. J. Chem. Theory Comput. 2010, 6 (10), 3153-3161. 75. Li, L.; Li, C.; Zhang, Z.; Alexov, E., On the dielectric “constant” of proteins: Smooth dielectric function for macromolecular modeling and its implementation in delphi. J. Chem. Theory Comput. 2013, 9 (4), 2126-2136. 76. Sharma, V.; Jambrina, P. G.; Kaukonen, M.; Rosta, E.; Rich, P. R., Insights into functions of the h channel of cytochrome c oxidase from atomistic molecular dynamics simulations. Proc. Natl. Acad. Sci. U. S. A. 2017, 114 (48), E10339. 77. Leontyev, I. V.; Stuchebrukhov, A. A., Polarizable mean-field model of water for biological simulations with amber and charmm force fields. J. Chem. Theory Comput. 2012, 8 (9), 3207-3216.

30 ACS Paragon Plus Environment

Page 30 of 37

Page 31 of 37 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

Table of Contents Image

31 ACS Paragon Plus Environment

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

Figure 1. (A) The structure of cytochrome c oxidase (PDB code: 1V54). The catalytic BNC is located near the center of the high spin heme. The z-axis shown is the same coordinate axis used to construct B and C. (B) The hydrophobicity of the standard protein residues calculated as a function of their depth in the membrane. The Whimley-White hydrophobicity scale36 was used to calculate individual amino-acid hydrophobicities and these hydrophobicities where averaged as a function of their corresponding z-coordinate. (C) The density of buried water molecules in cytochrome c oxidase plotted as a function of the depth in the membrane. 152x87mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 32 of 37

Page 33 of 37 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

Figure 2. (A) Comparison between predictions of Dowser++ and experimental crystal structure water molecules (PDB Code: 1V54) in the D and K channels and near the Heme-a3. This PDB structure corresponds to fully oxidized enzyme. Dowser++ predictions are shown as red and white water molecules while crystal structure waters as shown as blue spheres. (B) Comparison between the z-coordinate water density of Dowser++ predictions and experimental water molecules. The same z-coordinate was chosen as is illustrated in Figure 1A. A detailed analysis of water in D- and K-channels, as well as in the catalytic cavity of BNC, is given in SI. 127x137mm (300 x 300 DPI)

ACS Paragon Plus Environment

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

Figure 3. Comparison of MSD between Dowser++ and Experiment. Comparisons were drawn both amongst all buried water molecules and those only within 10Å of the two heme-groups (i.e. near the catalytic cavity). The green D++ - Exp curves represents the set of Dowser++ waters which do not match with experimental molecules within a distance of 2.8Å. The black dotted lines represent the lines of best fit for each respective colored curve. 152x110mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 34 of 37

Page 35 of 37 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

Figure 4. Comparison of binding energies between Dowser++ and Experiment. Binding energies were calculated as sum of total of non-bonded potential energies between each water molecule and the surrounding protein environment. The numbers displayed to the right of each colored curve represent the time averaged binding energy values taken over the last nanosecond of simulation. Once again the green D++ - Exp curves represent the set of Dowser++ waters which do not match with experimental molecules within a distance of 2.8Å. 152x110mm (300 x 300 DPI)

ACS Paragon Plus Environment

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

Figure 5. Snapshot of Dowser++ water clusters identified in CcO taken from fully oxidized state simulation at timestep 1350. The red spheres represent the center of a particular cluster. The cluster numbers shown correspond to the cluster numbers in table 3. 82x103mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 36 of 37

Page 37 of 37 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

Table of Contents FIgure 44x55mm (300 x 300 DPI)

ACS Paragon Plus Environment