Heme Cavity Dynamics of Photodissociated CO from ba3-Cytochrome

Aug 7, 2009 - Massimiliano Porrini,*,† Vangelis Daskalakis,† Stavros C. Farantos,†,‡ and. Constantinos Varotsis‡. Institute of Electronic St...
0 downloads 0 Views 3MB Size
Download PDF
J. Phys. Chem. B 2009, 113, 12129–12135

12129

Heme Cavity Dynamics of Photodissociated CO from ba3-Cytochrome c Oxidase: The Role of Ring-D Propionate Massimiliano Porrini,*,† Vangelis Daskalakis,† Stavros C. Farantos,†,‡ and Constantinos Varotsis‡ Institute of Electronic Structure and Laser, Foundation for Research and Technology-Hellas (FORTH), P.O. Box 1527, Vasilika Vouton, Heraklion 71110, Crete, Greece, and Department of Chemistry, UniVersity of Crete, P.O. Box 2208, Vasilika Vouton, Heraklion 71305, Crete, Greece ReceiVed: May 13, 2009; ReVised Manuscript ReceiVed: July 1, 2009

Intracavity molecular dynamics studies of photodissociated carbon monoxide from ba3-cytochrome c oxidase have been performed by sampling the phase space with several hundreds of trajectories each integrated up to 100 ps time interval. It is shown that the cis conformation of protonated ring-D propionate of heme a3 and its trans conformation for the deprotonated species control the CO location by creating two distinct equilibrium states for CO confined in a cavity internal to the distal heme pocket. Thus, these cis (closed gate) and trans (open gate) conformations of heme a3 propionate D play the role of a switch, opening or closing a gate for confining CO in a cavity internal to the heme pocket or releasing it to a bigger outer cavity. The geometry of the inner cavity and the validity of the potential function employed are further investigated by Density Functional Theory calculations for the active site, potential of mean force curves along the copper-CO bond, as well as with Quantum Mechanics/Molecular Mechanics calculations. In the light of the present study, trajectory scenarios for the dissociation of CO previously suggested from time-resolved infrared spectroscopy are re-examined. Introduction Cytochrome c oxidase (CcO) is the last enzyme in the respiratory electron transport chain of mitochondria and bacteria. Receiving four electrons from cytochrome c and binding four protons from the mitochondrial matrix, it operates the reduction of one molecule of oxygen to two of water1 and, in the meanwhile, it pumps four protons across the inner mitochondrial membrane so that it helps keeping the transmembrane gradient of the proton electrochemical potential2,3 necessary for synthesizing ATP. The enzymatic complex ba3-CcO from Thermus thermophilus contains two hemes, b and a3, and two copper centers, the bimetallic CuA and CuB. In fact, the iron of heme a3 and CuB form a binuclear center that is the site of oxygen reduction;4 a scheme of this enzyme and a magnification of the active site region are shown in Figure 1. As it has already been demonstrated in the literature,5-11 to follow the dynamics of whatever ligand (like CO, O2, and NO), it is very important to consider the channels that the ligand can go through, the cavities where the ligand can spend significant portions of time, and the diffusion of the ligand itself. Carbon monoxide has a high oscillator strength in the infrared region and its frequency and bandwidth depend a lot on the environment,12 so that it can be used as a probe for studying ligation reactions. Many experiments regarding the photodissociation of CO have been performed in the last decades and they have mainly been addressed to the understanding of the binding reactions in the binuclear center. For example, the studies on the intermediates of CcO + O2 reaction13,14 have gotten large * Corresponding author. Telephone: +30 2810 545058. Fax: +30 2810 1305. E-mail: [email protected]. † Institute of Electronic Structure and Laser, Foundation for Research and Technology-Hellas (FORTH). ‡ University of Crete.

advantages from flow-flash experiments of CO photodissociation.15,16 Furthermore, the CO reactions give us information about the possible paths followed by the ligand to reach and to leave the active site and about the cavities where the ligand can be trapped. Besides, we have to consider that the dynamic and kinetic processes in protein-ligand systems cover a large range of time scales: from the ultrafast and fast dynamics, which span a range of a few hundreds of femtoseconds to picoseconds, up to the slow dynamics, from microseconds to milliseconds. A lot of experimental and theoretical studies have been conducted in the last decades and they were aimed at rationalizing both the dynamics and kinetics of photodissociated small ligands and the roles played by the protein internal structure. For CO ligand it was experimentally proved that it binds CuB in less than 1 ps after dissociation from heme a3, whereas a very slow kinetic process takes place in milliseconds time scale, where CO comes back to its initial position; this process is the so-called geminate recombination. Vos and co-workers17-20 accomplished studies about the fast and ultrafast dynamics of CO. They observed that, in the binuclear center of several cytochrome c oxidases, the CO ligand transfer process from heme a3 to CuB atom lasts around 300 fs and it occurs predominantly via a ballistic mechanism rather than a diffusive one. Koutsoupakis et al.21-23 analyzed experimentally the microsecond time scale kinetics, and apart from the transfer from Fe to Cu, they identified the presence of a docking site, internal to heme pocket, which shelters around 20% of CO ligand (the other 80% is bound to copper). The investigators applied timeresolved step-scan Fourier transform infrared (TRSS-FTIR) difference spectroscopy, a technique capable of probing from nanosecond to millisecond time scales. Molecular dynamics (MD) calculations were performed by Vos and co-workers19 of two picoseconds duration and employ-

10.1021/jp904466n CCC: $40.75  2009 American Chemical Society Published on Web 08/07/2009

12130

J. Phys. Chem. B, Vol. 113, No. 35, 2009

Porrini et al.

Figure 1. Main components of ba3-cytochrome c oxidase from Thermus thermophilus with a magnification of the active site on the right side. The two hemes b and a3 (in the latter ring-A and ring-D propionates are shown) and CuB atom of the binuclear center coordinated by three histidines and cross-linked Tyr237 are depicted. The green helices form chain-A of the protein, the orange helix is chain-B, and the blue helix is chain-C: Green ball is iron, purple is copper, cyan balls and sticks are carbons, red balls and sticks are oxygens, and blue sticks are nitrogens. The figure was created using VMD software.65

ing a model of the Paracoccus denitrificans CcO. The main purpose of these simulations was to illustrate the binding of CO to CuB. To the best of our knowledge, longer time dynamics of CO trapped in the cavity internal to the distal heme a3 pocket of CcO have not been studied yet. On the other hand, there are numerous MD investigations of ligand migration in myoglobin up to nanosecond time scales.24 The purpose of the present work is to investigate the dynamics of CO in the cavity inside the distal heme pocket after the photodissociation of CO from iron of ba3-cytochrome c oxidase of Thermus thermophilus; hereafter, we shall call this cavity “inner cavity” to distinguish it from an “outer cavity” (external and adjacent to heme pocket and defined more accurately later on). Because we are not able to cover by MD the experimental microseconds time scales we study the dynamics of well specified conformations of the active site such as protonated or deprotonated chemical species, inducing cis to trans conformational changes, the formation of which could be triggered in microseconds time scale. In particular, the role played by the cis/trans isomers of the ring-D propionate (due to protonation events) in trapping CO is revealed. To achieve that, we construct potential models optimized by Density Functional Theory (DFT) electronic structure calculations for the binuclear center and the cavity above heme a3. The electronic excitation of heme a3CO from the bound ground state to the unbound excited one is modeled in the same way as in previous studies by the “sudden” approximation.25-27 We assume that the experimental results of Koutsoupakis et al.21,22 are referred to thermodynamic equilibrium states. Thus, by making the appropriate potential model pertinent to each equilibrium state, such as CO bound to iron, copper, or inside the cavity, we run trajectories up to 100 ps long to investigate isomerization reactions. In this way, we have identified the inner and outer cavities, the docking site geometry and the orientations of CO in it. The rest of the article is divided into three sections: the first one deals with the definition of the system, computational techniques, and the construction of the potential functions. In the second one, the results are presented and discussed with respect to experimental observations. Finally, in the third section, we summarize and draw the main conclusions of this study.

Potential Functions and Computational Methods The calculations were based on the X-ray structure of ba3CcO from Thermus thermophilus (protein data bank coded 1XME28). It consists of three main chains (named A, B, and C, see Figure 1), two hemes (b and a3), three atoms of Cu (two CuA and one CuB), and several crystallization water molecules. This structure was saturated with hydrogen atoms by employing the Pymol program29 and a molecule of CO was bound to Fe of the heme a3. The total system contains 12336 atoms. The initial crystallographic configuration had the ring-D propionate of heme a3 in a trans conformation and hydrogen bonded to a glycerol molecule, added to stabilize the crystal structure of the protein.30 We removed this molecule and minimized the energy of the protein. If propionate D is protonated, it forms a hydrogen bond to the nearby glycine (Gly232) residue in a cis conformation. Deprotonation of heme a3 ring-D propionate results in a trans optimized geometry for propionate D itself. In Figure 1, the active site shows the trans conformation of the propionate D. Ring-A propionate of heme a3 is hydrogen bonded to an aspartic acid (Asp372) residue and it conserves the trans configuration. These findings have led us to investigate the dynamics of ba3-CO system for different conformations of propionate D: cis and trans. We used the Molecular Dynamics package of programs Tinker31 to construct the potential energy surface (PES). We employed the parameters from Amber99 force field32 and the TIP3P model33 for water molecules. The active site made up of 94 atoms consists of the heme a3 with the residue His384 at proximal and CO at distal coordination scheme, and CuB ligated to the residues His282, His283, and His233, while the latter is cross-linked to the Tyr237 residue. Parameters for the metal atoms, Fe and Cu, and the CO molecule or the heme groups were taken from Charmm27 force field34 and further optimized by reproducing geometrical and spectroscopic results from DFT calculations of the active site. Specifically, geometric parameters for the active site, with CO either bound to Fe(II) of heme a3 or to CuB(I), were based on DFT studies of CcO models containing the binuclear center.35 Vibrational frequencies for CuB-C-O moiety, such as the vibrational stretching frequencies (ν) of CO at 2031 cm-1, the frequency of CuB-C bond at 414

Heme Cavity Dynamics of Photodissociated CO cm-1, and the bending frequency ν(CuB-C-O) at 326 cm-1 were obtained from DFT calculations using the DMOL3 program36 with the BLYP/DND functional. To obtain vibrational frequencies for heme Fe-C-O moiety, we designed a model molecule in which the axial coordination of ferrous iron contains an imidazole and a CO ligand. Heme a3 is taken without the propionate groups, and it is represented by an unsubstituted iron-porphyrin. A full geometry optimization (restricted singlet state) was performed by using the density functional RI-BLYP method. No restrictions were applied on the geometry during optimization. We used the triple-ζ valence basis set (tzvp), as implemented in Turbomole software package.37 Subsequently, evaluation of the Hessian matrix at the optimized geometry and a computation of the harmonic vibrational spectrum gave the stretching frequencies of ν(CO) at 1993 cm-1, ν(Fe-C) at 440 cm-1 and the bending frequency ν(Fe-C-O) at 542 cm-1. The different DFT methodologies used give results in good agreement to the experimental values, when available, justifying the employment of diverse level of theory for different parts of the system. The interaction potentials of the carbon atom of CO with Fe and Cu in the ground electronic state of the protein are described by Morse type functions, so is the bond of CO itself. Whereas an exponential function of Fe-CO distance was implemented to describe the dissociation of CO. All the values of the parameters and related descriptions are given, respectively, in Table 1 and in the first section of Supporting Information (SI). Although, the fitting of the potential parameters of the active site was achieved by employing Tinker software, the calculations for the total protein were performed with the Molecular Dynamics package of programs DL_POLY.38 Care was taken to transfer the force field (FF) from Tinker to DL_POLY format. We properly modified DL_POLY program to ensure exactly the same potential function in both programs. Most of the calculations were carried out in the EGEE European Grid computational infrastructure39 through HellasGrid,40 which is part of EGEE. We submitted batches of 500 trajectories each to the Grid as members of the virtual organizations CompChem41 and SEE (South East Europe Grid).42,43 After thermalizing the system with a 20 ps run trajectory at 300 K, a simulation for 1 ns at constant energy was run to collect initial coordinates and velocities of the atoms in the protein. These data together with the excitation energy of Fe-CO bond determine the initial conditions of the trajectories run. Results and Discussion Vibrational Frequencies. The optimized potential energy surface for the ground electronic state was initially tested by calculating characteristic vibrational frequencies. These are obtained by Fast Fourier Transforms (FFT) of time series of Fe-C-O and Cu-C-O species bonds and angles. The results are depicted in Figure 2a for the active site alone and in Figure 2b for the whole protein (details of the computation are provided in the caption of the figure). Red lines are used for the spectra of Cu-C-O and black lines for the Fe-C-O moiety. In the active site and for the Fe-C-O moiety, we find the stretching frequency of Fe-CO bond to be equal to 489 cm-1 and the bending frequency equal to 568 cm-1. These frequencies can be compared with the experimental values of 507 cm-1 for the stretching and 568 cm-1 for the bending modes given in ref 44. For Cu-C-O species the peaks at 419 cm-1 and 361 cm-1 are assigned to the stretching and bending vibrational modes, respectively, whereas the corresponding frequencies from the DFT calculations are 414 and 326 cm-1. Generally, comparing

J. Phys. Chem. B, Vol. 113, No. 35, 2009 12131

Figure 2. Power spectra obtained from Fe-CO (black line) and Cu-CO (red line) bond time series FFT. (a) Active site: time series were collected every 3 fs along a simulation time of 225 ps. (b) Whole enzyme: time series were collected every 2 fs along a simulation time of 1 ns.

the spectral peaks of the active site with those of the total protein we observe a downshift in all lines. The vibrational frequencies of C-O molecule bound to iron/ copper are 1993/2031 cm-1, respectively, and they can be compared with the experimental values44 of 1973 cm-1 and 2053 cm-1, respectively. Over all, the frequencies are found to be in reasonable agreement to experimental ones. Dissociation Dynamics. According to the “sudden approximation” model of the dissociation process, the Morse type interaction of iron atom of heme a3 with the carbon atom of CO, valid in the ground electronic state, is replaced by a repulsive force described with an exponential function, whereas at the same time the intermolecular interactions of carbon with the environment of iron are switched off. The nonbonded interactions of carbon are turned on again to that of the ground electronic state after CO has dissociated. Two types of dissociating trajectories are found: those in which carbon monoxide is attracted by CuB in less than 1 ps, it loses its energy and it remains trapped in the Cu-CO Morse potential well, and those in which CO either escapes to the heme cavity immediately after dissociation or it first approaches copper but it is not trapped by its attractive potential. In the latter type of trajectories the nonbonded interactions of CO with the environment of CuB are turned on and the dynamics of CO is followed inside the cavity for 30 ps and in several cases up to the maximum of 100 ps. Representative trajectories from each type are shown in Figure 3 by plotting the distance of Cu-CO for the first 3 ps (left) and longer times up to 30 ps (right). From the left panel it is clearly visible how CO is trapped either from copper (red line) in less than 1 ps or sheltered into the cavity (black line) after bouncing a few times on the residues of the heme pocket (see Movies 1 and 2 provided in SI). Trajectories trapped in the cavity have larger oscillation amplitudes than trajectories trapped in CuB well. This means that it takes longer for CO to dissipate its energy when it dissociates into the cavity. Nevertheless, we have found that for all trajectories CO reduces its total energy (internal plus translational) to about 2.5 kcal/mol in less than 2 ps (see Figure 1 in SI). Integration of the trajectories for longer times and up to 100 ps reveals that CO may return to copper visiting it for short time intervals. It is beyond our present computer capabilities to simulate ba3-CcO at micro and milliseconds time scales studied spectroscopically and keeping at the same time the accuracy of the employed PES. Hence, here we focus on

12132

J. Phys. Chem. B, Vol. 113, No. 35, 2009

Porrini et al.

Figure 3. CuB-CO bond time series for two typical dissociating trajectories of the ba3-CcO enzyme with propionate D in cis conformation. The black line is for a trajectory trapped in the cavity internal to distal heme a3 pocket; the red line is a trajectory trapped in CuB minimum.

investigating the picoseconds dynamics of CO for specific geometries of the cavity. As it was discussed in the previous section there are different minimum energy structures when ring-D propionate is protonated, a conformation which we call closed, and deprotonated, a conformation called open. These are due to the hydrogen bonding of propionate D. Rotation around the dihedral angle of propionate D, δ(CA-CT-CT-CA), by 180° produces the cis and trans conformations. The symbols are defined as follows: the first CA is the carbon atom with sp2 hybridization that propionate shares with heme a3, the two CTs are the sp3 carbon atoms (around whose internuclear axis the rotation takes place) and the last CA is the sp2 carbon atom of the carboxylic group (see also Figure 2 in SI, where more DFT model calculations are presented). Batches of 500 trajectories were run and the location of CO inside the heme cavity was identified. In Figure 4a we present the position of the center of mass of COs as orange spheres after evolution of 30 ps for the protonated propionate D. The closed geometry of propionate and its hydrogen bond with Gly232 can be seen. In Figure 4b the projections of CO center of mass on the heme plane (defined by the iron atom and two nitrogen atoms) are shown. Similarly, for the open geometry of the deprotonated propionate D the distribution of COs obtained from 500 trajectories of 30 ps long are depicted in Figure 5. Apart from the obvious confinement of CO distribution in the closed structure of the cavity and the wide one in the open conformation, distinct differences are found when the projections of CO center of mass on the heme plane are examined. For the closed structure, we can see that CO is accommodated between rings A and D of heme a3 and close to the two metal atoms of the binuclear center. Contrary to that, in the open geometry more trajectories exhibit CO escaping from the heme pocket to a larger cavity external to the heme pocket and the majority of COs are concentrated above the ring-D of heme a3 (see Movie 3 provided in SI). Batches of trajectories have been executed at much higher excitation energies and with protonated propionates. Details and results are given in the SI (Figures 3 and 4). The outer cavity represents the lower stretch of the characteristic Y-shaped hydrophobic channel of CcO from Thermus thermophilus, made up by seven Xe-binding sites and proposed in ref 45. This Y-shaped channel is supposed to permit oxygen to get into the active site and water to escape outside of the CcO.46,47 Nevertheless, the results of the deprotonated propionate lead one to infer that ring-D propionate acts as a gate for

Figure 4. (a) Distribution of dissociated CO molecules inside heme cavity for the cis (closed gate) ring-D propionate. Orange balls depict the center of mass of CO molecules after 30 ps time evolution and for 500 trajectories. (b) Projections of CO center of mass on a plane defined by iron and two nitrogen atoms of heme a3 and for the closed configuration of ring-D propionate.

controlling the ligand (and likely the water molecules resultant from the oxygen reduction) to escape from the active site. A similar function of ring-D propionate has also been proposed in ref 48, where the propionate D of CcO from bovine heart mitochondria forms a gate with arginine through a hydrogen bond. More quantitatively, the distributions (in bins) of 500 trajectories of the mean value of the Fe-CO and Cu-CO distances averaged over each trajectory of 30 ps run are depicted in Figure 6 for both cases: (a) cis and (b) trans ring-D propionate. From these plots one can infer that, in the closed gate configuration, COs are clustered in an area located around 3.3 Å from both Fe and Cu atoms and with carbon atoms being in majority slightly closer to iron than to copper, whereas with the open gate, COs are quite spread in the inner and outer cavity and equidistant from both metals, at about 5.3 Å. In this figure it is fairly evident that the propionate can control the dynamics of the ligand (in general and not CO in specific) and, besides, its conformation would depend on the pH value. It has been argued that pH and the proton pumping mechanism are interrelated in controlling the movement of a ligand in and out the active site.49-53 The above remarks are in accord to proposed mechanisms for proton pumping in CcOs, which involve H+ donors and acceptors in an area proximal to both ring-D propionates of the hemes b and a3, see refs 54-57. Koutsoupakis et al.23 have studied the implications of protonation of ring-A propionate in association to the Q-proton pathway in the ba3CcO enzyme. As far as the orientation of COs with respect to iron and copper atoms is concerned, we plot the distributions of the mean

Heme Cavity Dynamics of Photodissociated CO

J. Phys. Chem. B, Vol. 113, No. 35, 2009 12133

Figure 7. (a) Histograms of the mean values over 30 ps time evolution of the Fe-C-O and Cu-C-O angles: (a) cis ring-D propionate, (b) trans ring-D propionate. The x-axis represents bins of 10° of width and starting from 20°.

Figure 5. (a) Distribution of dissociated CO molecules inside the heme cavity for the trans (open gate) ring-D propionate. Orange balls depict the center of mass of CO molecules after 30 ps time evolution and for 500 trajectories. (b) Projections of CO center of mass on a plane defined by iron and two nitrogen atoms of heme a3 and for the open configuration of ring-D propionate. The black line denotes an approximate border of inner and outer cavities.

Figure 6. Histograms of the mean values over 30 ps time evolution of the Fe-CO and Cu-CO distances: (a) cis conformation of propionate D, (b) trans conformation of propionate D. The x-axis represents bins of 1 Å of width and starting from 1.811 Å.

values of Fe-C-O and Cu-C-O angles over trajectories of 30 ps long in Figure 7a, for the cis structure of propionate D, and Figure 7b, for the trans isomer. As one can see, the Fe-C-O angle has an equilibrium value around 75° in both cis and trans cases, but with a broader distribution for the closed structure, which means that with propionate open, CO can be docked in the cavity (“above” the ring-D) more firmly with respect to iron. Conversely, for the Cu-C-O angle distribution we have two different equilibrium values for the cis/trans cases: for the closed propionate the peak of the angles band is located around 100° with a pretty broad bandwidth, whereas in the open propionate case the mean value drops to 85° and the distribution

is again narrower, confirming that with the opening of propionate the orientation of CO in the cavity is conserved. In spectroscopic studies of myoglobin58-60 B states are assigned to CO just after dissociation to be distinguished from the A states of CO bound to iron. In more detail, states A0 and B0, A1 and B1, and so on have been assigned according to spectroscopic and kinetic studies. Koutsoupakis et al.22 ascribe the CO sheltered in the cavity to a B1 state with a vibrational frequency (2131 cm-1) close to the free gas carbon monoxide (2143 cm-1), and they argue that CO is oriented with oxygen toward iron and close to ring-A propionate. In our calculations we found that for both closed and open configurations about 62.5% of CO molecules have the oxygen atom oriented toward Fe, and in the case of the open structure, the expected position of CO is “above” ring-D. Our results are similar to those found for myoglobin,58 where the B1 state is with the oxygen toward iron atom of heme a3. However, we point out that in ref 61, the authors claim that this state corresponds to the CO with carbon oriented toward iron. To further test the empirical potential function and the geometric characteristics of the cavity, we have performed QM/ MM calculations for the open gate conformation of propionate D by employing the Schro¨dinger suite of programs.62 Details for geometry optimization and frequency calculations are given in Figure 5 of Supporting Information. The results demonstrate that CO does interact with its environment producing different vibrational frequencies as it was found in the experimental spectra,21 where the B1 and B0 states were attributed to different CO orientations or protein matrix conformations. Koutsoupakis et al.21 have observed with TRSS-FTIR difference spectroscopy two states of CO attributed to the B1 and B0 found in Mb. They are revealed with two peaks at 2131 and 2146 cm-1 in the spectra of ba3-cytochrome c oxidase at room temperature and pH equal to 5.25. At 10 and 35 µs measurements subsequent to CO photolysis only the 2131 cm-1 mode is present and only after 50 µs the 2146 cm-1 peak becomes clearly visible, both of them to disappear later at 110 µs. These two characteristic peaks imply high concentrations of CO, and it is therefore tempting to assign them to the two equilibrium states of heme’s inner cavity, with protonated and deprotonated propionate D, a mechanism more likely to occur with oxygen molecule as a ligand-substrate. However, the diffusion of CO from the inner to outer cavity may also take microsecond time scales, thus, the two peaks could be due to concentrations of CO in these two cavities. We can not exclude this mechanism, but we must point out that the wide distribution of COs in the

12134

J. Phys. Chem. B, Vol. 113, No. 35, 2009

Porrini et al. the confinement of COs in the cavity. From the crossings of the two pmf curves in both plots we might conclude that the potential barrier is about 2.4 kcal/mol for the cis conformation and 1.5 kcal/mol for the trans structure. However, these can not be real barriers, because we have found that there are pathways from which CO may rebind to CuB. Conclusions

Figure 8. Potential curves of mean force obtained by translating CO between copper and ring-D of heme a3: (a) cis and (b) trans ring-D propionate. The black curves denote CO approaching copper and red curves moving away from copper atom toward ring-D with the intermolecular interactions switched off (see text for details).

outer cavity found in our calculations do not support a sharp spectral band as has been seen in the FTIR spectra. Another alternative for the 2146 cm-1 peak is to be attributed to energy transfer between CO vibrational modes. The QM/ MM vibrational frequency analysis showed a few nearby CO harmonic vibrational frequencies. However, such a mechanism implies a very weak coupling term to satisfy a 50 microseconds time interval for energy transfer. Obviously, an unequivocal assignment of these two spectral features (2131/2146 cm-1) of ba3-CcO has to wait until MD calculations at microseconds time intervals as well as the study of quantum effects are feasible. Potential of Mean Force Calculations. Further insight in the B state of CO is provided by calculating potential of mean force (pmf) curves for cis and trans ring-D propionate shown in Figure 8: (a) for the closed gate structure and (b) for the open gate one. The curves are the result of thermodynamic integrations63 of values of mean force, calculated by freezing the CO carbon atom by means of two constrained distances: one corresponding to the pmf coordinate (between CO carbon and copper) and the other one between CO carbon and one reference atom of the protein matrix. For each isomer, two curves are computed: one starts from a configuration with CO bound to copper and it proceeds by translating CO along the direction of the internuclear axis of Cu-C bond and toward ring-D. The potential includes the Morse function, but no intermolecular interactions of carbon with copper environment. The second curve is calculated by placing CO in a position very close to ring-D of heme a3 and by translating it toward CuB but taking into account the intermolecular forces. Because in pmf calculations we do not compute absolute energies, we rely on the experimental result to estimate the energy gap between the two minima of the two pmf curves: being this gap equal to the free energy difference (∆G ) -RT ln Kx) between the two states in equilibrium (CO bound to copper and sheltered in the cavity), and from the results of ref 44 we set ∆G ) 0.82 kcal/mol. The mole fraction equilibrium constant, Kx, is evaluated from the concentrations of CO in the cavity (20%) and that one in CuB (80%) of the total concentration of dissociated CO. The curves in Figure 8 justify our previous conclusions. The flatness of the potential in the trans case of propionate D is in accordance to widely distributed COs, whereas for the cis case the deeper minimum in the long-range interaction explains

Different spectroscopic techniques covering a wide range of time scales have revealed several sites of cytochrome c oxidases, which can capture carbon monoxide after its photodissociation from heme a3. Vos and co-workers19 using the mid-infrared chirped-pulse upconversion technique time-resolved the dissociation process of CO from heme iron and its association to CuB by recording the vibrational signature of the ligand. They demonstrated that the total process is happening in less than 1 ps. On the other hand, Varotsis and co-workers21,22 using TRSSFTIR difference spectroscopy demonstrated that 15-20% of CO concentration is captured in a docking site inside the inner cavity of heme a3. The spectroscopic signature of the docking site is a spectral peak at 2131 cm-1 attributed to the vibrational frequency of CO in this equilibrium position. This peak survives up to 35 µs and it disappears in less than 110 µs with the simultaneous appearance of a new peak at 2146 cm-1, almost the frequency of the free CO (2143 cm-1). The two intracavity equilibrium states of CO assigned to 2131 cm-1 and 2146 cm-1, respectively, have been named B1 and B0. We have carried out molecular dynamics calculations following hundreds of trajectories up to 100 ps. We confirmed that the carbon monoxide is associated to CuB in less than 1 ps. Most interestingly, we do find that the equilibrium position of CO in the inner cavity is above the ring-D of heme a3 when the attached propionate is in trans conformation, whereas it is between ring-A and ring-D and closer to the binuclear center when propionate D is in the cis conformation. Because the protonation/deprotonation procedures in CcO enzymes are microseconds events, it is tempting to assign B1 and B0 states of CO to these two intracavity equilibrium geometries, although migration of the ligand from the inner to the outer cavity of heme a3 can not be excluded. Our calculations have mostly been performed on the EGEE European Grid, which provides thousands of computers and allows us to efficiently sample the phase space of the enzyme. In this work, we vary the potential function with a discontinues way by protonated/deprotonated species and switching off/on intermolecular interactions. Future work is planned to produce a smooth continuous potential model.64 Acknowledgment. We are grateful to Dr. Stamatis Stamatiadis and Manos Giatromanolakis for their computational assistance. V.D. thanks Dr. Victor Guallar for his hospitality at the Supercomputing CentersCentro Nacional de Supercomputacio´n, Barcelona, Spain, during which the QM/MM calculations with the Schro¨dinger suite of programs were carried out. The present work was financially supported by the European Union ToK Grant GRID-COMPCHEM (MTKD-CT-2005029583). We also acknowledge the EGEE European project for providing to us the necessary computational facilities, as well as the COST CMST D37 Action, GRIDCHEM. Supporting Information Available: Parameters for the CO interactions, movies of several trajectories, energy disposal after CO photodissociation, DFT calculations, dynamics at higher excitation energy, and QM/MM geometry optimization and

Heme Cavity Dynamics of Photodissociated CO calculation of CO vibrational frequencies. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Babcock, G. T.; Wikstro¨m, M. Nature 1992, 356, 301–309. (2) Michel, H.; Behr, J.; Harrenga, A.; Kannt, A. Annu. ReV. Biophys. Biomol. Struct. 1998, 27, 329–356. (3) Verkhovsky, M. I.; Jasaitis, A.; Vekhovskaya, M. L.; Morgan, J. E.; Wikstro¨m, M. Nature 1999, 400, 480–483. (4) Soulimane, T.; Buse, G.; Bourenkov, G. P.; Bartunik, H. D.; Huber, R.; Than, M. E. EMBO J. 2000, 19, 1766–1776. (5) Olson, J. S.; Phillips, G. N., Jr. J. Biol. Chem. 1996, 271, 17593– 17596. (6) Gibson, Q. H.; Regan, R.; Elber, R.; Olson, J. S.; Carver, T. E. J. Biol. Chem. 1992, 267, 22022–22034. (7) Quilli, M. L.; Li, T.; Olson, J. S.; Phillips, G. N., Jr.; Dou, Y.; Ikeda-Saito, M.; Elber, R.; Li, H.; Regan, R.; Gibson, Q. H. J. Mol. Biol. 1995, 245, 416–436. (8) Nutt, D. R.; Meuwly, M. Biophys. J. 2003, 85, 3612–3623. (9) Brunori, M. Biophys. Chem. 2000, 86, 221–230. (10) Frauenfelder, H.; Fenimore, P. W.; McMahon, B. H. Biophys. Chem. 2002, 98, 35–48. (11) Cohen, J.; Arkhipov, A.; Braun, R.; Schulten, K. Biophys. J. 2006, 91, 1844–1857. (12) Pinakoulaki, E.; Yoshimura, H.; Daskalakis, V.; Yoshioka, S.; Aono, S. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 14796–14801. (13) Blackmore, R. S.; Greenwood, C.; Gibson, Q. H. J. Biol. Chem. 1991, 266, 19245–19249. (14) Varotsis, C.; Woodruff, W. H.; Babcock, G. T. J. Am. Chem. Soc. 1989, 111, 6439–6440. (15) Gibson, Q. H.; Greenwood, C. Biochemistry J. 1963, 86, 541–554. (16) Greenwood, C.; Gibson, Q. H. J. Biol. Chem. 1967, 242, 1782– 1787. (17) Lambry, J. C.; Vos, M. H.; Martin, J. L. J. Phys. Chem. A 1999, 103, 10132–10137. (18) Liebl, U.; Lipowski, G.; Ne´grerie, M.; Lambry, J. C.; Martin, J. L.; Vos, M. H. Nature 1999, 401, 181–184. (19) Treuffet, J.; Kubarych, K. J.; Lambry, J.-C.; Pilet, E.; Masson, J.B.; Martin, J.; Vos, M. H.; Joffre, M.; Alexandrou, A. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 15705–15710. (20) Vos, M. H. Biochim. Biophys. Acta 2008, 1777, 15–31. (21) Koutsoupakis, C.; Soulimane, T.; Varotsis, C. J. Biol. Chem. 2003, 278, 36806–36809. (22) Koutsoupakis, C.; Soulimane, T.; Varotsis, C. J. Am. Chem. Soc. 2003, 125, 14728–14732. (23) Koutsoupakis, C.; Stavrakis, S.; Pinakoulaki, E.; Soulimane, T.; Varotsis, C. J. Biol. Chem. 2002, 277, 32860–32866. (24) Ruscio, J. Z.; Kumar, D.; Shukla, M.; Prisant, M. G.; Murali, T. M.; Onufriev, A. V. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 9204–9209. (25) Elber, E. R.; Karplus, M. J. Am. Chem. Soc. 1990, 112, 9161– 9175. (26) Petrich, J. W.; Lambry, J. C.; Kuczera, K.; Karplus, M.; Poyart, C.; Martin, J. L. Biochemistry 1991, 30, 3975–3987. (27) Meller, J.; Elber, E. R. Biophys. J. 1998, 74, 789–802. (28) A Resource for Studying Biological Macromolecules: http:// www.rcsb.org/. (29) PyMOL: http://pymol.sourceforge.net/. (30) Hunsicker-Wang, L. M.; Pacoma, R. L.; Chen, Y.; Fee, J. A.; Stout, D. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2004, 61, 340–343. (31) Tinker-Software Tools for Molecular Design: http://dasher.wustl. edu/tinker/.

J. Phys. Chem. B, Vol. 113, No. 35, 2009 12135 (32) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049–1074. (33) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Kleins, M. L. J. Chem. Phys. 1983, 79, 926–935. (34) Foloppe, N.; MacKerell, A. D. J. Comput. Chem. 2000, 21, 86– 104. (35) Daskalakis, V.; Pinakoulaki, E.; Stavrakis, S.; Varotsis, C. J. Phys. Chem. B 2007, 111, 10502–10509. (36) DMOL is a first-principles (ab initio) quantum chemistry software package. (37) Ahlrichs, R.; Ba¨r, M.; Ha¨ser, M.; Horn, H.; Ko¨lmel, C. Chem. Phys. Lett. 1989, 162, 165–169. (38) DL_POLY: http://www.cse.scitech.ac.uk/ccg/software/DL_POLY/. (39) Enabling Grid for E-sciencE (EGEE): http://www.eu-egee.org/. (40) HellasGrid: http://www.hellasgrid.gr/. (41) CompChem: http://compchem.unipg.it/. (42) South East Europe (See): http://www.egee-see.org/. (43) Daskalakis, V.; Giatromanolakis, M.; Porrini, M.; Farantos, S. C.; Gervasi, O. 2009, to be submitted. (44) Pinakoulaki, E.; Ohta, T.; Soulimane, T.; Kitagawa, T.; Varotsis, C. J. Biol. Chem. 2004, 279, 22791–22794. (45) Luna, V. M.; Chen, Y.; Fee, J. A.; Stout, C. D. Biochemistry 2008, 47, 4657–4665. (46) Koutsoupakis, C.; Stavrakis, S.; Soulimane, T.; Varotsis, C. J. Biol. Chem. 2003, 278, 14893–14896. (47) Koutsoupakis, C.; Pinakoulaki, E.; Stavrakis, S.; Daskalakis, V.; Varotsis, C. Biochim. Biophys. Acta 2004, 278, 347–352. (48) Wikstro¨m, M.; Ribacka, C.; Molin, M.; Laakkonen, L.; Verkhovsky, M.; Puustinen, A. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 10478–10481. (49) Branden, G.; Branden, M.; Schmidt, B.; Mills, D. A.; FergusonMiller, S.; Brzezinski, P. Biochemistry 2005, 44, 10466–10474. (50) Hosler, J. P.; Ferguson-Miller, S.; Mills, D. A. Annu. ReV. Biochem. 2006, 75, 165–187. (51) Mills, D. A.; Geren, L.; Hiser, C.; Schmidt, B.; Durham, B.; Millett, F.; Ferguson-Miller, S. Biochemistry 2005, 44, 10457–10465. (52) Seibold, S. A.; Mills, D. A.; Ferguson-Miller, S.; Cukier, R. I. Biochemistry 2005, 44, 10475–10485. (53) Xu, J.; Sharpe, M. A.; Qin, L.; Ferguson-Miller, S.; Voth, G. A. J. Am. Chem. Soc. 2007, 129, 2910–2913. (54) Kannt, A.; Roy, C.; Lancaster, D.; Michel, H. Biophys. J. 1998, 74, 708–721. (55) Behr, J.; Hellwig, P.; Ma¨ntele, W.; Michel, H. Biochemistry 1998, 37, 7400–7406. (56) Wikstro¨m, M.; Verkhovsky, M. I.; Hummer, G. Biochim. Biophys. Acta 2003, 1604, 61–65. (57) Siegbahn, P. E. M.; Blomberg, M. R. A.; Blomberg, M. L. J. Phys. Chem. B 2003, 107, 10946–10955. (58) Lim, M.; Jackson, T. A.; Anfinrud, P. A. Nat. Struct. Biol. 1997, 4, 209–214. (59) Sagnella, D. E.; Straub, J. E.; Jackson, T. A.; Lim, M.; Anfinrud, P. A. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 14324–14329. (60) Anfinrud, P. A.; de Vivie-Riedle, R.; Engel, V. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 8328–8329. (61) Nienhaus, K.; Olson, J. S.; Franzen, S.; Nienhaus, G. U. J. Am. Chem. Soc. 2004, 127, 40–41. (62) Molecular Modeling Platform: http://www.schrodinger.com/. (63) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: New York, 1996. (64) Mavrandonakis, A.; Farantos, S. C.; Froudakis, G. J. Comput. Theor. Nanosci. 2009, 6, 880–885. (65) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33–38.

JP904466N