Cytochrome c oxidases (Ccoxs) are the terminal enzymes of the respiratory chain in mitochondria and most bacteria. These enzymes couple dioxygen (O2) reduction to the generation of a transmembrane electrochemical proton gradient. Despite decades of research and the availability of a large amount of structural and biochemical data available for the A-type Ccox family, little is known about the channel (s) used by O2 to travel from the solvent/membrane to the heme a3-CuB binuclear center (BNC). Moreover, the identification of all possible O2 channels as well as the atomic details of O2 diffusion is essential for the understanding of the working mechanisms of the A-type Ccox. In this work, we determined the O2 distribution within Ccox from Rhodobacter sphaeroides, in the fully reduced state, in order to identify and characterize all the putative O2 channels leading towards the BNC. For that, we use an integrated strategy combining atomistic molecular dynamics (MD) simulations (with and without explicit O2 molecules) and implicit ligand sampling (ILS) calculations. Based on the 3D free energy map for O2 inside Ccox, three channels were identified, all starting in the membrane hydrophobic region and connecting the surface of the protein to the BNC. One of these channels corresponds to the pathway inferred from the X-ray data available, whereas the other two are alternative routes for O2 to reach the BNC. Both alternative O2 channels start in the membrane spanning region and terminate close to Y288I. These channels are a combination of multiple transiently interconnected hydrophobic cavities, whose opening and closure is regulated by the thermal fluctuations of the lining residues. Furthermore, our results show that, in this Ccox, the most likely (energetically preferred) routes for O2 to reach the BNC are the alternative channels, rather than the X-ray inferred pathway.
Cytochrome c oxidases (Ccoxs) are the terminal enzymes of the respiratory chain in eukaryotes and in aerobic prokaryotes (reviewed in [1]). These integral membrane proteins belong to the heme-copper oxidases superfamily and couple dioxygen (O2) reduction to the translocation of protons across the membrane. Ccox takes up four electrons from cytochrome c (cyt c) in the positively charged side of the membrane (the inter-membrane space in mitochondria or the periplasm in bacteria) and eight protons from the negatively charged side (eq. 1) [2], [3]: (1) where the subscripts P and N refer to the positive and negative sides of the membrane, respectively. Four of the eight protons reported in equation 1 are used to reduce one O2 molecule and form two water molecules [2], [3], whereas the remaining protons are pumped from the negative to the positive side of the membrane. This overall process contributes to the generation and maintenance of a transmembrane electrochemical proton gradient, which can be further utilized for several energy-requiring processes, such as ATP synthesis [4]. Based on structural and phylogenetic analysis, the heme-copper oxidases superfamily is currently divided into three major subfamilies [5]: A, B and C. The main differences between the three families are the pathways and mechanisms of proton transfer/pumping. The A-type Ccoxs, which are the subject of this work, are widespread through all kingdoms of life [5] and among them are the most thoroughly explored Ccoxs [3], [6], such as the bovine heart mitochondria, the Paracoccus (P.) denitrificans and the Rhodobacter (R.) sphaeroides enzymes. These Ccoxs contain, in the catalytic subunit (subunit I), a low spin heme a and a heterodinuclear center named binuclear center, BNC (Fig. 1A). The BNC is deeply buried in the core of the protein and it is formed by a high-spin heme a3 and a copper ion (CuB). In subunit II, these Ccoxs contain only one redox center, a binuclear copper center named CuA, which accepts electrons from the soluble cyt c and then transfers them to the BNC via heme a. It is believed that protons (both chemical and pumped) are transported from the N-side of the membrane to the BNC via two special proton conducting pathways [3]: the D- and K-pathways (Fig. 1A). A third putative proton-conducting pathway, the H-pathway, was proposed for the mammalian Ccox only [6], [7], and it was suggested to be exclusively used for the transfer of the pumped protons [8]. Several high-resolution crystallographic structures of the A-type family are nowadays available in the literature (e. g. mammalian [7], [9]–[11] and bacterial Ccoxs [12]–[15]) and, based on these structures, it is known that all A-type members share a remarkable structural similarity of the core functional unit formed by subunits I and II (Fig. 1A). Subunit I consists of twelve transmembrane α-helices and contains the BNC and the heme a center. Subunit II is formed by a solvent exposed globular β-sheet domain (which functions as a docking surface for cyt c) and two transmembrane α-helices. It contains only one redox center, the binuclear copper center (CuA). Moreover, at the interface between subunits I and II, Ccox has one Mg+2 ion whose function is still not well understood, but it was suggested to be part of the exit pathway for the pumped protons and for water formed in the BNC [16], [17]. Subunit III, although not considered to be part of the core functional unit, is also highly conserved among the A-type subfamily. Nevertheless, its absence significantly increases the probability of suicide inactivation [18], [19] and thereby reduces the catalytic lifespan of Ccox (in 600-fold or more) [19]. Based in the X-ray data available (eg. [7], [12], [13]), a putative O2 channel for the A-type family was proposed (Fig. 1B). Iwata and co-workers, after pressurizing R. sphaeroides Ccox crystals with xenon, were able to identify a continuous hydrophobic channel that starts in the membrane region of subunit I [13]. This putative O2 channel has two possible entrances that merge together in a region close to the proton-gating residue, E286I (the residues are numbered according to the R. sphaeroides Ccox sequence and the subscript indicates the subunit number). This pathway presents a constriction point which does not allow the access of O2 to the BNC, at least without the occurrence of some conformational change in the protein. Unfortunately, until now, none of the mutagenesis and biochemical studies performed in this channel [20]–[22] was able to clearly demonstrate that it serves as an O2 route into the BNC. All the tested mutations were located too close to the BNC [20], [21], which made the interpretation of the results difficult and did not allow to unambiguously distinguish between the structural obstruction of the O2 channel and the perturbation of the BNC binding kinetics. However, and contrary to the A-type family, in the B-type family the channel used by O2 to reach the BNC is nowadays considered to be well established. The crystallographic studies (with xenon pressurization) performed in the Thermus (T.) thermophilus ba3 enzyme [23]–[25], lead to the identification of a “Y-shaped” hydrophobic channel that runs from the membrane region towards the BNC. This channel, although located roughly at the same position of the putative O2 channel in the A-type Ccox, does not possess a constriction point close to the BNC. In the A-type Ccoxs, the narrowing of the O2 channel is mainly caused by two conserved bulky residues (W172I and F282I in R. sphaeroides [13]), whereas in the B-type Ccox, smaller residues occupy these positions (Y133I and T231I in T. thermophilus [23], [24]). The differences between the A- and B-type regarding the O2 channel are thought to reflect the different functional environments of each type of Ccox. Although the static crystal structures have been a valuable tool for providing insights into the O2 diffusion and for identifying potential O2 channels in Ccox, the elucidation of the molecular basis of O2 diffusion requires the knowledge of the Ccox conformational dynamics. Transiently formed cavities and openings inside the protein (frequently regulated by side chain rotation or by water movements) are not visible in the static X-ray structures, but have already been shown to be very relevant for ligand diffusion (see for example [26]). In this context, molecular dynamics (MD) simulation techniques (with sufficient simulation time and conformational sampling) appear as an alternative for studying the dynamic behavior of proteins and to determine their ligand occupation probabilities inside the protein. In the last decade, computational methods have been widely used to study gas migration in a number of proteins and MD simulations have successfully allowed the identification of several alternative routes for ligand diffusion (e. g. hydrogenase [26]–[28], myoglobins [29], [30], oxidases [31], [32] and laccases [33]). Moreover, the combination of MD simulations with Implicit Ligand Sampling (ILS) [29] calculations allows the calculation of the energy cost of transferring any small, apolar molecule (like O2 or H2) from the solvent to the protein and consequently to compute a 3D free-energy landscape for that specific ligand molecule (e. g [29], [33], [34]). Over the last decades, most of the Ccox research using computational methods focused on the mechanisms and energetics of reduction and/or proton pumping (e. g [17], [35]–[56]). In the A-type Ccox, little emphasis has been given to the identification of the routes used by O2 to move from the solvent towards the BNC, a question only addressed, to our knowledge, by Hofacker and Schulten [31] and by Farantos and co-workers [31]. In the first work, Hofacker and Schulten [31] used MD simulations to study O2 diffusion in the vicinity of the BNC in a bacterial Ccox from P. denitrificans and in the bovine CcOx enzyme. Their simulations revealed a unique, well-defined O2 diffusion channel starting in the membrane-spanning surface of subunit I, close to the interface with subunit III. More recently, Farantos and co-workers [32] have applied the ILS method in order to study the binding of several small gas molecules around the BNC region in the A-type Ccox from P. denitrificans and in the B-type Ccox enzyme from T thermophilus. From these calculations, the authors were able to identify several cavities around the heme a3 region that are conserved in both the A-type and B-type enzymes. This study is however limited to the BNC region, not including other parts of the protein and, consequently, not allowing the analysis of the whole O2 permeation process. The main objective of this work is to identify the O2 channels in the fully reduced Ccox from R. sphaeroides [15] using a combination of MD simulations (with and without explicit O2) and ILS calculations. Our results revealed the existence of three putative O2 diffusion channels. One of channels correlates very well with the channel inferred from the X-ray data available, whereas the other two are alternative routes for O2 to reach the BNC, and were not observed in the X-ray structures pressurized with xenon. Both alternative channels start in the membrane phase and terminate close to Y288I.
Although A-type Ccoxs have been widely studied during the last four decades, the details of the O2 diffusion mechanism are still very incomplete. In particular, the existence and the characteristics of the channel (s) used by O2 to travel from the solvent/membrane to the BNC are still unclear. In this study, we have used an integrated strategy of all-atom MD simulations (with and without explicit O2 molecules) and ILS calculations, designed to examine and characterize the O2 delivery channels in fully reduced Ccox from R. sphaeroides. Altogether, our results suggest that O2 does not diffuse unspecifically inside this protein and instead, uses three well-defined channels running from the interior of the membrane (where O2 solubility is higher than in the aqueous phase) towards the Ccox core. The first pathway has two entrance points, located between helices 5 and 8 and helices 11 and 13 of subunit I, which converges into the constriction point just before the BNC. This channel correlates very well with the channel inferred from the available X-ray structures. The second pathway has only one entry located between the transmembrane helices 13 and 16 of subunit I and it terminates close to Y288I. The third identified pathway approaches the BNC from the subunit II side. This channel runs parallel to the heme a3 hydroxylethylfarnesyl tail and also terminates just below Y288I. According to our observations, the hydrophobic channel detected in the X-ray structures does not constitute the most likely (energetically preferred) entrance point for the O2 molecules in this Ccox. From the O2 affinity map, O2 accesses the BNC via the alternative dynamic channels formed by transient hydrophobic cavities, whose opening and closure is regulated by the thermal fluctuations of the protein. This may be the reason why these channels were not visible in the static X-ray structures. In summary, our results suggest that the original hypothesis (based on static X-ray structures and mutational studies on A-type Ccox) that proposed, that O2 permeation occurs via a unique, continuous and permanently open channel, is indeed a simplification. Our current work does not rule out the role of the X-ray inferred channel, but suggests other alternative routes to the BNC. Furthermore, it emphasizes the need to take into account the dynamic behavior of the protein in order to obtain a more complete description of the O2 putative channels and a more detailed picture of the mechanisms underlying O2 diffusion in these Ccoxs.
The 2. 15 Å resolution crystal structure of the fully reduced Ccox from R. sphaeroides (pdb code: 3FYE) [15] was used as the starting point for this work. This X-ray structure only contains the minimum functional unit (subunits I and II) for Ccox. Only the water molecules with a relative accessibility to the solvent lower than 50% were kept. The relative accessibility of water was computed using the program ASC [66], [67], resulting in the selection of 240 water molecules. Since the GROMOS 54A7 force-field [68] lacks the proper parameterization for the Ccox redox centers, the atomic partial charges for reduced CuA, heme a and BNC centers were calculated using quantum mechanical calculations with the software Gaussian09 [69] and RESP fitting [70], as described in detail in S1 Text in section 1. The van der Waals parameters for the iron atom (located in the two heme groups) were taken from the universal force field [71] whereas the remaining bonded and van der Waals parameters for the metal centers were adapted from the GROMOS 54A7 force field [68]. The protonation state of each individual protonable group at pH 7. 0 was determined using a combination of Poisson-Boltzmann calculations, performed with the package MEAD (version 2. 2. 5) [72]–[74], and Metropolis Monte Carlo simulations, using the program PETIT (version 1. 3) [75]. These calculations were performed using the methodologies described in [75], [76]. For details related with the determination of the protonation state of the protonatable residues, see section 2 in S1 Text. Subunits I and II of Ccox were inserted in a pre-equilibrated dimysristoylphosphatidylcholine (DMPC) lipid membrane (for details related with the membrane construction, equilibration and characterization see [77]). The optimal position of the protein relative to the membrane was determined using the location of the charged residues in the transmembrane helices as a reference. After Ccox insertion into the membrane, all the DMPC molecules located within a cut-off distance of 1. 2 Å from the protein atoms were removed, as described in detail elsewhere [77], [78]. Subsequently, the system (protein, membrane and crystallographic waters) was hydrated in a orthorhombic box using a pre-equilibrated box of SPC water molecules [79]. The water molecules misplaced in the center of the membrane (formed by the highly hydrophobic lipid tails), were removed upon visual inspection. The final system contained the reduced Ccox embedded in a 175 DMPC lipid membrane surrounded by 19,645 water molecules, in a total of 75,178 atoms. All MD simulations were performed using the software package GROMACS 4. 0. 4 [80] together with the united atom GROMOS 54A7 force-field [68] for the protein and lipids and the previously described atomic partial charges and parameters for the redox centers. The simple point charge (SPC) water model was used [79]. Periodic boundary conditions were applied to all simulations. Non-bonded interactions were calculated using a twin range method [81] with short and long-range cut-offs of 8 and 14 Å, respectively. A reaction field correction [82], [83] was applied for the truncated electrostatic interactions, considering a dielectric constant of 62 [84]. The SETTLE algorithm [85] was used to constraint the bond lengths and angle in water molecules, while the LINCS algorithm [86] was used to keep all remaining bonds constrained. The time step for integrating the equations of motion was 0. 002 ps and the neighbor list was updated every 5 steps. The simulations were performed at the constant temperature of 310 K, which is above the phase transition temperature for the DMPC lipids (Tm = 296–297 K) in order to ensure that the membrane is in the liquid crystalline state [87]. A Berendsen heat bath [88] was used, with separate couplings for the protein, membrane and solvent, using a relaxation time constant of 0. 1 ps. The pressure was coupled semi-isotropically (coupling constant of 5. 0 ps and isothermal compressibility of 4. 6×10−5 bar−1 [84]), resulting in an independent coupling of the lateral (Px+y) and perpendicular (Pz) pressures. For all simulations, the x+y and z pressure components were kept at 1 atm and no surface tension was applied [84]. These simulation conditions were shown by Poger et al. [84], [89] to correctly reproduce several experimental measurements for this type of membranes. The system was energy minimized with the steepest-descent method in order to remove excessive strain by performing 5000 steps of steepest-descent minimization with harmonic restraints applied to all non-hydrogen atoms (protein and lipids), followed by further 5000 steps restraining the non-hydrogen atoms of the protein, ending with 5000 steps with restraints applied to the Cα atoms only. After the minimization procedure, and in order to allow proper repacking of the lipids around the protein, a 20 ns MD relaxation was executed in three steps. First, a 0. 5 ns simulation was performed with position restraints to all non-hydrogen atoms of the protein and solvent, at constant temperature and pressure. Afterwards, an additional 0. 5 ns simulation was performed, with position restraints applied to the non-hydrogen atoms of the protein only. Finally, only the Cα atoms were restrained for a period of 19 ns. A force constant of 1000 kJ mol−1nm−2 was used for all the steps that included harmonic position restraints. The unrestrained simulations started after these 20 ns of restrained simulation. In order to reduce the well known sampling problems in membrane-protein simulations, five MD simulations, 100 ns each, were performed, resulting in 0. 5 µs of total simulation time. All replicates were initiated with different sets of random velocities. These simulations will be hereafter designated as O2-free simulations. After 20 ns of restrained simulations, we randomly added 84 molecules of dioxygen (O2) in the solvent zone of each system. No O2 was placed inside the protein nor inside the hydrophobic core of the membrane (see S4 Figure in S1 Text). This new set of simulations will be, hereafter, designated as O2 simulations. The water molecules within a 2 Å distance from the O2 molecules were deleted, similarly to the procedure described in [33]. In order to allow the solvent to adapt to the newly added O2 molecules, a 0. 5 ns MD simulation with position restraints on all non-hydrogen atoms (force constant of 1000 kJ mol−1nm−2) was performed. After this initialization procedure, unrestrained MD simulations were carried out and the simulation conditions and parameters were similar to the ones described previously for the MD simulations without O2, except for the temperature coupling groups used. In this set of simulations, the O2 molecules were included in the same group as the protein. 5 MD simulations, 100 ns each, were performed. The parameters for the O2 molecules were taken from the previously published work of Victor et al [90]. The 84 O2 molecules added to the system corresponds to an O2 concentration of ∼0. 235 M, which is higher than the experimental solubility of this gas in water. However, this high O2 concentration does not affect the structural properties of the protein as shown in S5 Figure in S1 Text. Moreover, the use of this high number of O2 molecules is necessary to obtain reliable statistics within a reasonable simulation time. Sites with high O2 affinity were determined using the ILS method, as previously described in [29]. In this method, the potential of mean force for placing an O2 molecule in any position inside the protein is calculated according to: (2) In equation 2, the implicit ligand potential of mean force, , is an average over a finite number of protein and solvent configurations () and over a number of different equally probable orientations of the ligand (). Moreover, is the Boltzmann constant, is the absolute temperature, and is the interaction energy between the protein and solvent configuration () with the ligand located at position with the orientation. In our case, the O2 free energy map was constructed using the last 50 ns (for each replicate) of the O2-free simulations. For the calculations, all 50005 conformations (= 10001 conformations x 5 replicates in equation 2) were fitted to the X-ray structure using the Cα atoms. A grid of 51×55×87 dimensions was used with a grid spacing of 1 Å and 400 O2 insertions were performed per grid point (= 400 in equation 2). All calculations were carried out using a version of GROMACS 4. 0. 4 Widom TPI algorithm, modified almost in the same way as described in [33]. The only difference is that the ligand insertions here were made within the whole space of the grid cube (the grid cube is centered at the insertion point and with edge length equal to the grid spacing), while in the previous work (described in [33]) the insertions were only possible within the inscribed sphere on the grid cube. The 3D free energy map obtained describes the Gibbs free energy of moving an O2 molecule from vacuum to a given position in the system, ΔGvac→prot (O2). This map was then converted into the ΔGwat→prot (O2) map of interest using a ΔGwat→prot (O2) calculated as described in [33]. The secondary structure assignment was performed with the program DSSP [91]. To determine the percentage of secondary structure loss relative to the X-ray structure, the secondary structure classes considered were: α-helix, 310-helix, 5-helix, β-sheet and β-bridge. For the energy landscape analysis, we used the method described in [33]. In short, this method classifies the energy landscape into energy basins through a steepest-descent tessellation and, afterwards, identifies the lowest-energy point within the boundaries between each pair of neighboring basins, i. e. the saddle point between those basins. After this procedure, a network of paths between all energy minima of the landscape can be constructed using the steepest-descent paths from the saddle points to the minima. A cutoff of 20 kJ·mol−1 was used for the network construction. The errors of the free energy profiles were calculated using two blocks: the first block corresponds to the frames ranging from 50 ns to 75 ns (for all five replicates) whereas the second block contains all the frames ranging from 75 ns to 100 ns. The errors were determined as half the difference of the energies observed between the two blocks for each minina and for each transition. The method used for error calculation assumes that similar minima and pathways can be identified in the two blocks. However, for the O2 channel 1, the pathways connecting M6 to M8 and M10 to M11 were not visible in one of the blocks and by this reason their associated errors could not be calculated.
Cytochrome c oxidases (Ccoxs), the terminal enzymes of the respiratory electron transport chain in eukaryotes and many prokaryotes, are key enzymes in aerobic respiration. These proteins couple the reduction of molecular dioxygen to water with the creation of a transmembrane electrochemical proton gradient. Over the last decades, most of the Ccoxs research focused on the mechanisms and energetics of reduction and/or proton pumping, and little emphasis has been given to the pathways used by dioxygen to reach the binuclear center, where dioxygen reduction takes place. In particular, the existence and the characteristics of the channel (s) used by O2 to travel from the solvent/membrane to the binuclear site are still unclear. In this work, we combine all-atom molecular dynamics simulations and implicit ligand sampling calculations in order to identify and characterize the O2 delivery channels in the Ccox from Rhodobacter sphaeroides. Altogether, our results suggest that, in this Ccox, O2 can diffuse via three well-defined channels that start in membrane region (where O2 solubility is higher than in the water). One of these channels corresponds to the pathway inferred from the X-ray data available, whereas the other two are alternative routes for O2 to reach the binuclear center.
biomacromolecule-ligand interactions
biochemistry
proteins
electron transport chain
biology and life sciences
bioenergetics
biophysics
biophysical simulations
Exploring O2 Diffusion in A-Type Cytochrome c Oxidases: Molecular Dynamics Simulations Uncover Two Alternative Channels towards the Binuclear Site