Exploring O2 Diffusion in A-Type Cytochrome c Oxidases: Molecular Dynamics Simulations Uncover Two Alternative Channels towards the Binuclear Site

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.


Introduction
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 (O 2 ) 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]: 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 O 2 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 hemecopper 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 a 3 and a copper ion (Cu B ). In subunit II, these Ccoxs contain only one redox center, a binuclear copper center named Cu A , 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 Hpathway, 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][10][11] and bacterial Ccoxs [12][13][14][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 ahelices and contains the BNC and the heme a center. Subunit II is formed by a solvent exposed globular b-sheet domain (which functions as a docking surface for cyt c) and two transmembrane ahelices. It contains only one redox center, the binuclear copper center (Cu A ). 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 O 2 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 O 2 channel has two possible entrances that merge together in a region close to the proton-gating residue, E286 I (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 O 2 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][21][22] was able to clearly demonstrate that it serves as an O 2 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 O 2 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 O 2 to reach the BNC is nowadays considered to be well established. The crystallographic studies (with xenon pressurization) performed in the Thermus (T.) thermophilus ba 3 enzyme [23][24][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 O 2 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 O 2 channel is mainly caused by two conserved bulky residues (W172 I and F282 I in R. sphaeroides [13]), whereas in the B-type Ccox, smaller residues occupy these positions (Y133 I and T231 I in T. thermophilus [23,24]). The differences between the A-and B-type regarding the O 2 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 O 2 diffusion and for identifying potential O 2 channels in Ccox, the elucidation of the molecular basis of O 2 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][27][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 O 2 or H 2 ) 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,). In the Atype Ccox, little emphasis has been given to the identification of the routes used by O 2 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 O 2 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 O 2 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 a 3 region that are conserved in both the A-type and B-type enzymes. This study is

Author Summary
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 O 2 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 O 2 delivery channels in the Ccox from Rhodobacter sphaeroides. Altogether, our results suggest that, in this Ccox, O 2 can diffuse via three well-defined channels that start in membrane region (where O 2 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 O 2 to reach the binuclear center.
however limited to the BNC region, not including other parts of the protein and, consequently, not allowing the analysis of the whole O 2 permeation process.
The main objective of this work is to identify the O 2 channels in the fully reduced Ccox from R. sphaeroides [15] using a combination of MD simulations (with and without explicit O 2 ) and ILS calculations. Our results revealed the existence of three putative O 2 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 O 2 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 Y288 I .

Results/Discussion
In this work, we investigate the diffusion of O 2 molecules from the solvent to Ccox using MD simulations. We started our study by performing simulations with the enzyme (in the fully reduced state) with explicit O 2 molecules, as described in detail in the Materials and Methods. Although all O 2 molecules were initially placed randomly in the solvent, as time progresses, the gas enters the membrane and concentrates in the lipid tails region (see S8 Figure  in Text S1), similarly to what has been described experimentally [57,58]. After 100 ns of simulation, 46 O 2 molecules were internalized in the membrane, which corresponds to more than 50% of the O 2 placed originally inside the simulation box. Text). In general, before entering Ccox, the O 2 molecules explore the protein's surface and bind briefly to the cavities and niches formed mainly by hydrophobic residues. However, after 100 ns, none of the internalized O 2 molecules was able to reach the BNC in any replicate. Nonetheless, and in order to determine which regions of the protein are more populated by the O 2 molecules during our simulations, we calculated the O 2 probability density maps [59] over the 100 ns (for all five replicates) and the results are depicted in Fig. 2.
As it can be seen, most of the high-affinity regions are located in the membrane phase ( Fig. 2A) and some of them show a good correlation with the putative O 2 channel inferred from the X-ray data pressurized with xenon [13] (Fig. 2B). This channel starts at the membrane spanning region in direct contact with the hydrophobic tails of the lipids and it is formed predominately by uncharged and aromatic residues. It has a Y shape with the two entrances located between helices 5 and 8 and helices 11 and 13 in subunit I. The side chains of I104 I , L105 I , A153 I , L157 I , F108 I , L174 I , V194 I , L246 I , and I250 I contribute to form this channel. The two possible entrances merge together into a constriction point located close to E286 I . This reduction in the diameter of the channel is caused by the phenyl ring of F282 I and the indole ring of W172 I (see Fig. 2B). It has been suggested that during the O 2 Fig. 1. Structure of the membrane bound Ccox from R. sphaeroides [15]. For simplicity, only subunits I and II are represented. A) Architecture of the Ccox core functional unit formed by subunits I (colored in green) and II (colored in cyan). The small blue spheres represent the water molecules resolved in the D-and K-pathways in the X-ray structure [15]. During an O 2 reduction cycle, the electrons are sequentially donated from cyt c and delivered to the BNC, via Cu A and heme a, like in the scheme: cyt c R Cu A R heme a R BNC (heme a 3 /Cu B ) R O 2 . B) Putative O 2 channel inferred from the X-ray data available. The grey mesh represents the putative O 2 channel and it was calculated with the program HOLLOW [ Recent computational studies have suggested that the hydration level of the internal hydrophobic cavity located after the constriction point (formed by F282 I , W172 I and E286 I ) could be important for the regulation of proton transfer in Ccoxs (e.g [55,56]). It was also suggested that the water distribution inside this cavity is modulated by the protonation of the heme a 3 D-propionate and that these hydration changes strongly affect the E286 I proton affinity [56]. This hydrophobic cavity bridges between the end of the D-pathway and the BNC region and it has been predicted to contain several water molecules, at least transiently (e.g. [38,[60][61][62]). However, in our simulations, no water molecules were observed in this region and no.direct and visible connection between E286 I and the BNC was identified.
Since the simulations with explicit O 2 were not able to properly sample the diffusion events (probably due to the simulation timescales), and in order to identify possible alternative O 2 diffusion channels in fully reduced Ccox, the ILS method [29] was used. This method uses the protein conformations obtained from a ligand-free MD trajectory and calculates, for all the positions inside the protein, the potential of mean force for placing a small, apolar, low-interacting molecule at that point. The 3D energy map A) The probability density contours at 0.00015 Å 23 are depicted as a grey surface. The X-ray structure is represented as a ribbon with subunit I colored in green and subunit II in cyan. The hemes are depicted as green sticks. The yellow, light blue and green spheres represent the Cu, Fe (from the heme groups) and the Mg atoms, respectively. The residues forming the putative O 2 channel inferred from the X-ray structure (according to [13]) are colored in orange and the residues forming the constriction point in this channel are highlighted in magenta. B) Zoom image of the putative O 2 channel inferred from the X-ray structure pressurized with xenon [13]. The red arrow identifies the two possible entrance points for this O 2 channel. doi:10.1371/journal.pcbi.1004010.g002 O2 Diffusion in A-Type Cytochrome c Oxidases PLOS Computational Biology | www.ploscompbiol.org generated in this case represents the Gibbs free energy of moving a O 2 molecule from water into any place inside the protein (DG watRprot (O 2 )) and it can be used to infer about the energetically favourable diffusion paths inside a given protein. In these maps, the regions where DG watRprot (O 2 ) is low represent the positions where O 2 has a high probability of residing, meaning that the O 2 affinity in that position is high. The ILS method reduces significantly the sampling problems and it takes into account the dynamic conformational changes of the protein and all the transiently formed cavities, which can be combined to form transitory diffusion pathways.
From the 3D affinity map for Ccox obtained from the ILS calculations (Fig. 3A), we can see that Ccox possesses a complex free energy landscape with several possible O 2 binding cavities, most of them located at the protein surface and that do not directly connect to the BNC. Moreover, the calculated affinity map from ILS not only correlates well with some of the high probability regions found in the explicit O 2 simulations (see Fig. 3B) but also provides a more complete description of the free energy landscape for O 2 inside Ccox due to the sampling of lower affinity zones, not sufficiently sampled by the normal MD simulations. Indeed, the ILS approach allowed us to identified three major routes (named channel 1, channel 2 and channel 3) with high probability of O 2 occupancy, interconnecting the protein surface to the BNC. Two of the channels (channel 1 and channel 2) approach the BNC from the subunit I side whereas the third one (channel 3) approaches the BNC from the subunit II side (Fig. 3A). Interestingly all three channels start in the membrane spanning region where the O 2 concentration is higher than in aqueous phase, which makes physical sense.
Furthermore and in order to identify the energetically preferred routes for O 2 to access the BNC, we extracted, from the ILS affinity map, the lowest free energy pathways connecting the exterior with the O 2 high affinity sites (basins in the O 2 free energy landscape) inside Ccox, using the same methodology as Damas et al [33] (for details, see the Data Analysis section of Materials and Methods). From the analysis of the O 2 free energy landscape (Fig. 3C), we observe that there are many energy minima and low free energy pathways connecting the solvent/membrane region to the interior of the protein. Nonetheless, all these low free energy entrance channels converge into three distinct pathways as we approach the BNC.
The O 2 channel 1 approaches the BNC from the subunit I side and corresponds to the channel inferred from the X-ray data (for R. sphaeroides [13] or for bovine [7]). This pathway has two entry points that are fused together in a free energy minimum located in the constriction point just before the BNC (M 10 in Fig. 4A). The free energy profile for this pathway (Fig. 4A) is characterized by a very high permeation free energy barrier in the constriction point (associated to the bulky side chain of W172 I ) and three deep local minima (M 8 , M 9 and M 10 in Fig. 4A) located between W172 I , F282 I and E286 I for M 10 , L243 I , F282 I and L283 I for M 8 , and I250 I , V194 I and F108 I for M 9 . The local minima observed just below the Ccox surface (M 1 and M 6 ) can probably act as scavengers for the O 2 freely diffusing in the membrane and help to create an O 2 reservoir inside the protein. In this pathway, W172 I and E286 I seem to act as the gateway residues that control O 2 access to the catalytic site. Nevertheless, the passage from the constriction point to the BNC (moving from M 10 to M 11 in Fig. 4A) implies the overcoming of a very high free energy barrier of 39.5 kJ?mol 21 (see Fig. 4A and Fig. 5A), which makes O 2 diffusion via this pathway very slow and difficult.
Over the last 20 years, several biochemical studies [20][21][22] have tried to demonstrate that this pathway is the one used for O 2 diffusion in the A-type Ccoxs, but, until now, no direct measurement proved this unequivocally. All the tested mutations (e.g V287 I I [20,21] and G283 I V [22]) were located too close to the two metals (Fe in heme a 3 and Cu in Cu B ) in the BNC, which made the interpretation of the results very difficult. It was impossible to unambiguously distinguish between the steric hindrance of the O 2 channel (introduced by the mutations) and the perturbation in the structure and local binding kinetics at the BNC. Moreover, until now, only one X-ray structure of the A-type Ccoxs has been pressurized with xenon (considered to be an O 2 analog) [13]. In this structure [13], two hydrophobic xenon binding sites were identified, none of them located beyond the constriction point of channel 1. One xenon was observed at the entrance of the X-ray inferred channel while the other was close to E286 I . Since no xenon was observed after the constriction region, there is still no consensus whether this pathway is indeed a functional channel for O 2 diffusion in these Ccoxs. Additionally, Hofacker and Schulten [31], with the objective of studying O 2 diffusion towards the BNC, performed several MD simulations of Ccoxs with several explicit O 2 molecules. The authors used the locally enhanced sampling (LES) technique [63,64] to identify the pathway used by O 2 both in a bacterial (P. denitrificans) and in the bovine Ccox enzyme. In their simulations, some of the O 2 molecules reached the BNC via a unique and well-defined channel located in the same region as the X-ray inferred channel from [13]. Given the very short duration of their simulations (picoseconds), the observation of such events was only possible due to the lowered free energy barriers experienced by the O 2 molecules when using the LES technique. Nevertheless, a single protein trajectory is still effectively employed in the LES approach, making the protein conformational sampling of that study more limited than the one in the present study. Finally, Farantos et al. [32] by performing ILS calculations in the A-type Ccox from P. denitrificans and in the B-type Ccox enzyme from T thermophilus, were able to calculate the free energy surfaces for the interaction of several small polar (CO, NO) and apolar (O 2 and Xe) around the BNC region. Regarding the X-ray inferred channel, their results correlate very well with the ones reported in the present study and the high affinity cavity Xe 1 observed by the authors in [32] coincides with minima M 11 in the O 2 channel 1.
The channel 2 (Fig. 3A, Fig. 4C and 4D) also starts at the membrane region and has only one possible entrance point located between the transmembrane helices 13 and 16 of subunit I, around residues V234 I , L238 I and V325 I . Moreover, this channel terminates close to Y288 I , a tyrosine residue covalently bonded to H284 I , one of the Cu B -coordinating histidines. Y288 I is located in the end of the K proton-conducting pathway and it has been suggested to be an essential residue for the catalytic mechanism, supplying the proton used in the cleavage of the O-O bond (see [65] for a review). The energy profile for this alternative channel (Fig. 4C and 4D) is characterized by one very low local free energy minimum (M 19 in Fig. 4C) located between Y288 I and T359 I and by a high free energy minimum (M 13 in Fig. 4C) around L292 I and I323 I . The highest energy barriers for this small pathway are 20.4 kJ?mol 21 (Fig. 5B), which are substantially smaller than the energy barrier in the constriction point of the O 2 channel 1 (39.5 kJ?mol 21 ). However, the high energy minimum M 13 located in the middle of the channel (see Fig. 4C) may hinder O 2 diffusion via this pathway.
Lastly, the O 2 channel 3 approaches the BNC from the subunit II side (Fig. 3A, Fig. 4C and 4D) and its entrance is located between the transmembrane helices 28 and 30 of subunit II, around residues F71 II and V167 II . This channel runs parallel to the heme a 3 hydroxylethylfarnesyl tail and also terminates just bellow Y288 I . It is generally formed by hydrophobic and aromatic residues, such as I363 I , F391 I and V272 II (see Fig. 4D). The energy profile for this alternative channel (Fig. 4C) is characterized by only one deep local energy minimum (M 19 in Fig. 4B) located between Y288 I and T359 I and by several smaller (in comparison with the O 2 channel 1) energy barriers, which suggest a higher probability of O 2 passage. In this pathway, the largest energy barriers (19.9-22.3 kJ?mol 21 in Fig. 5C) are related with the zone of the hydroxylethylfarnesyl tail of heme a 3 and the side chains of I363 I and F391 I . This alternative O 2 channel was firstly suggested by Tsukihara et al. [7], based on the X-ray structure of the bovine heart Ccox, as one of the three possible O 2 entrance points in the A-type oxidases. More recently, Farantos et al. [32] were also able to observe the final part of this channel (located close to heme a 3 ) in the A-type Ccox from P. denitrificans.
Even though only subunits I and II were considered in this work, the entry points for all the channels described above are not obstructed by the missing subunits.
Moreover, it is also interesting to notice that all O 2 channels overlap in space with canonical proton pathways. The O 2 channel 1, which corresponds to the X-ray inferred channel, overlies with the end of the D-pathway, whereas the two alternative channels (O 2 channel 2 and channel 3), which end close to Y288 I , overlap with the final part of the K-pathway.
For the two alternative channels identified here, the O 2 pathway is not permanently open and by this reason escaped detection Fig. 3. O 2 free energy landscape obtained from the ILS calculations and lowest free energy pathways obtained from the ILS energy landscape. A-O 2 free energy landscape obtained from the ILS calculations. Two isosurface contours are shown, with free energy levels at -5 (blue inner contour) and 5 (cyan contour) kJ?mol 21 . B-Comparison of the O 2 free energy landscape obtained from the ILS calculations (blue maps) and high-affinity regions from the O 2 explicit MD simulations (red maps). The probability density contours at 0.00015 Å 23 are represented as a red isosurface. C-O 2 lowest free energy pathways obtained from the ILS energy landscape. The spheres represent the local free energy minima while the tubes connecting the minima represent the pathways. The size of the spheres representing the local energy minima scales linearly along the displayed free energy range (small spheres indicates high free energy while large spheres indicate low free energy). The diameter of the pathways scales linearly with the free energy and thus with the oxygen affinity (thinner radius represent high free energies and low oxygen affinity). The free energy at the minima and along the pathways follows the same color scale. The channel 1 coincides with the channel inferred from the X-ray data [7,13] whereas the channel 2 and channel 3 are alternative routes to reach the BNC. doi:10.1371/journal.pcbi.1004010.g003  residues lining the pathways, it is allowed to diffuse further into the protein core until it reaches the BNC.
Deciphering the actual role of these O 2 channels will require further mutational and biochemical analysis. Several mutations sites can be suggested (see Table 1), on the three O 2 channels, in order to clarify the role of each pathway in O 2 diffusion towards the BNC in the A-type Ccoxs.
For the X-ray inferred channel, it is clear that alternative mutation sites, located farther away from the metal ions in the BNC, need to be constructed and studied in order to clarify the actual role of this pathway in O 2 diffusion. W156 I and F282 I are excellent candidates for mutation experiments and their mutation for less bulkier residues, such as tyrosine and threonine (similar to what is observed in the T. thermophilus ba 3 enzyme [23][24][25]) would give important insights into how O 2 diffusion occurs via this channel. Furthermore, for the alternative channels, the residues lining both pathways are also good candidates for mutational studies, because longer and/or bulkier side-chains could, in principle, obstruct the alternative channels and thus clarify the question of whether these channel are, indeed, used to supply O 2 to the BNC in the A-type Ccoxs. Nevertheless, until new experimental data is available, we cannot rule out the hypothesis that all three channels may be working under physiological conditions.

Concluding remarks
Although A-type Ccoxs have been widely studied during the last four decades, the details of the O 2 diffusion mechanism are still very incomplete. In particular, the existence and the characteristics of the channel(s) used by O 2 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 O 2 molecules) and ILS calculations, designed to examine and characterize the O 2 delivery channels in fully reduced Ccox from R. sphaeroides. Altogether, our results suggest that O 2 does not diffuse unspecifically inside this protein and instead, uses three welldefined channels running from the interior of the membrane (where O 2 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 Xray structures. The second pathway has only one entry located between the transmembrane helices 13 and 16 of subunit I and it terminates close to Y288 I . The third identified pathway approaches the BNC from the subunit II side. This channel runs parallel to the heme a 3 hydroxylethylfarnesyl tail and also terminates just below Y288 I . 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 O 2 molecules in this Ccox. From the O 2 affinity map, O 2 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 O 2 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 O 2 putative channels and a more detailed picture of the mechanisms underlying O 2 diffusion in these Ccoxs.

Starting structure
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 Cu A , 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][73][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.

Ccox insertion into a lipid membrane
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 preequilibrated 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.

MD simulations of Ccox inserted in a lipid membrane
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 (T m = 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.6610 25 bar 21 [84]), resulting in an independent coupling of the lateral (P x+y ) and perpendicular (P z ) 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 Ca 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 nonhydrogen 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 nonhydrogen atoms of the protein only. Finally, only the C a atoms were restrained for a period of 19 ns. A force constant of 1000 kJ mol 21 nm 22 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 ms of total simulation time. All replicates were initiated with different sets of random velocities. These simulations will be hereafter designated as O 2 -free simulations.

MD simulations of Ccox inserted in a lipid membrane with explicit O 2
After 20 ns of restrained simulations, we randomly added 84 molecules of dioxygen (O 2 ) in the solvent zone of each system. No O 2 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 O 2 simulations. The water molecules within a 2 Å distance from the O 2 molecules were deleted, similarly to the procedure described in [33].
In order to allow the solvent to adapt to the newly added O 2 molecules, a 0.5 ns MD simulation with position restraints on all non-hydrogen atoms (force constant of 1000 kJ mol 21 nm 22 ) 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 O 2 , except for the temperature coupling groups used. In this set of simulations, the O 2 molecules were included in the same group as the protein. 5 MD simulations, 100 ns each, were performed. The parameters for the O 2 molecules were taken from the previously published work of Victor et al [90].
The 84 O 2 molecules added to the system corresponds to an O 2 concentration of ,0.235 M, which is higher than the experimental solubility of this gas in water. However, this high O 2 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 O 2 molecules is necessary to obtain reliable statistics within a reasonable simulation time.

Implicit ligand sampling (ILS) calculations in Ccox
Sites with high O 2 affinity were determined using the ILS method, as previously described in [29]. In this method, the potential of mean force for placing an O 2 molecule in any position inside the protein is calculated according to: In equation 2, the implicit ligand potential of mean force, W (r), is an average over a finite number of protein and solvent configurations (M) and over a number of different equally probable orientations of the ligand (C). Moreover, k B is the Boltzmann constant, T is the absolute temperature, b~(k B T) {1 and DE(r,q m ,V k ) is the interaction energy between the protein and solvent configuration (q m ) with the ligand located at position r with the orientation V k .
In our case, the O 2 free energy map was constructed using the last 50 ns (for each replicate) of the O 2 -free simulations. For the calculations, all 50005 conformations (M = 10001 conformations x 5 replicates in equation 2) were fitted to the X-ray structure using the C a atoms. A grid of 51655687 dimensions was used with a grid spacing of 1 Å and 400 O 2 insertions were performed per grid point (C = 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 O 2 molecule from vacuum to a given position in the system, DG vacRprot (O 2 ). This map was then converted into the DG watRprot (O 2 ) map of interest using a DG watRprot (O 2 ) calculated as described in [33].

Data analysis
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: a-helix, 3 10 -helix, 5-helix, bsheet and b-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 21 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 O 2 channel 1, the pathways connecting M 6 to M 8 and M 10 to M 11 were not visible in one of the blocks and by this reason their associated errors could not be calculated.