Molecular Oxygen Activation and Proton Transfer Mechanisms in

May 13, 2009 - This amino acid, in concert with an active-site water network, catalyzes a facile protonation of the peroxo intermediate and offers a s...
1 downloads 11 Views 6MB Size
8170

J. Phys. Chem. B 2009, 113, 8170–8182

Molecular Oxygen Activation and Proton Transfer Mechanisms in Lanosterol 14r-Demethylase Catalysis Kakali Sen and John C. Hackett* Institute for Structural Biology and Drug DiscoVery, Virginia Commonwealth UniVersity, 800 East Leigh Street, Richmond, Virginia 23219 ReceiVed: March 31, 2009

The CYP51 lanosterol 14R-demethylases are evolutionarily ancient enzymes ubiquitously distributed throughout the biological domains. The experimental X-ray crystal structure of Mycobacterium tuberculosis (Mtb) CYP51 is the first of an enzyme capable of catalyzing inert C-C bond cleavage. Amino acid sequence comparisons of CYP51 family members with other members of the CYP superfamily reveal the almost universally conserved “acid-alcohol” pair, putatively involved in proton transport and O2 activation, is replaced with a His-Thr dyad. In this study, extended molecular dynamics (MD) simulations and hybrid quantum mechanics/molecular mechanics calculations (QM/MM) are applied to characterize reactive oxygen intermediates and to unravel mechanisms of O2 activation Vis-a`-Vis proton transport for this important enzyme. MD confirms stable His259δH+-Thr260OH-O2 (Mtb numbering) hydrogen bonding early in the simulations, suggesting these amino acids could function similarly to the Asp251-Thr252 pair in CYP101. QM/MM calculations support this dyad competently catalyzes the peroxo to Compound 0 (Cmpd 0) reaction, albeit an endothermic homolytic O-O scission mechanism affording Compound I (Cmpd I) was identified. Disruption of the His259H+Thr260OH hydrogen bond in MD simulation divulges a second previously unidentified hydrogen-bond network, including three water molecules linking Glu173 in the CYP51 F-helix to the distal O2 atom. Expansion of the QM region to contain these atoms unveils an unprecedented triradicaloid electronic structure of the peroxo intermediate characterized by spin polarization to the Glu173 side chain, attributable to the protein electrostatic environment. This amino acid, in concert with an active-site water network, catalyzes a facile protonation of the peroxo intermediate and offers a series of redundant heterolytic and homolytic mechanisms, affording exothermic formation of the ultimate oxidant Cmpd I. In summary, this study highlights the importance of the protein electrostatic environment to tune the electronic structure of CYP catalytic intermediates in addition to Cmpd I and illustrates the diversity of proton transport pathways available to these enzymes to drive catalysis. Introduction The cytochromes P450 (CYPs) use O2 and two reducing equivalents to perform regioselective and regiospecific oxygen insertion into a myriad of organic compounds.1-5 These enzymes are involved in the biosynthesis of endogenous compounds and are responsible for the oxidation of ∼70% of the therapeutic drugs and those in the development pipeline.6-8 The widely quoted CYP catalytic cycle involves many intermediates, several of which have been elusive to experimental characterization. On substrate binding, the pentacoordinate Fe3+ enzyme (Scheme 1, A) transitions from a low to high spin, increasing its redox potential. Upon reduction, the Fe2+ CYP ligates O2 resulting in oxyferrous CYP (Scheme 1, B). This is the final quasi-stable and observable intermediate in the catalytic cycle.9 A second single-electron reduction results in the formation of the peroxo intermediate (Scheme 1, C). Experimental observation of this species has proven difficult, and its existence as a stable intermediate is controversial. Denisov and co-workers have recently reported observation of this intermediate after radiolytic reduction in the Asp251Asn mutant of CYP101 by resonance Raman spectroscopy.10 A number of experimental11 and theoretical studies12-18 of the CYP101 system suggest the second * Corresponding author. Telephone: 804.828.5679. Fax: 804 0.827.3664. E-mail: [email protected].

electron and first proton transfer occur in a concerted manner resulting in the hydroperoxo species, also known as Compound 0 (Cmpd 0, Scheme 1, D). Further protonation of Cmpd 0 results in O2-bond scission affording Compound I (Cmpd I, Scheme 1, E). This species is widely accepted to be the reactive species in the CYP oxidation chemistry, although its direct experimental observation has proven elusive. Alignments of CYP sequence and structure reveal the acidalcohol pair of amino acids are conserved in the active sites of many of these enzymes.19 The amino acid occupying the alcohol position is typically threonine and less frequently serine. The Asp251-Thr252 pair of the prototypical CYP101 enzyme has been the subject of extensive mutagenesis studies. Biophysical characterization of these mutants supports a role for this pair of residues in the proton-delivery machinery involved in O2 activation.20,21 Whether the role of the “conserved” Thr is to directly protonate the peroxo species or to stabilize a network of water molecules serving this purpose is controversial. Proton inventory analysis of CYP101 catalysis indicates at least two protons are in flight in the proton-transfer transition state, supporting the Asp251-Thr252 dyad works, in concert with active-site water molecules, to drive the formation of catalytic intermediates.11 The O2-consumption rates of a Thr252Ala mutant are similar to a wild-type; however, 5-exo-hydroxy camphor formation is significantly diminished, apparently as a

10.1021/jp902932p CCC: $40.75  2009 American Chemical Society Published on Web 05/13/2009

CYP51 Molecular Oxygen Activation

Figure 1. Overlay of the Mtb CYP51 (green) and CYP101 (magenta) active sites highlighting the positions of residues putatively involved in O2 activation. Water molecules connecting the putative proton donor Glu366 to the distal O2 atom in CYP101 are represented by magenta spheres.

SCHEME 1

result of unproductive shuttling of reducing equivalents to form hydrogen peroxide and water.20 Kimata and co-workers replaced Thr252 with methoxythreonine, and the resulting mutant maintained monooxygenase activity, suggesting that the presence of a free hydroxyl group in this position is not a prerequisite for O2 activation.21 An alternative route for proton flow in the oxygen-activation process was revealed in the X-ray crystal structure of oxyferrous CYP101, characterized by a closed chain of four water molecules connecting the Thr252 hydroxyl group to Glu366.9 This crystal structure prompted a number of computational studies to assess roles for the acid-alcohol pair and crystallographically identified the water chain in O2 activation12,13,16 (Figure 1). These computational studies, complemented by decades of elegant experiments, have strongly invoked the acid-alcohol pair in facilitating CYP catalysis. However, questions regarding redundancy of proton-transfer networks and operational differences of these networks at different steps of the catalytic cycle remain to be resolved. The CYP51 enzymes are ancient members of the cytochrome P450 family with representatives from all biological domains.22 These enzymes catalyze the oxidative demethylation of 14Rmethylated steroids with concomitant introduction of a double bond in three oxidative steps, each consuming a single

J. Phys. Chem. B, Vol. 113, No. 23, 2009 8171 equivalent of NADPH and O223,24 (Scheme 2). The first two steps apparently proceed by the Cmpd I-mediated hydrogenatom abstraction/hydroxyl radical rebound mechanism originally proposed by Groves and co-workers.25,26 The final step results in cleavage of the C14-C32 bond with stereoselective removal of the 15R-hydrogen, resulting in the introduction of a 14,15 double bond and the release of formic acid (Scheme 2). Based on the interpretation of biochemical experiments utilizing isotopically labeled substrates, the peroxo species has been proposed to mediate the final catalytic step.24 However, the identity of the oxygen intermediate involved in this step is controversial. A QM study using truncated active-site models of the similar C-C bond cleavage reaction catalyzed by CYP19 suggested an alternative mechanism involving Cmpd I.27 Comparison of CYP51 primary structures reveals the acidalcohol pair is replaced by a His-Thr dyad.22 The X-ray crystal structure of Mycobacterium tuberculosis (Mtb) CYP51,28,29 the first of a sterol biosynthetic enzyme, when compared with other CYPs confirms this dyad occupies the same position as the acidalcohol pair in the active site. Furthermore, a nonionizable Gln residue occupies the same position as Glu366 in CYP101. Figure 1 illustrates the relative positions of the His-Thr dyad to residues implicated in proton transfer and the crystallographically observed water network. The absence of the putative proton relay machinery in the active site of CYP51 raises intriguing questions about its O2 activation mechanism and how the lifetimes of intermediates in the catalytic cycle are modulated. Experimental evidence supporting the use of multiple oxidants and an unusual active-site architecture in the context of proton delivery, conserved among the CYP51 family, have spawned our interest in the mechanisms of formation and the properties of CYP51 catalytic intermediates. In this study, computational methods are applied to unravel the mechanisms of proton transport and the reactive oxygen intermediate formation in Mtb CYP51 catalysis. Molecular dynamics (MD) simulations are used to study the evolution of hydrogen-bond networks and hydration of the oxyferrous CYP51 active site. Utilizing representative configurations derived from these simulations, hybrid quantum mechanics/molecular mechanics (QM/MM) methods are applied and illustrate a novel role for an acidic active-site residue (Glu173) in the CYP51 F-helix, apparently conserved in structure among many CYPs to tune the electronic structure of the peroxo intermediate and shuttle protons necessary to drive formation of CYP catalytic intermediates. Computational Methods System Preparation and Molecular Dynamics Simulation. The initial coordinates for Mtb CYP51 were derived from the X-ray crystallographic structure of this protein bound to the steroid ligand, estriol (PDB ID 1X8V).29 The propKa module of the PDB2PQR suite of programs30-32 and visual inspection of side chain local environments were used to adjust the protonation states of ionizable residues. The initial binding mode for the lanosterol substrate was generated by fitting atomic coordinates of the steroid nucleus to the corresponding crystallographic coordinates of estriol. CYP51 systems were solvated in a 78.7 × 91.2 × 102.5 Å3 box of TIP3P water including 36 sodium and 22 chloride ions. Explicit solvent molecular simulations utilizing Langevin dynamics and periodic boundary conditions were performed using the NAMD 2.6 suite of programs33 and the CHARMM27 force field.34 Molecular mechanics parameters for lanosterol were those derived by Cournia and co-workers.35 With the CYP51 protein and ligands restrained, the solvent and ions were subjected to 2500

8172

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

Sen and Hackett

SCHEME 2

minimization steps followed by constant volume (NVT) simulation at 100, 150, 200, 250, and 300 K each for 10 ps. The restraint on the protein and ligands was subsequently released, and the minimization and annealing protocol was repeated. This procedure was followed by a constant pressure (NPT) molecular dynamics simulation for 20 ns with a 2 fs time step. Constant pressure simulations were performed using the modified Nose´-Hoover method, in which Langevin dynamics is used to control fluctuations in the barostat.36,37 Partitioning into QM and MM Systems. Coordinates for the CYP51 protein, associated ligands, and water molecules within 5.0 Å of these species were extracted from MD trajectory points resulting in systems containing approximately 14 000 atoms for QM/MM treatment. The extracted coordinates were subjected to 11 000 geometry optimization steps (1000 steepest descent, 3000 conjugate gradient, and 7000 adopted basis Newton-Raphson) using the CHARMM molecular mechanics package.38 Systems were partitioned into a region for QM treatment containing the O2-coordinated iron porphyrin, relevant amino acid side chains, and water molecules. Propionate, methyl, and vinyl side chains of the iron protoporphyrin IX were excluded. The lanosterol substrate was excluded from the QM region; however, it maintained a catalytically competent binding mode with the 14R-methyl group oriented toward the O2 ligand, as shown in Figure 2. In the QM/MM calculations, the iron protoporphyrin IX, associated O2 ligand, and amino acids with their side chains in the QM region, as well as their N- and C-terminal neighbors, were unrestrained during the geometry optimization, while coordinates of the remaining atoms were frozen. Hybrid QM/MM Calculations. Hybrid QM/MM calculations were performed with the ChemShell package,39 which integrates the TURBOMOLE40 and DL-POLY41 packages for density functional theory and molecular mechanics calculations, respectively. Unless otherwise noted, the electrostatic embedding scheme with the charge shift correction was used to represent the electrostatic interaction between the QM region and protein environment.42 In this scheme, the MM region appears as a distribution of partial charges from the CHARMM force field to polarize the QM Hamiltonian. In QM/MM geometry optimizations, the QM region was treated with the unrestricted hybrid density functional B3LYP,43-45 while the CHARMM22 force field was used to treat the MM region. The Ahlrichs VTZ46 for Fe, 6-31+G* for N, O, and S, and 6-31G* for C and H basis sets were used in QM/MM geometry optimizations. Reactants, approximate transition states, and products on the reaction coordinate were subjected to single point energy calculations using the Wachters +f basis set47,48 for Fe and the TZVP49 basis set for the remaining atoms. The levels of theory used for geometry optimization and single point energy calcula-

tions are referred to as B1 and B2 hereafter, respectively. Population analyses were carried out with the natural population analysis (NPA) method50 in the TURBOMOLE 5.91 suite of programs. NPA spin densities and charges reported in the text are derived from B2 wave functions at the B1-optimized geometry. The size of the QM systems and their numerous degrees of freedom described herein made identification and characterization of authentic transition states computationally intractable. To map the potential energy surfaces and to calculate approximate barriers for possible proton-transfer pathways, a series of restrained geometry optimizations along defined reaction coordinates were performed. After identification of an approximate transition state, geometries of the corresponding reactant and product complexes were obtained by releasing restraints at the extremes of the reaction coordinate followed by full geometry optimization. Electron paramagnetic resonance (EPR) studies of reduced oxyferrous CYP101,51 as well as density functional theory calculations of Harris and Loew,52 using a truncated model of the active site, support a doublet ground state for the peroxo intermediate. A QM/MM study of this intermediate in CYP101 predicted only a stable doublet state, and the calculations

Figure 2. Active site of oxyferrous CYP51. The lanosterol molecule is shown in yellow.

CYP51 Molecular Oxygen Activation

J. Phys. Chem. B, Vol. 113, No. 23, 2009 8173

Figure 3. Ala256CO-Thr260OH (green), Thr260OH-O2 (red), and Thr260OH-His259 δH (blue) hydrogen-bond distances in production NPT MD simulations initiated with the neutral (A) and cationic (B) protonation states of His259 and the δNH oriented toward the Thr260OH.

supported higher spin states (quartet, sextet) that were dissociative, in accord with EPR studies of this intermediate in the same enzyme.53 With the exception of the singlet oxyferrous intermediate, all calculations were performed on the doublet potential energy surface. Results and Discussion Molecular Dynamics Simulation. A number of elegant experimental and theoretical studies have implicated a “conserved” acid-alcohol pair of amino acids as important in shuttling protons for the O2 activation in cytochrome P450 enzymes. In Mtb CYP51, as well as in CYP51 enzymes across the biological domains, the acidic residue in the I-helix is replaced with histidine. This amino acid maintains unique properties among the naturally occurring amino acids, including a pKa near physiological pH, tunable by the protein environment. The relative positions of the CYP51 F- and I-helices, as well as important residues in the oxyferrous form of this enzyme with lanosterol bound, are shown in Figure 2. In the Mtb CYP51 crystal structure with estriol, the plane of the His259 imidazole ring is nearly parallel to the heme plane. We hypothesized this side chain could deprotonate an upstream proton source, such as solvent water. Subsequently, the histidinium cation could protonate the peroxoferric intermediate via the Thr260 hydroxyl. For the success of this proton-transfer network, a His259Thr260OH hydrogen bond is a prerequisite. Therefore, the His259 side-chain orientation was manipulated to generate multiple initial positions in which the δNH and εNH were oriented toward the Thr260 oxygen. Manipulation of the His259 side-chain orientation and protonation state resulted in four configurations of the enzyme, each of which was subjected to 20 ns of NPT MD simulation with O2 and lanosterol bound. To assess the stability of the His259(δ/ε)H-Thr260OH-O2 networks, hydrogen-bond distances were monitored over the NPT MD simulations. Distances for simulations initiated with the δNH oriented toward Thr260OH are shown in Figure 3. Equivalent distances for simulations initiated with the His259εH oriented toward the Thr260OH are provided in the Supporting Information. In both protonation states of His259, the His259δHThr260OH-O2 hydrogen bonds are maintained early in the MD simulation. During this time, the average His259δH-Thr260OH hydrogen-bond distances are 2.0 and 2.4 Å for the cationic and neutral states, respectively. In the former simulation, transient hydrogen bonding with the backbone carbonyl of Ala256 in the CYP51 I-helix is observed (Figure 3B). Harris and Loew made a similar observation in MD simulation of CYP101 in which transient, bimodal hydrogen bonding of Thr252 with Gly248 and the distal O2 atom occurred.54 In CYP51, this phenomenon is manifested by correlated changes in the

Thr260OH-O2 and Thr260OH-Ala256CO hydrogen-bond distances. In contrast, when His259 is neutral, no transient hydrogen bonding to Ala256 is observed. Instead, the Thr260OH continuously hydrogen bonds to the distal O2 atom over ∼5 ns (Figure 3A). It is evident from Figure 3 that the His259δH-Thr260OH hydrogen bond is disrupted beyond 4 ns when His259 is cationic and beyond 5 ns when this residue is neutral. From these points through the end of the simulations, there are significant increases in the average His259δH-Thr260OH distances to 6.1 and 5.2 Å in those initiated with cationic and neutral His259, respectively. Moreover, the Thr260OH-O2 distance increases to approximately 4.0 Å in each case. Correlated with these increases in the average Thr260OH-O2 distance is transient hydrogen bonding of the Thr260OH with the Ala256 backbone carbonyl. To our surprise, disruption of the His259δH+Thr260OH-O2 (Figure 3B) occurs in concert with reorientation of the histidinium side chain and influx of water molecules into the active site. This influx results in the formation of a hydrogenbond network consisting of an average of three or four water molecules connecting His259, Thr260, O2, and Glu173 (Figure 2). Three water molecules bridge the distal O2 atom with Glu173, a fourth water molecule bridges the His259 δH and central water molecule of the aforementioned triad. In the analogous simulation initiated with neutral His259, influx and organization of a comparable number of water molecules are not observed. In MD simulations initiated with the His259 εNH oriented toward Thr260OH, no influx of water is observed. Alternatively, pairwise side-chain interactions preventing water influx develop. In the cationic case, the histidinium side chain forms a stable ion pair with Glu173, blocking access of water to the active site. In the neutral case, the His259 εNH forms an apparent NH-π interaction with the Phe255 side chain, also preventing solvent entry into this cavity. Distances between the Glu173/ His259 and Phe255/His259 pairs over the relevant MD trajectories supporting evolution of these interactions are reserved for Supporting Information. The MD simulation initiated with a His259H+ cation oriented favorably for hydrogen bonding to the Thr260 hydroxyl (Figure 3B) provides two configurations of the CYP51 active site with the potential to shuttle protons to the peroxo intermediate. To confirm the viability of these proton-transfer networks, QM/ MM methods were applied to investigate the potential energy landscapes connecting CYP51 catalytic intermediates and their properties in each configuration of the active site. Due to the computational expense of the QM/MM method applied in this study and the relative stability of the hydrogen-bond networks established over the two time frames, we have extracted a single

8174

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

Sen and Hackett

Figure 4. B1:CHARMM-optimized geometries of the Mtb CYP51 active site characterized by the His259H+-Thr260OH (A) and Glu173-(H2O)4 (B) hydrogen-bond networks. Bond distances for the corresponding oxyferrous and peroxo species (in parentheses) are listed in Å.

representative frame from each of the 0-4 ns and 4-20 ns regions of the MD simulation. The geometries associated with these points on the MD trajectory were prepared as described in Computational Methods and used as the initial geometries to study the transformation of oxyferrous CYP51 to the Cmpd I species. QM/MM Evaluation of His259-Thr260-Catalyzed Proton Transfer. It is widely suspected that initiation of proton flow is essentially concerted with reduction of the oxyferrous species in CYP enzymes. Resonance Raman and EPR studies support reduction of the oxyferrous complex is a prerequisite for protonation of the distal O2 atom.55 Moreover, utilizing radiolytic reduction at cryogenic temperatures, it was shown that oxyferrous CYP101 is directly converted to Cmpd 0 as a result of single-electron reduction.56-58 Geometries of the oxyferrous and peroxo species of each active site configuration were fully optimized at the B1 level of theory and are shown in Figure 4. Considering the His259δH+-Thr260OH configuration shown in Figure 4A, reduction of the oxyferrous to the peroxo species results in lengthening of the Fe-O (0.19), Fe-S (0.11), and O-O (0.04 Å) bonds. This event also influences hydrogen bonds distant from the catalytic core. Consistent with stronger hydrogen bonding in the peroxo species, there are marked decreases in the Thr260OH-O2 (0.30 Å) and His259δH+Thr260OH (0.08 Å) hydrogen bonds (Figure 4A). Changes in geometry and electronic structure observed in the peroxo species are consistent with, and natural orbital occupations support, single-electron occupation of an antibonding orbital primarily localized to the O2 unit (π*OO) (see Supporting Information). These results are similar to those first described in a study by Harris and others utilizing truncated active-site models, demonstrating reduction of the oxyferrous intermediate to the peroxo species, which fills a π*yz orbital leading to longer Fe-O and Fe-S bonds.52,59 Moreover, the localization of the excess electron to the O2 unit is consistent with the electronic structure of this intermediate recently predicted by QM/MM computation in CYP101.53 Defining the reaction coordinate (ξThr-Od) as the distance between the Thr260 hydroxyl proton and the distal O2 atom (Od), a series of restrained geometry optimizations were

performed to map the reaction coordinate for conversion of the peroxo intermediate to Cmpd 0 with His259H+ as the terminal proton donor. Fully optimized geometries of the approximate transition state, corresponding reactant, and product species, as well as their spin densities and charges, are shown in Figure 5. The approximate transition state (ξThr-Od ) 1.3 Å) lies 4.1 (4.8) kcal/mol above the corresponding peroxo species. Geometry optimization of species with ξThr-Od < 1.3 Å collapses to Cmpd 0 with concomitant reprotonation of Thr260 by the His259 δH+. Formation of Cmpd 0 is -24.0 (-22.8) kcal/mol exothermic, relative to the peroxo species. Full geometry optimization of species with ξThr-Od > 1.3 Å collapses to the corresponding peroxo species. The energetics described here for His259Thr260-mediated protonation are similar to those obtained by Wang and co-workers in a QM/MM study of the Asp251-H2OThr252 proton shuttle in CYP101.15 These calculations demonstrated that the dyad in CYP101 could catalyze a facile, concerted proton transfer proceeding with a somewhat smaller barrier of 1.2 kcal/mol and a similar -22.0 kcal/mol exothermicity. These results support the His259H+-Thr260OH structural variation in the CYP51 also offers a concerted, energetically competitive pathway for Cmpd 0 formation in the CYP51 catalytic cycle. A similar series of restrained geometry optimizations were performed to assess the ability of the His259H+-Thr260OH to protonate Cmpd 0, hypothetically resulting in heterolytic scission of the O-O bond and formation of Cmpd I. Transfer of the Thr260 hydroxyl proton to the distal O2 atom through a series of restrained geometry optimizations results in a continuous energy increase to 32 kcal/mol relative to the Cmpd 0 equilibrium geometry at the B1 level of theory. Release of the applied restraints followed by full geometry optimization failed to produce Cmpd I and water. Instead, geometry optimization returned these to the corresponding Cmpd 0 species. Similar problems were encountered in QM/MM calculations focused on the Asp251-H2O-Thr25216 and the Glu36613 protonation channels. However, in agreement with QM/MM studies of the analogous proton channel in CYP101, these theoretical methods are unable to support the existence of a heterolytic mechanism

CYP51 Molecular Oxygen Activation

J. Phys. Chem. B, Vol. 113, No. 23, 2009 8175

Figure 5. B1:CHARMM fully optimized geometries, energies (relative to the peroxo species; kcal/mol), spin densities, and charges (parentheses) along the reaction coordinate for His259H+-Thr260OH mediated transformation of the peroxo species to Cmpd 0. B2:CHARMM//B1:CHARMM energies are in parentheses. Important bond distances are listed in Å.

for the Cmpd 0 to Cmpd I transformation catalyzed by the His259H+-Thr260 dyad. We also explored the possibility this transformation could be achieved via a mechanism initiated by homolytic cleavage of the bond between the proximal and distal O2 atoms (Op-Od). Although these efforts initially identified an intermediate complex (IC) consisting of Cmpd I and an OH moiety 13.0 (13.5) kcal/mol above the corresponding Cmpd 0 species, the ultimate formation of Cmpd I was 2.3 (4.8) kcal/mol endothermic. In light of these results, we have reserved those data for the Supporting Information and focused the remainder of the paper on the apparent proton shuttle evolving beyond ∼4 ns in MD simulations. QM/MM Evaluation of Glu173-(H2O)4-Catalyzed Proton Transfer. The hydrated active-site configuration evolving beyond 4 ns of MD simulation offered a second putative protontransfer network, employing Glu173 in the CYP51 F-helix as a potential proton donor. The presence of this acidic residue appropriately positioned for reprotonation during catalysis, combined with the evolution of a stable hydrogen-bonded network of water molecules, all in the context of a unique activesite architecture, led us to consider the possibility for a novel proton shuttle in CYP51. In pursuit of QM/MM evaluation of this presumed shuttle, expansion of the QM region was necessary. In addition to the Cys394 methylene thiolate, iron porphyrin, and O2, the Glu173 side chain, complete His259H+, and Thr260 were included (Figure 4B). Considering the primary focus of these calculations was to evaluate Glu173 as the proton source, the carboxylate side chain was protonated, and four water molecules linking this amino acid to the distal O2 atom and His259δH+ were also included in the QM region. The propKa pKa estimate for Glu173 is 8.3, consistent with ∼90% of these residues in the unionized form at physiological pH. Aware of the possible operation of two His259 charge states in the presence of neutral Glu173, we fully optimized the oxyferrous and peroxo systems with the His259 positively charged, and in the neutral protonation state with the remaining proton located on the imidazole εN. In both the oxyferrous and peroxo states of the enzyme, Glu173 remained protonated in the presence of the histidinium cation. However, deprotonation of the His259

TABLE 1: 1AFeO2 (AFeO2) B2:CHARMM//B1:CHARMM NPA Group and Atomic Spin Densities of the Peroxo Intermediate for the Glu173-(H2O)4 and Expanded QM Region

Fe Op Od porphyrin Glu173

E173

E173 + K167, Y169, R174, N436

-0.94 (0.87) 0.47 (0.50) 0.51 (0.47) -0.01 (0.00) -0.99 (-0.92)

0.94 (-0.91) -0.47 (-0.50) -0.51 (-0.48) 0.01 (-0.03) 0.97 (0.97)

at the δNH resulted in a spontaneous shift of the Glu173 proton via the water network to reprotonate the His259 δN. We concluded that, in this configuration of the enzyme active site, the existence of the nearby histidinium cation is a prerequisite for Glu173 to be a competitive source for distal O2 protonation. As a consequence, the protonated, cationic state His259H+ was used in all subsequent calculations. In the B1-optimized geometry, with the exception of the Fe-O bond, geometric changes observed in the His259H+Thr260OH configuration upon single-electron reduction are abrogated in the hydrated active site (Figure 4B). In contrast to the His259H+-Thr260OH configuration, initial calculations converged to a triradicaloid peroxo species characterized by two electrons localized to Fe and O2 in a triplet configuration (3AFeO2), while a third electron of antiparallel spin was localized to the protonated Glu173 side chain. A spin-flip strategy was used to converge a second, nearly degenerate, state characterized by an open-shell singlet configuration on Fe and O2 (1AFeO2) 0.7 kcal/mol below the 3AFeO2 state. Relevant spin densities for these states are listed in Table 1. The triradicaloid electronic structure for the peroxo intermediate is novel, so in an effort to verify this observation was not an artifact of the B3LYP functional; single point calculations with the 3AFeO2 state were performed with several additional functionals and the B2 basis sets. The spin densities for these levels of theory are listed in Table 2, all of which reflect the original B3LYP result (Table 1). The only variation occurred with the BHLYP60 functional. While maintaining the triradicaloid configuration, spin density is shifted

8176

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

Sen and Hackett

TABLE 2: 3AFeO2 B2:CHARMM//B1:CHARMM NPA Group and Atomic Spin Densities of the Peroxo Intermediate for the Glu173-(H2O)4 System Using Several Density Functionalsa Fe Op Od porphyrin Glu173 His259 a

B3LYP

BLYP

PBE

TPSS

BHLYP

PBE0

TPSSh

0.87 (0.70) 0.50 (0.51) 0.47 (0.41) 0.00 (-0.12) -0.92 (-0.24) 0.00 (-0.44)

0.75 (0.63) 0.46 (0.45) 0.35 (0.31) 0.00 (-0.06) -0.70 (-0.25) 0.00 (-0.36)

0.78 (0.65) 0.45 (0.45) 0.35 (0.31) -0.01 (-0.09) -0.70 (-0.25) 0.00 (-0.36)

-0.54 (0.19) 0.45 (0.46) 0.42 (0.33) 0.01 (-0.03) 0.74 (-0.22) 0.00 (0.31)

-0.02 (0.01) 0.55 (0.60) 0.45 (0.38) 0.93 (0.00) -1.00 (0.00) 0.00 (0.27)

0.94 (0.83) 0.49 (0.52) 0.49 (0.44) 0.02 (-0.19) -0.99 (-0.23) 0.00 (-0.48)

0.80 (-0.45) 0.48 (0.52) 0.44 (0.41) 0.00 (0.11) -0.82 (0.25) 0.00 (0.40)

Values in parentheses were computed in the absence of the partial charge distribution.

from Fe to the porphyrin π system. To evaluate the role of the protein electrostatic environment represented by the CHARMM force field, the partial charge distribution was removed, and the single point calculations were repeated with this collection of density functionals. At the B3LYP, BLYP,44,61-63 PBE,64,65 PBE0,66 and TPSSh67 levels of theory, the QM system remains triradicaloid, although some spin density shifts to the His259 imidazolium (Table 2). This spin density shift supports the necessity of the protein electrostatic environment to maintain the triradicaloid electronic structure of the peroxo intermediate. The protonated Glu173 carboxyl is adjacent to two positively charged amino acids (Lys167 and Arg174) and accepts hydrogen bonds from the Tyr169 hydroxyl and the Gln436 amide side chain (Figure 2). To verify accumulation of the odd electron was not an artifact of the QM boundary selection or an inappropriate treatment of the local electrostatic environment, the QM region around the Glu173 was expanded to include the side chains of these residues. Single point calculations at the B2 level of theory converged to the 1AFeO2 state with a third electron exclusively localized to Glu173. The spin-flip procedure was used to converge upon the corresponding 3AFeO2 state, which lies 0.9 kcal/mol above 1AFeO2. Relevant group and atomic spin densities derived from calculations containing these additional residues are listed in Table 1. The spin density distribution and 1 AFeO2-3AFeO2 energetic spacing characteristic of these triradicaloid states are unperturbed by introducing the local electrostatic and hydrogen-bond environment. Futhermore, the lack of significant spin polarization beyond the Glu173 side chain validates our initial selection of the QM boundary of the hydrated system. Unless noted otherwise, the immediate Glu173 local environment (Lys167, Tyr169, Arg174, and Gln436) is partitioned into the MM layer. To evaluate the possible role of Glu173 as the proton source for conversion of the peroxo species to Cmpd 0 via the hydrogen-bonded network of water, an analogous series of restrained geometry optimizations were performed to map the potential energy surface for this process. Realizing the ultimate proton donor to the distal O2 atom could either be water or Thr260, two reaction coordinates were considered. First, assuming the water molecule C to be the ultimate proton donor (See Figure 6), a series of restrained geometry optimizations were performed along ξH2O-Od. The fully optimized approximate transition state, reactant, and product complex in the peroxo to Cmpd 0 transformation are shown in Figure 6 with their spin densities and charges listed. On the 1AFeO2 surface, the Glu173(H2O)4-mediated formation of Cmpd 0 is facile, with the approximate transition state (ξH2O-Od ) 1.5 Å) lying 1.8 kcal/ mol above the peroxo species at the B1 level of theory. Release of the restraint when ξH2O-Od > 1.5 Å, followed by full geometry optimization, returned to the peroxo intermediate. Conversely, restraint release and full geometry optimization of the species with ξH2O-Od < 1.5 Å collapsed to the Cmpd 0 species. This occurred with concerted reprotonation of the terminal water

molecule by the Glu173 proton via water molecules A and B. This process is considerably more exothermic, greater than -100 kcal/mol below the peroxo intermediate, than when mediated by the His259H+-Thr260 dyad. The relative energetics for these species on both the 1AFeO2 and 3AFeO2 surfaces were evaluated with single point calculations at the B2 level of theory (Figure 6, parentheses). Protonation of the peroxo intermediate proceeds with a negligible barrier on both spin surfaces and is greater than -100 kcal/mol exothermic. If the QM region is expanded to include the residues in the Glu173 local environment, the reaction also proceeds with a negligible barrier and is -90 kcal/ mol exothermic (see Supporting Information). If the partial charge distribution contributed by the MM force field is removed, resulting in an environment akin to the gas phase, the reaction is -85 kcal/mol exothermic. Taken together, the results support the predominant contribution to this large exothermicity arises from altered bonding and interactions within the QM region. The remaining energetic contribution was mapped to the local environment surrounding Glu173. If the partial charges of residues immediately surrounding Glu173 are removed, the exothermicity is -86 kcal/mol, which is very similar to the gas-phase result. The energetic discrepancy between the gas phase and complete QM/MM treatments likely is due to the increased negative charge on glutamate after deprotonation, which would be stabilized by the hydrogen bond originating from Tyr169 and Gln436 as well as by the electrostatic interactions with Arg174 and Lys167. We next considered the feasibility of direct proton transfer from the Thr260 hydroxyl to the distal O2 atom as a route for Cmpd 0 generation. With ξThr-Od defined as the distance between the Thr260OH and distal O2 atom, restraint of the hydroxyl proton to ξThr-Od e 1.35 Å results in initiation of proton transfer from Glu173 via water molecules A-C. Subsequent release of this restraint, followed by full geometry optimization, results in the formation of Cmpd 0, deprotonation of Glu173, and return of the restrained proton to the Thr260 oxygen atom. When ξThr-Od ) 1.4 Å, the Glu173 is not deprotonated. However, this point lies 5.8 (5.5) kcal/mol above the 1AFeO2 peroxo species, slightly higher than when water mediates Cmpd 0 formation. These values are at least a lower bound for the energetic barrier necessary for Thr260 to behave as the ultimate proton donor. Thus, given the small barrier necessary for water alone to mediate the proton transfer in Cmpd 0 formation, these results suggest that under these circumstances the Thr260 hydroxyl functions to maintain the orientation of the O2 ligand for successful protonation of the distal O2 atom. Alternatively, protonation of the proximal O2 would be a first step toward H2O2 production and idle consumption of reducing equivalents by the enzyme (Scheme 1, DfG). This observation is entirely consistent with the experimental results obtained for the Thr252Ala and Thr268Ala mutants of CYP101 and CYP102, respectively. Each suffers from uncoupling of electron transfer from monooxygenation.20,68,69 As demonstrated above, regardless

CYP51 Molecular Oxygen Activation

J. Phys. Chem. B, Vol. 113, No. 23, 2009 8177

Figure 6. B1:CHARMM fully optimized geometries, energies (relative to the peroxo species; kcal/mol), spin densities, and charges (parentheses) along the reaction coordinate for the Glu173-(H2O)4-mediated transformation of the peroxo species to Cmpd 0. B2:CHARMM//B1:CHARMM energies are in parentheses in the order (1AFeO2 and 3AFeO2). Important bond distances are listed in Å.

of the defined reaction coordinate, only three (A, B, and C) of the four water molecules directly mediate proton flow between the Glu173 and distal O2 atom. To dissect any energetic role for water molecule D, this was deleted, and the restrained optimizations were performed in its absence (see Figure 6). When water is considered to be the ultimate proton donor, the approximate transition state (ξThr-Od ) 1.45 Å) lies (2.4 (1.5) kcal/mol) above the corresponding peroxo species. This result supports that the role of the fourth water molecule is primarily structural, maintaining the water network connecting Glu173 to the distal O2 atom. Following evaluation of Glu173 as the proton source catalyzing the peroxo to Cmpd 0 transformation, several mechanisms for the conversion of Cmpd 0 to Cmpd I were considered. First, mechanisms in which the Thr260 hydroxyl and adjacent water molecule serve as the terminal donors to catalyze heterolytic scission of the O-O bond were explored. These were followed by a consideration of mechanisms in which Cmpd I formation is initiated by homolytic scission of the Op-Od bond. In the latter mechanisms, roles of the aforementioned proton donors, as well as that of the Glu173, are discussed. Initially considering the terminal water molecule (C, Figure 6) as the proton donor in heterolytic cleavage of the Op-Od bond, a series of restrained geometry optimizations along ξH2O-Od were performed. An approximate transition state (ξH2O-Od ) 1.15 Å) was identified 12.2 (12.5) kcal/mol above the corresponding Cmpd 0 species. The fully optimized geometries for this approximate transition state and the corresponding reactant and product complexes are shown in Figure 7, as well as relevant charges and spin densities. Unlike the peroxo to Cmpd 0 transformation utilizing this pathway, no excess spin density accumulates on the Glu173 side chain. Full optimization of the species with ξH2O-Od > 1.15 Å returned to the Cmpd 0

intermediate. As hypothesized, when ξH2O-Od < 1.15 Å, optimization results in scission of the Op-Od bond, followed by formation of Cmpd I and water. This transformation of Cmpd 0 to Cmpd I is -58.5 (-53.7) kcal/mol exothermic. Shortening of the Thr260OH-H2O hydrogen bond by 0.19 Å near the transition state supports it may stabilize the proton-transfer reaction coordinate. The Op-Od distance is consistent with an early transition state and a very exothermic reaction, changing very little in the transition state relative to the Cmpd 0 reactant. Furthermore, there are negligible increases in the spin densities and charges on these atoms in the transition state. Recalculation of the potential energy surface after deletion of the water molecule D again results in a negligible increase in the energetic barrier (ξH2O-Od ) 1.10 Å) to 14.8 (15.1) kcal/mol. The Thr260 hydroxyl was next considered the ultimate proton donor involved in heterolytic cleavage. Fully optimized geometries of these species with their relevant charges and spin densities are shown in Figure 8. The approximate transition state (ξThr-Od ) 1.10 Å) is 15.1 (15.9) kcal/mol above the corresponding Cmpd 0 species. Full geometry optimization of species with ξThr-Od < 1.1 Å collapses to a complex containing Cmpd I with an exothermicity of -64.1 (-59.1) kcal/relative to Cmpd 0. Compared to when a water molecule is the donor, there are marked differences in the geometry and electronic structure of this transition state. Examination of the Op-Od bond reveals it is essentially broken (2.15 Å), consistent with a relatively later transition state. Furthermore, there is a notable accumulation of charge (-0.54 f -0.80) and spin density (0.00 f -0.25) on the incipient OH moiety derived from the Od atom. These geometric and electronic structure characteristics are more akin to the transition states characterized for homolytic dissociation of the Op-Od bond (Vide infra). Moreover, the transition state is characterized by a short hydrogen bond (1.44 Å) between

8178

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

Sen and Hackett

Figure 7. B1:CHARMM fully optimized geometries, energies (relative to the Cmpd 0 species; kcal/mol), spin densities, and charges (parentheses) along the reaction coordinate for the Glu173-(H2O)4 heterolytic cleavage of the Cmpd 0 Op-Od bond, affording Cmpd I. B2:CHARMM//B1: CHARMM energies are in parentheses. Important bond distances are listed in Å.

Figure 8. B1:CHARMM fully optimized geometries, energies (relative to the Cmpd 0 species; kcal/mol), spin densities, and charges (parentheses) along the reaction coordinate for the Glu173-(H2O)4-Thr260OH heterolytic cleavage of the Cmpd 0 Op-Od bond affording Cmpd I. B2:CHARMM// B1:CHARMM energies are in parentheses. Important bond distances are listed in Å.

the Thr260 oxygen atom and the adjacent water molecule, consistent with the observed concerted reprotonation of this oxygen en route to Cmpd I.

Similar proton-assisted heterolytic Op-Od bond cleavage of Cmpd 0 in CYP101 has been explored with protonated Glu366 and Asp251 as the proton donors. Zheng and co-workers

CYP51 Molecular Oxygen Activation

J. Phys. Chem. B, Vol. 113, No. 23, 2009 8179

Figure 9. B1:CHARMM fully optimized geometries, energies (relative to the Cmpd 0 species; kcal/mol), spin densities, and charges (parentheses) along the reaction coordinate for the Glu173-(H2O)4 homolytic cleavage of the Cmpd 0 Op-Od bond affording Cmpd I. B2:CHARMM//B1: CHARMM energies are in parentheses. Important bond distances are listed in Å.

predicted both pathways require inaccessible energy barriers or intermediates and disqualified their operation.16 Alternatively, Glu173 in concert with an ordered water network, regardless of the ultimate proton donor, CYP51 offers two heterolytic Op-Od cleavage pathways affording exothermic formation of Cmpd I. Using each of the Cmpd 0 geometries identified during exploration of the proton-assisted heterolytic cleavage, mechanisms for Cmpd I formation, resulting from initial homolytic cleavage of the Op-Od bond, were explored. Geometries, approximate transition states, charges, and spin densities for the water and Thr260 surfaces are shown in Figures 9 and 10,

respectively. In each case, the Op-Od bond distance was designated to be the reaction coordinate (ξOp-Od), and a series of restrained geometry optimizations were performed. From the Cmpd 0 complex in which water is the ultimate proton donor (Figure 7), restrained optimization along ξOp-Od results in an energetically shallow region on the potential energy surface in the vicinity of ξOp-Od ∼ 1.9-2.0 Å. A subsequent unconstrained optimization in this region results in an IC of Cmpd I and an OH moiety 12.2(12.6) kcal/mol above Cmpd 0, an energetic spacing less than 1 kcal/mol different from the barrier computed for the analogous heterolytic mechanism (Figure 7). In the IC, the O-O bond is completely broken (ξOp-Od ) 2.0) and

8180

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

Sen and Hackett

Figure 10. B1:CHARMM fully optimized geometries, energies (relative to the Cmpd 0 species; kcal/mol), spin densities, and charges (parentheses) along the reaction coordinate for the Glu173-(H2O)4-Thr260OH homolytic cleavage of the Cmpd 0 Op-Od bond affording Cmpd I. B2:CHARMM// B1:CHARMM energies are in parentheses. Important bond distances are listed in Å.

hydrogen bonds to an active site water molecule. Unlike geometries obtained with truncated models of the CYP active sites70 and QM/MM studies of the analogous mechanism in CYP101,16 the OH moiety does not “somersault” resulting in a FeO · · · HO hydrogen bond. Despite the lack of the somersault geometry, the OH moiety adopts an electronic structure similar to that characterized in CYP101, in which the spin density supports the notion that the OH moiety is not a pure OH radical. The corresponding charge dictates this species may react more like a hydroxide anion. Characteristic of this species, the OH moiety undergoes facile reprotonation with a 3.0 (3.4) kcal/

mol barrier by the hydrogen-bonded water molecule in concert with reprotonation by Glu173 with an exothermicity of -45.4 (-40.7) kcal/mol. Finally, a mechanism initiated by homolytic Op-Od bond cleavage, followed by protonation of the incipient OH moiety by the Thr260OH, was considered (Figure 10). The IC is 14.3 (14.6) kcal/mol above the corresponding Cmpd 0 species and has similar electronic structure to the aforementioned IC. Of notable interest in this case is the geometry of the approximate transition state for facile reprotonation of the OH moiety by Thr260. Geometric parameters including ξOp-Od, ξThr-Od, and the

CYP51 Molecular Oxygen Activation length of the strong hydrogen bond between the Thr260 hydroxyl and the adjacent water molecule are very similar to those in the approximate transition state for Thr260OH-mediated heterolytic O-O cleavage (Figure 8). As discussed earlier, the size of the QM regions described herein has prevented our ability to locate authentic transition states for these reactions. However, the geometric similarities in the vicinity of bond making and bond breaking, as well as their essentially identical energetic spacing relative to Cmpd 0, lead us to speculate that the Thr260OH-mediated “heterolytic” and “homolytic” reaction mechanisms may encounter a common transition state toward formation of Cmpd I. Each of the reaction paths results in exothermic formation of a triradicaloid Cmpd I species. Population analyses confirm the triplet electronic configuration of the FeO unit and the presence of a third electron of antiparallel spin localized primarily to the porphyrin π system (Figures 7-10).12,71-74 These exothermicities deserve comment in light of the results of Zheng and co-workers’ description of thermoneutral formation of CYP101 Cmpd I via a mixed homolytic/heterolytic mechanism. An important difference between the electrostatic environment surrounding Glu173 as compared to Glu366 in CYP101, for example, is that the former residue is surrounded by several positively charged amino acids, including Arg174 and His259. Upon deprotonation of Glu173, it is likely the resulting anion develops favorable electrostatic interactions with the surrounding environment contributing to the exothermicity. Dissecting the electrostatic contributions of individual amino acids to the electronic structure of CYP51 catalytic intermediates and their relative energetics is an area of future study. Conclusions We present the application of extended, implicit solvent MD simulations and hybrid QM/MM techniques to elucidate the mechanism of O2 activation and characterize reactive oxygen intermediates in this important enzyme. Simulations of His259H+ Mtb CYP51 support evolution of two distinct hydrogen-bond networks with the propensity to shuttle protons to the peroxo species in the CYP51 catalytic cycle. The first network is characterized by hydrogen-bonding between the His259, Thr260OH, and distal O2 atom. Hybrid QM/MM calculations support this hydrogen-bond network can catalyze the formation of Cmpd 0 with a slightly higher barrier and comparable exothermicity to that previously described in QM/ MM studies of the Asp251-H2O-Thr252 shuttle of CYP101.15 However, pursuit of a heterolytic O-O cleavage mechanism for the subsequent formation of Cmpd I was unsuccessful, and exploration of a mechanism initiated by Op-Od bond homolysis realized an endothermic reaction. Disruption of the His259H+-Thr260OH hydrogen bond in MD simulations was followed by the influx of water into the active site and the evolution of an apparent second protontransfer network, connecting the distal O2 atom to His259H+ and Glu173 via four water molecules. In this configuration of the active site, QM/MM calculations revealed the peroxo intermediate has an unprecedented triradicaloid electronic structure with either two parallel or antiparallel electrons localized to the FeO2 unit, while a third resides on the protonated Glu173 side chain. Further calculations demonstrated this electronic structure is a direct result of the local hydrogen bond and electrostatic environment contributed by the enzyme interior, illustrating an important role for the “protein environment” to tune the electronic structure of the peroxo intermediate. The computational results reported herein support that protonation

J. Phys. Chem. B, Vol. 113, No. 23, 2009 8181

Figure 11. Positions of the I-helices, F-helices, and structurally conserved acidic residues of CYP101 (PDB 1DZ8, magenta), CYP119 (PDB 1IO7, yellow), and CYP107A1 (PDB 1Z8O, blue) relative to the Mtb CYP51 structure (PDB 1X8V, green).

of this residue gives rise to a peroxo intermediate with an electronic structure reflecting a potentially less reactive oxyferrous species in the catalytic core. Hence, protonation of this distant residue (∼9 Å) can tune the reactivity of this intermediate. Unlike hurdles encountered with observation of Cmpd I in CYP enzymes, characterization of the peroxo intermediate using EPR, UV-vis, and resonance Raman spectroscopies of radiolytically reduced oxyferrous CYP enzymes at cryogenic temperatures have been successful.10,51,55-57 It will be of great interest to the field to validate the results of these calculations using similar experimental means as well as to probe the sensitivity of the triradicaloid character to experimental conditions. Since this catalytic intermediate has been implicated in the final C-C bond cleavage step catalyzed by CYP51 and other members75-77 of this enzyme family, this is an intriguing observation. Owing to the low barrier for protonation of the peroxo intermediate, it is unlikely this intermediate accumulates during catalytic turnover. For the peroxo intermediate to participate in the final step, nucleophilic addition of the aldehyde intermediate (Scheme 2) must proceed with a negligible barrier or this species may disrupt the proton shuttle, extending the peroxo lifetime for participation in substrate oxidation. Glu173 is located in the CYP51 F-helix adjacent to other polar residues with apparent access to bulk water, and in this milieu, its transient protonation is a likely scenario. These calculations were unsuccessful in mapping the reaction for protonation of the peroxo species by the Thr260 hydroxyl; however, they do support that active-site water may take this role. In Glu173-(H2O)4 and Glu173-(H2O)4-Thr260OH, proton shuttles are capable of catalyzing the transformation of Cmpd 0 to Cmpd I. Each provides this enzyme with a repertoire of concerted, energetically competitive “heterolytic” and “homolytic” reaction pathways utilizing both active site water and the Thr260 hydroxyl as the proton source for exothermic O2bond scission. Comparison of CYP51 with the experimental structures of other CYPs reveals an acidic residue is structurally conserved in the F-helices of many of these enzymes. The structural conservation of this residue in CYP101,9 CYP107A1,78 CYP119,79 and CYP5129 is shown in Figure 11, indicating similar proton transport mechanisms may have evolved in other members of this important enzyme family. Acknowledgment. The authors are grateful to the Ohio Supercomputer Center for a generous allocation of computational resources.

8182

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

Supporting Information Available: Cartesian coordinates of atoms in the QM layer, absolute energies, complete references, and additional computational data discussed herein. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Sono, M.; Roach, M. P.; Coulter, E. D.; Dawson, J. H. Chem. ReV. 1996, 96, 2841–2888. (2) Meunier, B.; de Visser, S. P.; Shaik, S. Chem. ReV. 2004, 104, 3947. (3) Denisov, I. G.; Makris, T. M.; Sligar, S. G.; Schlichting, I. Chem. ReV. 2005, 105, 2253–2277. (4) Guengerich, F. P. Curr. Drug Metab. 2001, 2, 93–115. (5) In Cytochrome P450: Structure, Mechanism, and Biochemistry; Ortiz de Montellano, P. R., Ed.; Kluwer Academic/Plenum Publishers: New York, U.S.A., 2005. (6) Rendic, S. Drug Metab. ReV. 2002, 34, 83. (7) Miller, W. L. Endocr. DeV. 2008, 13, 1. (8) Brueggemeier, R. W.; Hackett, J. C.; Diaz-Cruz, E. S. Endocr. ReV. 2005, 26, 331. (9) Schlichting, I.; Berendzen, J.; Chu, K.; Stock, A. M.; Maves, S. A.; Benson, D. E.; Sweet, R. M.; Ringe, D.; Petsko, G. A.; Sligar, S. G. Science 2000, 287, 1615. (10) Denisov, I. G.; Mak, P. J.; Makris, T. M.; Sligar, S. G.; Kincaid, J. R. J. Phys. Chem. A 2008, 112, 13172. (11) Vidakovic, M.; Sligar, S. G.; Li, H.; Poulos, T. L. Biochemistry 1998, 37, 9211. (12) Shaik, S.; Kumar, D.; de Visser, S. P.; Altun, A.; Thiel, W. Chem. ReV. 2005, 105, 2279. (13) Guallar, V.; Friesner, R. A. J. Am. Chem. Soc. 2004, 126, 8501. (14) Kumar, D.; Hirao, H.; de Visser, S. P.; Zheng, J.; Wang, D.; Thiel, W.; Shaik, S. J. Phys. Chem. B 2005, 109, 19946. (15) Wang, D.; Zheng, J.; Shaik, S.; Thiel, W. J. Phys. Chem. B 2008, 112, 5126. (16) Zheng, J.; Wang, D.; Thiel, W.; Shaik, S. J. Am. Chem. Soc. 2006, 128, 13204. (17) Zheng, J.; Altun, A.; Thiel, W. J. Comput. Chem. 2007, 28, 2147. (18) Kamachi, T.; Yoshizawa, K. J. Am. Chem. Soc. 2003, 125, 4652. (19) Nebert, D. W.; Nelson, D. R.; Coon, M. J.; Estabrook, R. W.; Feyereisen, R.; Fujii-Kuriyama, Y.; Gonzalez, F. J.; Guengerich, F. P.; Gunsalus, I. C.; Johnson, E. F. DNA Cell Biol. 1991, 10, 1. (20) Imai, M.; Shimada, H.; Watanabe, Y.; Matsushima-Hibiya, Y.; Makino, R.; Koga, H.; Horiuchi, T.; Ishimura, Y. Proc. Natl. Acad. Sci. U.S.A. 1989, 86, 7823. (21) Kimata, Y.; Shimada, H.; Hirose, T.; Ishimura, Y. Biochem. Biophys. Res. Commun. 1995, 208, 96. (22) Lepesheva, G. I.; Waterman, M. R. Mol. Cell. Endocrinol. 2004, 215, 165. (23) Akhtar, M.; Alexander, K.; Boar, R. B.; McGhie, J. F.; Barton, D. H. Biochem. J. 1978, 169, 449. (24) Shyadehi, A. Z.; Lamb, D. C.; Kelly, S. L.; Kelly, D. E.; Schunck, W. H.; Wright, J. N.; Corina, D.; Akhtar, M. J. Biol. Chem. 1996, 271, 12445. (25) Groves, J. T.; McClusky, G. A. J. Am. Chem. Soc. 1976, 98, 859. (26) Groves, J. T.; McClusky, G. A. Biochem. Biophys. Res. Commun. 1978, 81, 154. (27) Hackett, J. C.; Brueggemeier, R. W.; Hadad, C. M. J. Am. Chem. Soc. 2005, 127, 5224. (28) Podust, L. M.; Poulos, T. L.; Waterman, M. R. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 3068. (29) Podust, L. M.; Yermalitskaya, L. V.; Lepesheva, G. I.; Podust, V. N.; Dalmasso, E. A.; Waterman, M. R. Structure 2004, 12, 1937. (30) Dolinsky, T. J.; Czodrowski, P.; Li, H.; Nielsen, J. E.; Jensen, J. H.; Klebe, G.; Baker, N. A. Nucleic Acids Res. 2007, 35, W522. (31) Kieseritzky, G.; Knapp, E. W. J. Comput. Chem. 2008, 29, 2575. (32) Bas, D. C.; Rogers, D. M.; Jensen, J. H. Proteins: Struct., Funct., Bioinf. 2008, 73, 765. (33) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. J. Comput. Chem. 2005, 26, 1781. (34) MacKerell, A. D.; et al. J. Phys. Chem. B 1998, 102, 3586. (35) Cournia, Z.; Smith, J. C.; Ullmann, G. M. J. Comput. Chem. 2005, 26, 1383.

Sen and Hackett (36) Zhang, Y.; Feller, S. E.; Brooks, B. R.; Pastor, R. W. J. Chem. Phys. 1995, 103, 10252. (37) Martyna, G. J.; Tobias, D. J.; Klein, M. L. J. Chem. Phys. 1994, 101, 4177. (38) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (39) Sherwood, P.; et al. THEOCHEM 2003, 632, 1. (40) Ahlrichs, R.; Ba¨r, M.; Ha¨ser, M.; Horn, H.; Ko¨lmel, C. Chem. Phys. Lett. 1989, 162, 165. (41) Smith, W.; Forester, T. R. J. Mol. Graphics 1996, 14, 136. (42) Bakowies, D.; Thiel, W. J. Phys. Chem. 1996, 100, 10580. (43) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (44) Becke, A. D. Phys. ReV. A: Math. Gen. 1988, 38, 3098. (45) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B: Condens. Matter Mater. Phys. 1988, 37, 785. (46) Scha¨fer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 2571. (47) Wachters, A. J. H. J. Chem. Phys. 1970, 52, 1033. (48) Bauschlicher, C. W., Jr.; Langhoff, S. R.; Partridge, H.; Barnes, L. A. J. Chem. Phys. 1989, 91, 2399. (49) Scha¨fer, A.; Huber, C.; Ahlrichs, R. J. Chem. Phys. 1994, 100, 5829. (50) Reed, A. E.; Weinstock, R. B.; Weinhold, F. J. Chem. Phys. 1985, 83, 735. (51) Davydov, R.; Kappl, R.; Hu¨ttermann, J.; Peterson, J. A. FEBS Lett. 1991, 295, 113. (52) Harris, D.; Loew, G.; Waskell, L. J. Am. Chem. Soc. 1998, 120, 4308. (53) Wang, D.; Thiel, W. J. Mol. Struct. 2009, 898, 90. (54) Harris, D. L.; Loew, G. H. J. Am. Chem. Soc. 1994, 116, 11671. (55) Sjodin, T.; Christian, J. F.; Macdonald, I. D.; Davydov, R.; Unno, M.; Sligar, S. G.; Hoffman, B. M.; Champion, P. M. Biochemistry 2001, 40, 6852. (56) Davydov, R.; Makris, T. M.; Kofman, V.; Werst, D. E.; Sligar, S. G.; Hoffman, B. M. J. Am. Chem. Soc. 2001, 123, 1403. (57) Denisov, I. G.; Makris, T. M.; Sligar, S. G. J. Biol. Chem. 2001, 276, 11648. (58) Sligar, S. G.; Makris, T. M.; Denisov, I. G. Biochem. Biophys. Res. Commun. 2005, 338, 346. (59) Harris, D. L.; Loew, G. H. J. Am. Chem. Soc. 1998, 120, 8941. (60) Becke, A. D. J. Chem. Phys. 1993, 98, 1372. (61) Dirac., P. A. M. Proc. R. Soc. London 1929, 123, 714. (62) Slater, J. C. Phys. ReV. 1951, 81, 385. (63) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B: Condens. Matter Mater. Phys. 1988, 37, 785. (64) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (65) Perdew, J. P.; Wang, Y. Phys. ReV. B: Condens. Matter Mater. Phys. 1992, 45, 13244. (66) Perdew, J. P.; Ernzerhof, M.; Burke, K. J. Chem. Phys. 1996, 105, 9982. (67) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Phys. ReV. Lett. 2003, 91, 146401. (68) Yeom, H.; Sligar, S. G. Arch. Biochem. Biophys. 1997, 337, 209. (69) Martinis, S. A.; Atkins, W. M.; Stayton, P. S.; Sligar, S. G. J. Am. Chem. Soc. 1989, 111, 9252. (70) Bach, R. D.; Dmitrenko, O. J. Am. Chem. Soc. 2006, 128, 1474– 1488. (71) Ogliaro, F.; Cohen, S.; Filatov, M.; Harris, N.; Shaik, S. Angew. Chem., Int. Ed. 2001, 40, 647. (72) Ogliaro, F.; de Visser, S. P.; Cohen, S.; Kaneti, J.; Shaik, S. ChemBioChem 2001, 2, 848. (73) Ogliaro, F.; de Visser, S. P.; Groves, J. T.; Shaik, S. Angew. Chem., Int. Ed. 2001, 40, 3503. (74) Schoneboom, J. C.; Neese, F.; Thiel, W. J. Am. Chem. Soc. 2005, 127, 5840. (75) Akhtar, M.; Calder, M. R.; Corina, D. L.; Wright, J. N. Biochem. J. 1982, 201, 569. (76) Akhtar, M.; Corina, D.; Miller, S.; Shyadehi, A. Z.; Wright, J. N. Biochemistry 1994, 33, 4410. (77) Lambeth, J. D.; Kitchen, S. E.; Farooqui, A. A.; Tuckey, R.; Kamin, H. J. Biol. Chem. 1982, 257, 1876. (78) Nagano, S.; Cupp-Vickery, J. R.; Poulos, T. L. J. Biol. Chem. 2005, 280, 22102. (79) Park, S. Y.; Yamane, K.; Adachi, S.; Shiro, Y.; Weiss, K. E.; Maves, S. A.; Sligar, S. G. J. Inorg. Biochem. 2002, 91, 491.

JP902932P