Computational Study on Compound I Redox-Active Species in

May 3, 2010 - Materiali, UniVersita di L'Aquila, Via Vetoio 67100, L'Aquila, Italy. ReceiVed: February 3, 2010. Molecular dynamics simulations for the...
0 downloads 0 Views 3MB Size
Download PDF
J. Phys. Chem. B 2010, 114, 6817–6824

6817

Computational Study on Compound I Redox-Active Species in Horseradish Peroxydase Enzyme: Conformational Fluctuations and Solvation Effects Costantino Zazza,*,† Amedeo Palma,*,‡ Nico Sanna,† Simone Tatoli,†,§ and Massimiliano Aschi⊥ CASPUR, Consorzio InteruniVersitario per le Applicazioni di Supercalcolo per UniVersita` e Ricerca, Via dei Tizii, 6/b, 00185 Roma, Italy, Istituto per lo Studio dei Materiali Nanostrutturati, CNR-ISMN, Via Salaria Km. 29.3, Sez. Montelibretti, Monterotondo S.(RM), Italy, Dipartimento di Chimica, UniVersita` di Roma La Sapienza, P. le A. Moro 00185, Rome, Italy, and Dipartimento di Chimica, Ingegneria Chimica e Materiali, UniVersita di L’Aquila, Via Vetoio 67100, L’Aquila, Italy ReceiVed: February 3, 2010

Molecular dynamics simulations for the compound I species in horseradish peroxidase were carried out over a nanoseconds time-scale. Results indicate that the supramolecular assembly composed of compound I in interaction with highly conserved distal residues (His42 and Arg38) exists in two well-defined conformations basically differing in the local position of the distal histidine (i.e., His42). Furthermore, we observe the presence of a biological channel in the distal side of the heme cavity, between Arg38 and Pro139 residues, that represents a direct link connecting the compound I (CpdI) species to bulk molecules. Our investigation supports the idea that when CpdI is formed, the biological machinery relaxes the local electrostatic forces, opening a structural channel through which an exchange of water molecules with the bulk solvent takes place without any significant kinetic barrier. Interestingly, we also show that the combined effect of enzyme and solvent, modulated by thermal fluctuations, affects the order and the energy difference between the lowest doublet and quartet magnetic states of the CpdI-His42-Arg38 complex. Consequently, when passing from the gas phase to the biological environment, the doublet spin state becomes slightly more stable than the higher spin multiplicity state. The distribution of the perturbed low-lying states energy variation, induced by surrounding fluctuations, yields results rather close to that obtained by other state of the art quantum-mechanics/molecular mechanics calculations, and still in line with previous polarizable continuum calculations. 1. Introduction Horseradish peroxidase (HRP) is a monomeric heme containing enzyme (see Figure 1) that efficiently removes toxic hydrogen peroxide from the intracellular region in mammalian tissues and plants.1 HRP catalyzes indirectly the oxidation of a large variety of substrates in nature by means of the stepwise Poulos-Kraut (P-K) mechanism,2 through which starting from a hydroperoxoferric complex (called peroxy complex and characterized by the hydrogen peroxide bound to ironIII of the heme group) produces a ferric oxo radical cation (CpdI) and a newly formed water molecule. The process rapidly eliminates hydrogen peroxide, forming an inorganic species having the iron with a high oxidation state (FeIVdO) after a formal intramolecular charge transfer from the FeIIIO moiety of a ferric hydroperoxide anion intermediate compound 0 (Cpd0) to the porphyrin ring of the heme.2 As a final step, the hosting matrix plays a key role in solvating the compound I species that is known to have a strong oxidizing character;3 indeed, one of the main biological function performed by the HRP/CpdI system is connected to the oxidation of natural phenols, which subsequently polymerize, creating a natural material able to build and repair the cell walls in plants.1 * Corresponding authors. (C.Z.) Phone: +39 06 44486720. Fax: +39 06 4957083. E-mail: [email protected]. (A.P.) Phone: +39 06 90672343. Fax: +39 06 90672839. E-mail: [email protected]. † CASPUR. ‡ CNR-ISMN. § Universita di Roma La Sapienza. ⊥ Universita di L’Aquila.

Figure 1. The horseradish peroxidase (HRP) enzyme (right panel) and its own heme-containing active site (left panel).

Several experimental4 and theoretical5 investigations have identified such a redox-active species in HRP enzyme, revealing its own open shell electronic structure. In particular, it is clearly reported in the literature that CpdI is electrostatically stabilized by highly conserved distal residues such as His42 and Arg38.5d Furthermore, the explicit inclusion of distal residues in the theoretical modeling of the reaction pocket is found to be necessary to reproduce recent X-ray data regarding the oxo-iron distance,4i as well as Mo¨ssbauer parameters supporting the assignment of a spin coupling between the FeIVdO triplet moiety and a porphyrin cation radical.4f In this context, quantum mechanical studies regarding the CpdI-His42-Arg38 complex in HRP have confirmed the

10.1021/jp101033w  2010 American Chemical Society Published on Web 05/03/2010

6818

J. Phys. Chem. B, Vol. 114, No. 20, 2010

aforementioned open shell character both using density functional theory (DFT)5a-e and a wave-function-based5f description. In turn, from these studies, a further interesting question is emerged concerning the spin multiplicity of the electronic ground state of the CpdI species in HRP enzyme. Quantum mechanics/molecular mechanics (QM/MM) calculations, in which the CpdI-His42-Arg38 supramolecular assembly was treated in the framework of DFT, predicted almost degenerate doublet and quartet spin magnetic states basically differing for the alpha-beta occupation of the unpaired electron localized over the porphyrin group.5b,d On the other hand, a wave-functionbased approach using a quadratic configuration interaction procedure with single and double electronic excitations (QCISD)6 in conjunction with the conductor-like polarizable continuum method (CPCM),7 revealed a doublet spin magnetic state at lower energy (-3.61 kcal/mol) with respect to a quartet one.5f However, in light of recent results,5i,8 the idea that the electronic state of a portion of a complex system (e.g., the spin magnetic state of a prosthetic group within a solvated protein) might be heavily influenced not only by the chemical, but also by the mechanical and dynamical features of the perturbing surrounding environment is becoming widely accepted. Therefore, it seems rather interesting to investigate the conformational space sampled by the [HRP-CpdI] enzyme and its influence in modulating the energy difference between the lowest spin multiplicities of CpdI species. Present results show that the Cpd0 f CpdI covalent transition in HRP is accompanied by a conformational relaxation step which irreversibly opens in the HRP[CpdI] simulated ensemble a structural gate directly connecting the active species with the bulk solvent. It is worth noting that this finding is in line with photothermal studies of carbon monoxide photodissociation from HRP enzyme in which the fast release of CO from the distal side of the heme pocket would be attributed to a direct channel linking the heme to the bulk.9 Inspired by these open issues and on the basis of our previous studies on wild-type5i and R38L mutant in HRP,10 we decided to study the CpdI species in aqueous HRP by an extended nanoseconds time-scale molecular dynamics (MD) simulation with the specific aim of (i) characterizing the mechanical and dynamical behavior of HRP matrix upon the formation of the CpdI according to P-K mechanism2 and (ii) providing some additional indication on the actual magnetic state of CpdI in HRP and under physiological conditions (i.e., dilute water solution, room temperature). 2. Computational Details 2.1. Classical Molecular Dynamics Simulation Setup. To study CpdI active species within the HRP framework, MD simulations were initiated using as starting coordinates the latest configuration sampled at 298 K for the Cpd0 intermediate in the HRP enzyme.5i In turn, the covalent contacts of Cpd0 were interconverted into a CpdI-H2O complex postulated by the P-K mechanism.2 The same rectangular box of 7.45 × 8.72 × 6.97 nm3 dimensions filled with single point charge11 (SPC) water molecules at a typical density of 1000 kg/m3 was used for the present simulation. As already reported in our recent works, structural Ca2+ ions were also included in the simulated system because of their well-known importance in maintaining the structural integrity of the heme group.12 A proper number of counterions (1 Cl-) was added to ensure the overall electrical neutrality of the system. The standard Gromos96 force field13 was adopted in our simulations. Atomic point charges were

Zazza et al. recalculated over the CpdI active species by an electrostatic potential (ESP) fitting procedure14 in conjunction with density functional theory calculations15 using the hybrid Becke’s three parameters exchange and Lee-Yang-Parr correlation functional (B3LYP)15,16 in conjunction with the following atomic basis sets: (i) for iron, the LANL2DZ effective core potential (ECP) for the inner electrons and a double zeta Gaussian basis set of (5s, 5p, 5d)/[3s, 3p, 2d] quality for valence electrons;17 (ii) the standard 6-311G(d) Gaussian basis set for nitrogen and oxygen; and (iii) the 3-21G Gaussian basis set for carbon and hydrogen.18 After an initial optimization using the steepest-descent algorithm, the system was gradually heated from 50 to 300 K using short (20 ps) molecular dynamics runs. The trajectories were then propagated up to 20 ns in a canonical (i.e., NVT) ensemble applying an integration step of 2.0 fs with the rototranslational constraint19 applied to the HRP-[CpdI] system. The temperature was kept constant by isokinetic temperature coupling,20 and all bond lengths were constrained using the standard LINCS algorithm.21 Long range electrostatics was computed by the particle mesh Ewald method22 with 34 wave vectors in each dimension and a fourth-order cubic interpolation. The obtained MD trajectory was then used to investigate the effects of HRP conformational transitions and solvent dynamics on the relative stability of doublet and quartet spin multiplicity states of CpdI-His42-Arg38 complex. To this purpose, an approximated approach was adopted on the basis of the perturbation theory.5i Briefly, the perturbing electrostatic field E and potential V exerted by the HRP enzyme framework and solvent was evaluated at each step of the classical sampling, on the center of mass (termed as rcom) of CpdI-His42-Arg38 subsystem and was used to evaluate the perturbed electronic energy pert by the following equation,

pert ) 0 + qTV[rcom, x(t)] - E[rcom, x(t)]µ

(1)

where 0 is the ground state electronic energy previously studied at the QCISD6 level of theory, according to our recent theoretical/computational investigations;5f qT is the total charge of the active site, µ is its permanent electric dipole moment, and x(t) represents the atomic coordinates of the perturbing environment. The previous procedure carried out for both the spin states has allowed an estimation of the perturbed electronic energy difference between doublet and quartet spin magnetic states along the MD simulations at 298 K. All the MD simulations were carried out using a modified version of the Gromacs software23 using the rototranslational constraint to the whole solute,19 while quantum chemical calculations were performed using the Gaussian 03 package.24 2.2. Analysis Tools. The conformational transitions of the C-R skeleton in the enzyme matrix were characterized by means of the essential dynamics (ED) technique. Details of the method are described elsewhere;25 here, we report only some basic features. Briefly, the covariance matrix of the atomic positional fluctuations of the entire enzyme and of a subportion of it is built and diagonalized, providing a set of generalized atomic coordinates (i.e., the eigenvectors) along which atomic motions occur. The eigenvectors with the largest associated eigenvalues form the essential subspace, which characterizes the largest amplitude (concerted) conformational fluctuations.25 This procedure was adopted for investigating both the conformational transitions of the overall enzyme using C-R covariance matrix and the fluctuations of CpdI plus the related distal residues (i.e., Arg38, His42, Pro139) using an all-atom covariance matrix. In

CpdI Redox-Active Species in HRP Enzyme

J. Phys. Chem. B, Vol. 114, No. 20, 2010 6819 cases, the position of the solvent molecules is then analyzed by computing volumetric density maps, which represents wet molecular isosurfaces with a resolution of 0.05 nm.

Figure 2. Components of the first (upper panel) and second (lower panel) essential eigenvectors according to essential dynamics technique carried out on the C-R of the HRP enzyme in the presence of CpdI (black) and Cpd0 intermediate (red), respectively.

this latter case, to provide some evidence about the alteration of the fluctuation patterns upon Cpd0 f CpdI transition, a direct comparison between statistically relevant conformations in the two simulated ensemble (i.e., HRP[Cpd0] and HRP[CpdI]) is also analyzed and discussed. Time dependency of the distance between Arg38 and Pro139 side chain geometrical centers was also monitored along the simulation to plausibly define a structural gate responsible for the accessibility of solvent water molecules into the distal side of the heme cavity (as described in the next section). In this context, to better characterize the distal hydration shell (if any) and the exact position of the solvent channel in HRP enzyme, we have extended an algorithm already applied in our previous works.5i Here, we discriminate between distal and proximal regions of the active site by projecting the center of mass of each water molecule present into a spherical cutoff of 0.5 nm centered on the iron atom, onto a unit vector defined, at each MD frame, parallel to the vector interconnecting the iron atom and the nitrogen atom of the proximal histidine (i.e., His170). The spherical cutoff was then extended up to 1.2 nm to follow the exit pathways of water molecules from the previously defined distal region. In both

3. Results and Discussion 3.1. Mechanical Behavior of HRP[CpdI] System in Water at 298 K. As an initial step, we addressed the question as to whether and to what extent a change in the covalent framework of the active site, from Cpd0 to CpdI, may affect the structural and mechanical behavior of HRP hosting matrix. For this purpose, in analogy with previous analysis carried out on the R38L mutant of HRP enzyme,10 we have compared the largest amplitude fluctuations of the C-R atoms of the protein, as extracted by means of standard essential dynamics (ED) technique25 carried out considering the two different simulated ensembles (i.e., HRP[CpdI] and HRP[Cpd0]). First of all, we observe that, if compared to the analogous analysis for HRP[Cpd0], the overall enzyme slightly reduces its conformational flexibility upon formation of CpdI in the active site. The first and second essential eigenvectors of the HRP[CpdI], showing related eigenvalues of 0.6 nm2 and 0.2 nm2 as well as the whole trace of the covariance matrix (2.1 nm2) appear slightly lower than the corresponding values previously obtained for HRP[Cpd0] of 0.9, 0.3, and 2.6 nm2. However, the most interesting result is undoubtedly the change in the nature of the fluctuation patterns accompanying the Cpd0-CpdI transition. As a matter of fact, both the atomic composition of the first essential eigenvectors of the HRP[CpdI] and HRP[Cpd0] shown in Figure 2 and pictorially evidenced in Figure 3 as well as the scarce overlap of the projection of the two trajectories onto the essential plane defined by the two largest eigenvectors of the HRP[CpdI] (see Figure 4) witnesses the strong difference in the enzyme fluctuation patterns. 3.2. Local Rearrangements Induced by CpdI Formation in HRP. The similar ED-based analysis, carried out to highlight alterations in the local conformational rearrangements of CpdI with respect to Cpd0, evidences the presence of two different conformational states in the heme pocket. These states have been characterized by projecting the CpdI-His42-Arg38 trajectory, extracted from the overall system, onto the plane defined by the first two eigenvectors obtained by the diagonalization of the CpdI-His42-Arg38 all-atom covariance matrix. The result, reported in Figure 5 shows two distinct spots corresponding to two different conformational basins, termed as conf1 and conf2,

Figure 3. View of the first two essential eigenvectors motion of the C-R atoms of the HRP enzyme structure; the first two essential eigenvectors have been represented superimposing different molecular structures obtained projecting the trajectory onto each eigenvector.

6820

J. Phys. Chem. B, Vol. 114, No. 20, 2010

Figure 4. Projection of the C-R trajectories of the HRP enzyme onto the essential plane defined by the first two essential eigenvectors of the HRP[CpdI] system: the trajectory of the HRP[CpdI] system is reported in black; that of the HRP[Cpd0], in red.

Zazza et al.

Figure 7. Normalized probability density of the geometrical center distance between Arg38 and Pro139 amino acid residues in HRP[Cpd0] (solid line) and HRP[CpdI] (dashed line) as extracted from the corresponding MD trajectories at 298 K.

reported in Figure 6. Interestingly, the same fluctuating distal residue was responsible for the conformational changes in the Cpd0 species embedded in the framework of the HRP enzyme, although in that case, the His42 basically was found to undergo a torquelike motion.5i The conf1-conf2 relative stability was then evaluated by calculating the Helmholtz free energy using the standard relation:

∆Aconf1fconf2 ) -KbT ln χp

Figure 5. Projection of the all-atom trajectory of the CpdI-Arg38His42 supramolecular assembly onto the essential plane defined by the diagonalization of the corresponding covariance matrix. The blue basin identifies conformation 1; the black basin, conformation 2. In the same panel, we report the free energy difference between them as obtained by eq 2.

basically differing for the mutual orientation of the His42 side chain with respect to the Fe(IV)dO moiety. In other words, the CpdI species in HRP at 298 K rapidly interconverts between conf1 and conf2 simply changing the orientation of the distal histidine around the rotation axis along the C-C bond, as

(2)

in which Kb is the Boltzmann constant and χp represents the probability ratio of the two local conformations within the aforementioned local basins. In the simulated conditions at 298 K, a ∆Aconf1fconf2 value equal to 1.25 kcal/mol is calculated, thus suggesting a spontaneous character of the reversible conf2 f conf1 conformational transition in HRP enzyme, that is, the distal histidine preferentially interacts with the Fe(IV)dO moiety. We also analyzed the mutual flexibility and orientation of Arg38 and Pro139 side chains, and the results were compared with those previously obtained in the case of Cpd0 species in HRP.5i Results in Figure 7 indicate that the Cpd0 f CpdI transition is accompanied by a sharp elongation of the Arg38-Pro139 mutual distance (maximum of the distributions at 0.56 nm for HRP[Cpd0] and 0.77 nm for HRP[CpdI]). Interestingly, the two curves feature the same broadening (i.e., full-width at half maximum), indicating that, although at a larger distance, these distal residues conserve the mutual flexibility

Figure 6. CpdI-Arg38-His42 system. Three MD frames representing the local motion describing the pathway from conf1 to conf2 along the first essential eigenvector.

CpdI Redox-Active Species in HRP Enzyme

J. Phys. Chem. B, Vol. 114, No. 20, 2010 6821

Figure 9. Normalized probability density of the hydrogen bonding contacts between Arg38 (H2OsH-N), Pro139 (HOHsOdC) and CpdI (HOHsOdFe(IV)) and SPC water molecules found in the distal side of the heme pocket according to MD sampling at 298 K.

Figure 8. 3D volumetric density map (in light gray) of the solvent molecules having a center of mass at 0.5 nm from the heme iron atom; the isosurface has a threshold of 0.1 water molecule/nm3. The most statistically relevant positions are also shown for both the proximal and distal wet volumes by means of the water molecule itself.

when CpdI is formed from Cpd0. Inspired by recent results on R38L mutant,10 we plausibly related the above result to the hydration extent of the internal cavity. As a matter of fact, MD simulations have shown that the replacement of distal Arg38 with LEU residue results in a hydration of the reaction pocket. Such a finding has evidenced the key role of the Arg38 polar side chain that, closely interacting with [OOH]- ligand, provides a mechanical barrier toward the entrance of bulk molecules into the distal side.10 Thus, the loss of the local hydrogen bonding network, essentially due to absence of the [OOH]- hook bound to the Fe(III) upon CpdI formation, might result, also in this case, in a larger hydration extent of the distal region and the related formation of a supramolecular channel connecting bulk solvent and the reaction pocket. In fact, the investigation on the presence of water molecules in the distal region of the internal cavity, differently from previous MD samplings for HRP[peroxycomplex] and HRP[Cpd0]5i and somewhat in analogy with what happens to R38L mutant,10 has revealed, as summarized in Figure 8, two wet regions: (i) the first one localized in the proximal side of the HRP enzyme, as also already observed in our recent works,5i and (ii) the second one in the distal side, where solvent molecules mainly interact, by means of a cooperative hydrogen bonding network, with Pro139, Arg38, and the Fe(IV)dO moiety of the CpdI species (see Figure 9). A deeper inspection of the above results reveals a further intriguing aspect regarding the dynamical character of the distal region hydration. From Figure 10, it emerges that HRP enzyme in water solution can be found, with essentially the same probability, having either a wet or a dry distal cavity. In light of this findings, the overall simulated system (i.e., enzyme plus solvent) may be viewed as a sort of complex supramolecular switch modulating the solvation dynamics of the heme cavity. In the wet conformations, we also note, from the same Figure 10, that a single water molecule is found to interact with the

Figure 10. Normalized probability of dry and wet conformations into the distal side of the heme pocket as extracted from a classical MD trajectory. Note that representative conformations of the active site are also shown.

distal residues. On the other hand, conformations with two water molecules are rather unlikely, at least in the simulated ensemble. In conclusion, a rapid exchange between wet and dry conformations of distal pocket results from our simulations. Such a continuous alternation between these two (wet-dry) states suggests that, when CpdI species is formed, the HRP enzyme should irreversibly open a structural gate connecting the active species with bulk molecules. Considering such an intriguing issue, we investigated the presence of a channel linking the CpdI to solvent molecules. To this end, as already described in section 2.2, we defined a second spherical cutoff at 1.2 nm apart from the heme iron atom; this second layer defines a boundary region between the heme pocket and the bulk, where the hypothesized structural gate might be located. Therefore, the water molecules leaving the first layer (at 0.5 nm from the heme iron atom) are then followed, step by step, up to the second layer. The ensemble of the resulting trajectories summarizes the water molecules’ dynamics into and out of the distal side of the heme pocket, as shown in Figure 11. From this figure, where we fully characterized a 3D volumetric density map of the water molecules’ center of mass, we have evidenced the channel through which the solvent molecules can access or leave the internal cavity. It is

6822

J. Phys. Chem. B, Vol. 114, No. 20, 2010

Zazza et al.

Figure 12. Normalized distribution of the perturbed (doublet-quartet) electronic energy difference (in kcal/mol) of the CpdI-His42-Arg38 complex, surrounded by the HRP hosting matrix and solvent molecules, according to current MD/QCISD calculations.

Figure 11. Structural gate (in magenta) allowing solvent water molecules to exit or to reenter into the distal region of the HRP active site along the MD sampling at room temperature. Note that for the sake of clarity, the bulk solvent has been omitted.

worth noting that, as expected on the basis of previous observations (Figure 9), such a structural gate is found to be localized between the Arg38 and Pro139 side chains, and it can be clearly seen in Figure 11 (lower panel). This finding might support a recent photothermal study regarding the carbon monoxide photodissociation from the HRP enzyme, in which the fast release of CO from the distal side of the heme pocket was ascribed to the existence of a channel directly connecting the heme to the bulk.9 Summarizing, analysis of MD simulations strongly suggests that, differently from the case of Cpd0,5i the HRP molecular architecture allows the CpdI species to be in the condition of interacting with external species to be chemically processed. As a matter of fact, the enzyme mechanical fluctuations catalyze the above primary interaction, opening a structural gate connecting the distal side to the reaction partners present in the bulk solution. 3.3. Solvation Effects at 298 K on the CpdI-His42-Arg38 Magnetic State. The final topic of this work is focused on the effects of aqueous-HRP dynamics on the electronic energy difference between the lowest doublet and quartet spin states of the CpdI-His42-Arg38 assembly. In this context, since the

previous molecular geometries applied in recent QM/MM5b,d and QM/C-PCM works5f match with the conf1 conformational basin of the CpdI-His42-Arg38 subsystem (see Figure 5), we have filtered the classical trajectory, extracting the MD frames that fall in conf1 subspace. Then, in analogy with our recent work,5f we described the electronic degrees of freedom of the CpdI-His42-Arg38 complex at the QCISD6 level while treating the remaining part of the simulation box as atomic point charges that change their positions according to MD sampling. Consequently, the surrounding and fluctuating environment is modeled as an external and homogeneous electrostatic field providing, step by step, a perturbing trajectory on the CpdI-His42-Arg38 center of mass.5i Starting from the unperturbed description in the gas phase, the electronic energies of the doubled and quartet spin multiplicity are then recalculated along the MD simulation in the presence of the most stable local conformation (i.e., conf1) embedded in a classical-like HRP/solvent media. In this respect, we would like to clarify that, similarly to previous literature data,5a,b,d,f also in this study, we disregard the thermodynamically less stable conf2 rearrangement shown in Figure 6. The distribution of the electronic energy difference between doublet and quartet multiplicity, reported in Figure 12, clearly shows that the slightly more polar species, that is, the doublet (having a dipole moment |µ| ) 13.4 D vs |µ| ) 13.1 D of the quartet5f), is found to more strongly interact with the biological environment under physiological conditions (HRP matrix plus solvent molecules at 298 K). This behavior resembles the trend observed addressing the catalytic role of structural fluctuations on the Cpd0 formation in the HRP enzyme:26 Between the two possible alternative mechanistic routes considering only doublet spin multiplicity profiles, the fluctuating aqueous enzyme favors the one involving the more polar transition state by substantially lowering the activation free energy by about -3.0 kcal/mol. In the less polar one, the activation free energy was estimated to decrease by almost -0.4 kcal/mol.26 If we compare the energy range spanned by the curve reported in Figure 12 with the energy value calculated considering the same system model at the same level of theory but in a continuum dielectric media mimicking the biological environment (-3.61 kcal/mol),5f we observe that also in this case, the doublet is predicted to be at lower energies. It is worth noting that in the gas phase, the energy order is inverted: the quartet is at lower energy than the doublet state.5f Using the above distribution of perturbed electronic energy differences, we calculated the reversible work (i.e., the Helm-

CpdI Redox-Active Species in HRP Enzyme

J. Phys. Chem. B, Vol. 114, No. 20, 2010 6823

Figure 13. Comparison between unperturbed (ideal gas-phase condition,5f black) and perturbed doublet-quartet free energy profile; previous CPCM/QCISD profile according to ref 5f is reported in blue; current MD/QCISD estimation, in red.

holtz free energy) to switch between doublet and quartet changing the alpha/beta occupation of the unpaired electron widespread over the porphyrin ring in the ensemble of the lower spin component. To this end, we used the basic statistical mechanics relation:5i

∆ADfQ ) -KBT ln 〈exp[-β(Q - D)]〉Do

(3)

where Q - D is difference between quartet (Q) and doublet (D) perturbed electronic energies; the 〈〉°D term indicates that the average is evaluated within the doublet (D) ensemble and that we assume negligible reorganization energy upon D f Q transition. This latter approximation might be justified by considering the high similarity between quartet and doublet structures. The results schematically reported in Figure 13, in line with the energy distribution pattern in Figure 12, indicate that the presence of protein-solvent perturbation provides a different picture with respect to the model system in the gas phase at the QCISD level.5f In fact, the doublet state is 0.60 kcal/mol more stable than the quartet, resulting in a doublet molar ratio about three times larger than quartet molar ratio at the equilibrium. At the same time, as remarked above, CPCM/QCISD5f computations that disregard the solvent and the conformational flexibility of the enzyme (i.e., solvation dynamics) tend to overestimate the D f Q energy variations. Finally, it should also be noted that the D f Q free energy variation, as obtained by the energy difference distribution shown in Figure 12, is very close to recent accurate QM/MM estimations (-0.30 kcal/ mol) in which the same model system was treated at the DFT level.5b 4. Conclusions The current contribution deals with conformational changes and related solvation effects in a relevant biological system such as the horseradish peroxidase enzyme, following the formation of CpdI redox-active species according to the Poulos-Kraut mechanism. By applying the essential dynamics technique on the nanosecond time scale MD simulations under physiological conditions, we report that a change in the covalent framework of the active site, from Cpd0 intermediate to CpdI, affects not only the local structure of the distal site but, interestingly, also the protein hydration and dynamics of the HRP enzyme. In particular, in our simulating scenario upon CpdI formation, the

investigated system undergoes a local conformational relaxation process driven by a complex balance between electrostatic forces and thermal effects, resulting in an active site rather different when compared to that one optimized to host the Cpd0 in the enzyme cavity. Therefore, in the HRP[CpdI] system, the distal histidine (i.e., His42) in the active site, oscillating between two different rotamers, splits up the CpdI-His42-Arg38 model system into two distinct conformational basins. Furthermore, the mechanical response to CpdI formation also involves a rearrangement in the mutual position of the Arg38 and Pro139 side chains, which irreversibly opens a structural channel that connects the active species to bulk water molecules. In this respect, we do not exclude the possibility that this channel, through which water molecules can easily exit or reenter the active site, may also be the biological gate responsible for the access of molecular substrates to allow their interaction with CpdI species in nature. Here, we also report that solvation effects, as arising from protein and solvent dynamics, play a key role in fine-tuning the spin multiplicity of an active site model system. In fact, as far as the energy difference between the lowest doublet and quartet magnetic spin states is concerned in conjunction with the most relevant conformation of the CpdI-His42-Arg38 complex, we would like to point out that our investigation basically predicts, in agreement with previous data at the CPCM/ QCISD level, the doublet state at lower energies. In turn, it is worth noting that the doublet-quartet energy difference distribution is found to be in line with recent QM/MM data. Finally, our results show that the overall system, HRP plus solvent, is capable of discriminating (on the basis of its own conformational flexibility driven by entropic effects) between different magnetic eigenstates of the active site, slightly favoring the more polar one. In light of this finding, it would be extremely appealing to investigate, in the near future, whether and to what extent a change in the conformational space accessible by thermal fluctuations (for example anchoring the HRP enzyme on nanostructured surfaces) might modulate the electronic properties and reactivity of the CpdI redox unit. Acknowledgment. We gratefully acknowledge Prof. Sason Shaik for very helpful discussions about HRP enzyme. All the computations were performed at the Supercomputing Centre for University and Research, CASPUR (Rome). References and Notes (1) (a) Dunford, H. B. Heme Peroxidases; Wiley-VCH: New York, 1999. (b) Gajhede, M.; Schuller, D. J.; Henricksen, A.; Smith, A. T. Nat.

6824

J. Phys. Chem. B, Vol. 114, No. 20, 2010

Struct. Biol. 1997, 4, 1032–1038. (c) Watanabe, Y. High-Valent Intermediates. In The Porphyrin Handbook; Kadish, K. M., Smith, K. M., Guilard, R., Eds.; Academic Press: New York, 2000; Vol. 4, Chapter 30, pp 97118. (d) Marnett, L. J.; Kennedy, T. A. In Cytochrome P450: Structure, Mechanisms and Biochemistry, 2nd ed.; Ortiz de Montellano, P. R., Ed.; Plenum Press: New York, 1995; pp 49-80. (e) Watanabe, L.; de Moura, P. R.; Bleicher, L.; Nascimento, A. S.; Zamorano, L. S.; Calvete, J. J.; Sanz, L.; Perez, A.; Bursakov, S.; Roig, M. G.; Shnyrov, V. L.; Polikarpov, I. J. Struct. Biol. 2010, 169, 226–242, and quoted references cited therein. (2) Poulos, T. L.; Kraut, J. J. Biol. Chem. 1980, 255, 8199–8205. (3) Douzou, P. In Oxidases and Related Redox Systems; King, T. E., Mason, H. S., Morrison, M., Eds.; University Park Press, Baltimore, 1973, vol I. (4) (a) Dolphin, D.; Forman, A.; Borg, D. C.; Fajer, J.; Felton, R. H. Proc. Natl. Acad. Sci. U.S.A. 1971, 68, 614–618. (b) Dolphin, D.; Felton, R. H. Acc. Chem. Res. 1974, 7, 26–32. (c) DiNello, R. K.; Dolphin, D. H. J. Biol. Chem. 1981, 256, 6903–6912. (d) Moss, T. H.; Ehrenberg, A.; Bearden, A. J. Biochemistry 1969, 8, 4159–4162. (e) Rodriguez-Lopez, J. N.; Lowe, D. J.; Hernandez-Ruiz, J.; Hiner, A. N. P.; Garcia-Canovas, F.; Thorneley, R. N. F. J. Am. Chem. Soc. 2001, 123, 11838–11847. (f) Schulz, C. E.; Rutter, R.; Sage, J. T.; Debrunner, P. G.; Hager, L. P. Biochemistry 1984, 23, 4743–4754. (g) de Ropp, J. S.; Sham, S.; Asokan, A.; Newmyer, S.; Ortiz de Montellano, P. R.; La Mar, G. N. J. Am. Chem. Soc. 2002, 124, 11029–11037. (h) Kincaid, J. R.; Zheng, Y.; Al-Mustafa, J.; Czarnecki, K. J. Biol. Chem. 1996, 271, 28805–28811. (i) Berglund, G. I.; Carlsson, G. H.; Smith, A. T.; Szoeke, H.; Henriksen, A.; Hajdu, J. Nature 2002, 417, 463–468. (5) (a) Zazza, C.; Aschi, M.; Palma, A. Chem. Phys. Lett. 2006, 428, 152–156. (b) Derat, E.; Shaik, S. J. Phys. Chem. B 2006, 110, 10526– 10533. (c) Derat, E.; Shaik, S.; Rovira, C.; Vidossich, P.; Alfonso-Prieto, M. J. Am. Chem. Soc. 2007, 129, 6346–6347. (d) Derat, E.; Cohen, S.; Shaik, S.; Altun, A.; Thiel, W. J. Am. Chem. Soc. 2005, 127, 13611–13621. (e) Roth, J. P.; Cramer, C. J. J. Am. Chem. Soc. 2008, 130, 7802–7803. (f) Zazza, C.; Sanna, N.; Tatoli, S.; Aschi, M.; Palma, A. Int. J. Quantum Chem. 2010, 110, 352–357. (g) Harris, D. L. Curr. Opin. Chem. Biol. 2001, 5, 724–735. (h) Harris, D. L.; Loew, G. H. J. Porphyrins Phthalocyanins 2001, 5, 334–344. (i) Zazza, C.; Amadei, A.; Palma, A.; Sanna, N.; Tatoli, S.; Aschi, M. J. Phys. Chem. B 2008, 112, 3184–3192. (6) Pople, J. A.; Head-Gordon, M.; Raghavachari, K. J. Chem. Phys. 1985, 82, 299. (7) Cossi, M.; Rega, N.; Scalmani, G.; Barone, V. J. Comput. Chem. 2003, 24, 669–81. (8) (a) Lonsdale, R.; Harvey, J. N.; Mulholland, A. J. J. Phys. Chem. B 2010, 114, 1156–1162. (b) Shaik, S.; Kumar, D.; de Visser, S. P.; Altun, A.; Thiel, W. Chem. ReV. 2005, 105, 2279–2328. (c) Lu, H. P. Acc. Chem. Res. 2005, 38, 557–565. (d) Aschi, M.; Zazza, C.; Spezia, R.; Bossa, C.; Di Nola, A.; Paci, M.; Amadei, A. J. Comput. Chem. 2004, 7, 974–984. (e) Lodola, A.; Mor, M.; Zurek, J.; Tarzia, G.; Piomelli, D.; Harvey, J. N.; Mulholland, A. J. Biophys. J. 2007, 92, 20–22. (9) Mokdad, A.; Miksovska, J.; Larsen, R. W. Biochim. Biophys. Acta 2009, 1794, 1558–1565.

Zazza et al. (10) Tatoli, S.; Zazza, C.; Sanna, N.; Palma, A.; Aschi, M. Biophys. Chem. 2009, 141, 87–93. (11) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F. In Intermolecular Forces; Pullman, B., Ed.; D. Reidel Publishing Company: Dordrecht, 1981; pp 331-342. (12) Laberge, M.; Huang, Q.; Schweitzer-Stenner, R.; Fidy, J. Biophys. J. 2003, 84, 2542–2552. (13) van Gunsteren, W. F.; Billeter, S. R.; Eising, A. A.; Hunemberger, P. H.; Kruger, P.; Mark, A. E.; Scott, V. R. P.; Tironi, I. G. Biomolecular Simulation: The GROMOS96 Manual and User Guide; vdf Hochschlverlag AG an der ETH Zurich: Zurich, 1996. (14) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361– 367. (15) Parr, R. G.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989. (16) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (17) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 299–310. (18) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650–654. (19) Amadei, A.; Chillemi, G.; Ceruso, M.; Grottesi, A.; Di Nola, A. J. Chem. Phys. 2000, 112, 9–23. (20) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Di Nola, A. J. Chem. Phys. 1994, 81, 3684–3690. (21) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Frajie, J. G. E. M. J. Comput. Chem. 1997, 18, 1463–1472. (22) Darden, T. A.; York, D. M.; Pedersen, L. G. J. Chem. Phys. 1993, 98, 10089–10092. (23) Lindhal, E.; Hess, B.; Van der Spoel, D. J. Mol. Model. 2001, 7, 306–317. (24) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, reVision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (25) Amadei, A.; Linssen, A. B. M.; Berendsen, H. J. C. Proteins: Struct., Funct., Genet. 1993, 17, 412–425. (26) Zazza, C.; Palma, A.; Amadei, A.; Sanna, N.; Tatoli, S.; Aschi, M. Faraday Discuss 2010, 145, 107-120.

JP101033W