J. Phys. Chem. B 2005, 109, 1015-1022
1015
Thermodynamic Properties of Internal Water Molecules in the Hydrophobic Cavity around the Catalytic Center of Cytochrome c Oxidase Motomichi Tashiro and Alexei A. Stuchebrukhov* Department of Chemistry, UniVersity of California at DaVis, One Shields AVenue, DaVis, California 95616 ReceiVed: August 19, 2004; In Final Form: October 21, 2004
Cytochrome c oxidase is a redox-driven proton pump that creates a membrane proton gradient responsible for driving ATP synthesis in aerobic cells. The crystal structure of the enzyme has been recently solved; however, the details of the mechanism of its proton pumping remain unknown. The enzyme internal water molecules play a key role in proton translocation through the enzyme. Here, we examine the thermodynamic properties of internal water in a hydrophobic cavity around the catalytic center of the enzyme. The crystal structure does not show any water molecules in this region; it is believed, however, that, since protons are delivered to the catalytic center, where the reduction of molecular oxygen occurs, at least some water molecules must be present there. The goal of the present study was to examine how many water molecules are present in the catalytic center cavity and why these water molecules are not observed in the crystal structure of the enzyme. The behavior of water molecules is discussed in the context of redox-coupled proton translocation in the enzyme.
1. Introduction Cytochrome c oxidase (CcO) has been the subject of many experimental and theoretical studies because of its central role in energy metabolism in aerobic cells and its intriguing protonpumping function.1-7 This transmembrane enzyme reduces molecular oxygen to water and utilizes the free energy of the reduction reaction for the creation of a proton gradient across the membrane; the created proton gradient subsequently drives the ATP synthesis.8 The crystal structure of the enzyme has been recently solved;9-11 however, the mechanism of its proton pumping remains poorly understood.2,12-14 One of the puzzles of CcO concerns its internal water, which is crucial for proton conductivity through the enzyme. In particular, in the cavity surrounding the catalytic center of the enzyme, in which the reduction of oxygen occurs, the crystal structure reveals no water molecules at all, while, in other regions of the enzyme, plenty of internal water has been detected.11 The presence of water in this region is considered to be crucial for two reasons. First, since oxygen reduction in the catalytic center requires the supply of protons, there must be at least some water molecules present in this region to facilitate proton transfer from a well documented proton donor, Glu242,15,16 which is located some 10 Å away from the catalytic site, to the site of oxygen binding. (In this paper, the amino acid numbering of bovine cytochrome oxidase is used.11) Second, recent models of proton pumping17-23 assume that at least some of the pumped protons are conducted between Glu242 and propionate D of heme a3 of the enzyme, across the cavity. Without internal water in this region, such proton translocation would be impossible. Also, the catalytic site is the place in the enzyme where new water molecules are formed, and the fate of these new water molecules is not clear at present.25 Some theoretical models suggest that the new water molecules are not immediately removed from the enzyme but rather occupy metastable states in that cavity and thereby facilitate proton conductivity in this region of the protein; in
particular, the delivery of protons for the catalytic reaction that generates new water molecules can be explained in this way.26 In this paper, using methods of statistical mechanics, we investigate the behavior of water molecules in the catalytic cavity of cytochrome oxidase and address the question of how many, if any at all, water molecules are present there and why they are not seen in the crystal structure. We then discuss the properties of these water molecules in the context of proton conductivity in that region of the protein. The structure of the paper is as follows. In the next section, we describe some relevant details of the enzyme structure and function, review the previous theoretical work, and formulate the problem in greater detail. The computational methods are discussed next, followed by the sections with the results and conclusions. 2. CcO Details and the Formulation of the Problem A detailed review of the structure and function of CcO can be found, for example, in refs 1, 3, and 8; here, we only briefly summarize some relevant results. The CcO oxygen reduction reaction occurs in the binuclear catalytic center (BNC) of the enzyme, which is located in the middle of the protein and consists of an iron ion, Fe(+2/+3), of heme a3 and a copper ion, Cu(+1/+2), of the CuB complex, see Figure 1. During its catalytic cycle, four electrons and four protons from the opposite sides of the membrane are transferred to the BNC, where one oxygen molecule is converted to two water molecules. In addition, for each oxygen molecule reduced, four protons are pumped through the enzyme. Two distinct proton entrance paths leading from the negative side of the membrane to the BNC, the so-called D- and K-channels, are now experimentally well established.27,28 The main proton translocation path is the D-channel, which extends from Asp91 on the negative side of the membrane to Glu242 near the BNC; this channel is believed to conduct up to seven out of eight protons of the catalytic cycle.6,27,28, The crystal
10.1021/jp0462456 CCC: $30.25 © 2005 American Chemical Society Published on Web 12/09/2004
1016 J. Phys. Chem. B, Vol. 109, No. 2, 2005
Figure 1. Structure of CcO around the enzyme catalytic center. Only the important residues are shown. The key elements are heme a, heme a3, the CuB center, Glu242, His291, and the propionates of heme a3. The hydrophobic cavity discussed in the text extends from Glu242 to the ∆-propionate of heme a3 and from Glu242 to the Fe-Cu BNC. The black points show the favorable position of a water molecule in the region.
structure shows the presence of plenty of water in the D-channel, and the formation of a hydrogen-bonded network of water capable of conducting protons along this pathway appears to be possible up to Glu242. The proton path beyond Glu242 however is not clear at all. Extending from Glu242 further up through the protein, toward the oxygen binding site in the BNC and from Glu242 to the propionates of heme a3, is a region that forms the catalytic site cavity; see Figure 1. In this region, the crystal structure shows no water molecules present, and therefore, the proton conduction path seems to break up. Beyond the BNC, however, just above the propionates of heme a3, there is a large hydrophilic cavity where a lot of water is observed, and the hydrogen-bonded network of water appears again, which extends up to the positive side of the membrane. Thus, the protons have to traverse somehow the hydrophobic region around the BNC, for the chemical protons to reach the site where the reduction of molecular oxygen occurs and for the pumped protons to move further up through the protein, to the positive side of the membrane. The absence of water in the catalytic site cavity observed in the crystal structure of the enzyme and the obvious need for such water according to the current understanding of the enzyme function present one of the intriguing puzzles of the oxidase mechanism. Several studies have suggested that in fact some water molecules are present around the BNC, perhaps only transiently, and help the translocation of protons in this region. For example, Zheng et al.26 studied the behavior of water molecules produced by the enzyme and found that these water molecules are not easily removed from the protein and before leaving the catalytic site cavity they form two hydrogen-bonded chains, one connecting Glu242 to the ∆-propionate of heme a3 and the other connecting Glu242 to CuB. The former water chain was suggested to conduct pumped protons to the outer membrane via the ∆-propionate of heme a3, whereas the later provides a path for chemical protons needed in the catalytic site for the reduction of oxygen.17,18 Wikstro¨m et al.19 later extended that model and suggested that the orientation of the dipoles of hydrogen-bonded water molecules in the chain depends on the
Tashiro and Stuchebrukhov redox state of the enzyme, providing an unusual mechanistic coupling between electron and proton transfer in this protein.7 Recently, Olkhova et al.29 published a full atomistic MD simulation that included subunits I and II of cytochrome c oxidase, lipid membrane, and both internal and external water molecules. They investigated the lifetimes and patterns of formation of hydrogen-bonded networks and how water molecules move in the protein. A similar study was published by Cukier.30 They observed significant mobility of individual water molecules along the D- and K-channels and the formation of fluctuating hydrogen-bonded networks. In the simulations of ref 29, no water chain connecting Glu242 to the ∆-propionate of heme a3 was observed for the reduced state; however, such a chain existed in the oxidized state of the enzyme. In the simulation of Riistama et al.,31 the authors predicted a chain of three water molecules connecting Glu242 and one of the histidines (His291) ligated to CuB. In previous theoretical studies, the dynamics of water chains was well investigated only on a relatively short time scale accessible for MD (typically a few nanoseconds, or so), but their thermodynamic stability, a much-longer-time-scale property, was not examined thoroughly. In many cases (see, for example, refs 26, 29, and 32), approximate methods such as Dowser or Grid were used to predict the number and location of internal water molecules. These methods are based on a rough estimate of the interaction energy between water and protein and do not incorporate the effects of water entropy nor even in some cases include the interaction between water molecules properly. It is clear that, for a better description, one needs to calculate accurately the free energy of the system. The number of experimental studies addressing directly the question of water around the BNC is rather limited. Puustinen et al.,33 on the basis of FTIR studies, suggested the existence of a water chain connecting Glu242 and a histidine ligand (His291) coordinated to CuB. Computer simulations of the same group31 were reported to support the proposal and a related CcO proton-pumping model.20 Using EPR, Schmidt et al.25 studied water exit pathways from the binuclear active site, where water is produced, and concluded that water generated in the active site moves out via a distinct pathway, like a D- and K-channel, passing through a region around the magnesium ion located above the BNC. In published X-ray crystal structures9-11 so far, no water molecule is found around the BNC region. Following the logic of the protein function, however, one can conclude that some water molecules must be present around the BNC region. The goal of this paper is to calculate the free energy of different numbers of water molecules in the BNC cavity and to evaluate their thermodynamic properties. First, we examine the behavior of a single probe water molecule around the binuclear center and calculate the potential of mean force (PMF), which represent the energy landscape of the region. The statistical properties of the protein are taken into account by sampling different configurations of the protein from a molecular dynamics (MD) simulation. Using the free energy landscape, we find local minima and estimate energy barriers between them. We then study the behavior of water clusters in the region and calculate their free energies as a function of the cluster size. On the basis of these results, we discuss the behavior of water molecules in the region around the binuclear center of CcO, their thermodynamically favorable numbers, and implications for proton transport in this region of CcO.
Water Molecules in the Catalytic Center Cavity of CcO 3. Theoretical Methods The initial protein configuration was taken from the fully reduced bovine heart cytochrome c oxidase structure (Tsukishima et al.11 PDB code 1V55). As in Zheng et al.,26 subunits I and II were taken from this structure and hydrogen atoms were added using the Amber program.34 The partial charges of the redox centers, heme a, heme a3, CuA, and CuB, were obtained by fitting the electrostatic potential in a Hartree-Fock calculation with the 6-31G* basis set. The calculations were separately performed for each of the redox centers together with their respective ligating residues. For the heme a group, we included His61A and His378A. For the heme a3 group, His376A was included. For the binuclear CuA center, Cys196B, Cys200B, His161B, and His204B were calculated together with two copper atoms, and for the CuB complex, His290A, His291A, and Tyr244A were calculated with a copper atom, for all centers maintaining the geometry of the crystal structure. The truncated carbon atoms were substituted by methyl groups. In these calculations, CuA and heme a were oxidized and heme a3 and CuB were reduced. The protonation states of the titratable residues were determined by the electrostatic continuum method of Popovic and Stuchebrukhov.17 The titratable residues with a proton occupancy larger than ∼0.3 were treated as fully protonated, while all others as unprotonated, so as to make the system electrostatically neutral and to avoid partial proton occupancies which cause ambiguities in the MD simulations. The titration states of the important residues are as follows: all the propiones of hemes a and a3 were deprotonated, Arg438A and 439A were protonated, Asp364A and Glu242A were protonated, His290A and His291A were protonated in the δ-position. (The partial charges of the redox centers are available in the Supporting Information.) After minimization of the hydrogen positions in the system, a molecular dynamics simulation was performed using the Amber program34 to collect protein conformations for the subsequent free energy analysis of the internal water of the protein. In this simulation, we selected all residues which overlap with a sphere with a radius of 10 Å centered at the BNC; most of the atoms of these residues were allowed to move, whereas the rest of the system was fixed in the crystal structure position. Specifically, the sphere was centered at the middle point between the CuB atom and the CMD atom of heme a3. (Here, we use the standard notation for the atoms of the heme groups, and other groups of the protein; thus, CMD is the carbon atom of the methyl group of ring D of the porphyrin; this notation is also utilized in the PDB format.) Such a choice was made in order to include all the important groups around the BNC, such as heme a3, the CuB center, the hydrophobic cavity between Glu242 and the ∆-propionate of heme a3 with its internal water, Glu242, and other residues surrounding the BNC. All atoms of these residues, except for CR atoms, were allowed to move in the simulation. For hemes a and a3, the propionates and hydrogen atoms belonging to the CMD atoms were allowed to move but the other atoms were fixed. The crystal structure internal water molecules were included in the simulation, but only part of them, within 10 Å from the BNC, were allowed to move. In addition, we put nine water molecules around the BNC region, where three water molecules were placed in regions 1 and 2 shown in Figure 1, as they were identified in the previous paper.26 In the simulation, a constant temperature of 300 K was applied. After 200 ps of equilibration, 1 ns of trajectory was accumulated. We used the Sander module of Amber with the SHAKE algorithm, Amber ff99 force field,35 and TIP3P water model.36
J. Phys. Chem. B, Vol. 109, No. 2, 2005 1017 We start the free energy analysis by determining the potential of mean force for a single water molecule in the region near the BNC. Monte Carlo (MC) simulations are used to determine the free energy minima in the region and the free energy barriers between the minima. Different configurations of the protein, obtained in MD, are sampled in the calculations. The free energy minima together with the barriers between them represent the free energy landscape of the catalytic cavity surrounding the BNC. We then investigate the thermodynamic properties of several water molecules that can form different clusters in the region. By comparing the free energy of water in the protein and that in the bulk, we can determine the thermodynamic stability of water inside the cavity. The potential of mean force F for a single water molecule was obtained by integrating the relation dF/dR ) 〈∂H/∂R〉c along the reaction coordinate R, where H is the total Hamiltonian of the water molecules in the cavity, so that the free energy difference between any two points in the protein is
F(B) - F(A) )
∫ABdR 〈∂H ∂R 〉c
(1)
The probe water molecule was allowed to move in a plane perpendicular to the reaction coordinate path, and the ensemble average, 〈〉c, was constrained to that plane. We choose the reaction coordinate to be a distance along the straight line connecting two points in space. The free energy of a water cluster in the protein is calculated by thermodynamic integration (TI)
F(λ)1) - F(λ)0) )
∫01〈
〉
∂H(λ) dλ ∂λ λ
(2)
where 〈〉λ is the ensemble average with scaled Hamiltonian
H(λ) ) Hpw(λ) + Hww(λ) + Hres(λ)
(3)
Here, Hpw(λ) and Hww(λ) are the interaction energies between water and protein and water and water, respectively. They are zero at λ ) 0 and have normal interaction strength at λ ) 1. Thus, F(λ)1) in eq 2 corresponds to the free energy of a water cluster. In the actual calculation, we used the separation-shifted scaling method37 for Coulomb and Lennard-Jones interactions to prevent unstable numerical integration as water CcO
Hpw(λ) )
∑i ∑n
λqiqn
+ (r2in + RC(1 - λ)2)1/2 λBin λAin (4) (r2in + RLJ(1 - λ)2)6 (r2in + RLJ(1 - λ)2)3
In eq 3, we applied the harmonic restraint potential38,39,41,42 for each water in the cluster with the form
Hres(λ) )
1
∑i 2(1 - λ)2Kbr2i
(5)
where Kb is the force constant and ri is the distance of the ith water oxygen from the restraint center which is specified for each water molecule as a putative equilibrium position. Because of this restraint potential, MC sampling near λ ∼ 0 can be limited around the original position. Also, we can analytically evaluate the F(λ)0) term in eq 3, which corresponds to the
1018 J. Phys. Chem. B, Vol. 109, No. 2, 2005
Tashiro and Stuchebrukhov
Figure 2. Largest displacements of the mean atomic positions in the MD simulation from crystal structure 1V55. The displacements are shown against the residue number. The residues without a vertical bar were not allowed to move in the simulation.
free energy of isolated water molecules interacting only with the harmonic restraint potential
F(λ)0) ) -T ln
[
ZNrot
N
( )] [ ( ) Kb
∏∫dr3i exp - 2Tr2i
Λ3NN!
i
-T ln
ZNrot
)
2πT
Λ3NN! Kb
]
3N/2
(6)
where Zrot is the classical rotational partition function of a water molecule, Λ is the thermal de Broglie wavelength, and N is the number of water molecules in the cluster. Notice that we use a rigid-water-molecule model and therefore ignore the effect of intramolecular vibrations in this work. The free energy evaluation method employed in this work is similar to that of previous works of Zhang and Hermans,38 Helms and Wade,39,40 Roux et al.,41 Gilson et al.,46 Olano and Rick,47 and Boresch et al.,42 in the sense that the restraining potentials were employed and their effects on the free energy were estimated analytically. The technical details for our MC calculation are as follows. For parameters of separation-shifted scaling, we used 5.0 Å2 for both the Coulomb and Lennard-Jones potentials. For the harmonic restraint, a force constant of Kb ) 5.0 kcal/mol Å-2 was applied. Thermodynamic integration was performed over 20 discrete values of the parameter λ. For each λ, about 104 steps of equilibration was done and then the configurations were sampled over about 105 steps using the Metropolis method. When the number of water molecules was larger than one, the trial move was applied in each step for a water molecule selected randomly from the water cluster. Statistical errors were estimated to be
