Computational Study of Catalytic Reaction of Quercetin 2,4

version 3.5; Science & Technology Facilities Council: Swindon, UK; www.chemshell.org ..... For a more comprehensive list of citations to this arti...
0 downloads 0 Views 4MB Size
Article pubs.acs.org/JPCB

Computational Study of Catalytic Reaction of Quercetin 2,4Dioxygenase Toru Saito,*,†,‡ Takashi Kawakami,† Shusuke Yamanaka,†,‡ and Mitsutaka Okumura†,‡ †

Department of Chemistry, Graduate School of Science, Osaka University, 1-1 Machikaneyama, Toyonaka, Osaka 560-0043, Japan Core Research for Evolutional Science and Technology (CREST), Japan Science and Technology (JST) Agency, Kawaguchi, Saitama 332-0012, Japan



S Supporting Information *

ABSTRACT: We present a quantum mechanics/molecular mechanics (QM/ MM) and QM-only study on the oxidative ring-cleaving reaction of quercetin catalyzed by quercetin 2,4-dioxygenase (2,4-QD). 2,4-QD has a mononuclear type 2 copper center and incorporates two oxygen atoms at C2 and C4 positions of the substrate. It has not been clear whether dioxygen reacts with a copper ion or a substrate radical as the first step. We have found that dioxygen is more likely to bind to a Cu2+ ion, involving the dissociation of the substrate from the copper ion. Then a Cu2+-alkylperoxo complex can be generated. Comparison of geometry and stability between QM-only and QM/MM results strongly indicates that steric effects of the protein environment contribute to maintain the orientation of the substrate dissociated from the copper center. The present QM/MM results also highlight that a prior rearrangement of the Cu2+alkylperoxo complex and a subsequent hydrogen bond switching assisted by the movement of Glu73 can facilitate formation of an endoperoxide intermediate selectively. (Glu73).6−9 According to the enzyme−substrate (ES) complex solved under anaerobic conditions (PDB code: 1H1I),8 the proton of the C3-OH group of Que is considered to be transferred to a carboxyl oxygen atom of Gly73 when bound to the copper center (see 1 in Scheme 2). Hernandez-Ortega et al. showed that the initial proton abstraction from the substrate can be seen in some cofactor-free dioxygenases that catalyze a direct dioxygenation of a substrate without involvement of a metal ion.18 A hydrogen bond (H-bond) between the deprotonated substrate (Que−) and protonated Glu73 plays a role in keeping Que− in a monodentate coordination. It has been suggested that the enzymatic reaction involves an endoperoxide intermediate as a precursor of the product.7−17 Some synthetic model complexes with bidentate coordination of quercetin indeed may lead to a different product through a 1,2-dioxetane, showing chemiluminescence.15−17 Density functional theory (DFT) calculations revealed that both endoperoxide and 1,2-dioxetane intermediates can be obtained in a nonenzymatic reaction.19 Based on the X-ray structure, two possible reaction mechanisms for the first step have been proposed as shown in Scheme 2.8,9 Since the initial Cu2+ complex in the active site is incapable of binding to O2 and activating it, both mechanisms involve formation of a tautomer with a Cu+ ion and substrate

1. INTRODUCTION Quercetin (3,3′,4′,5,7-pentahydroxyflavone) (Que) belongs to flavonoids and shows antioxidant and antibacterial activities.1−3 Consumption of quercetin contained in vegetables and fruits is considered to be effective for prevention human diseases such as cancer and heart disease. Quercetin 2,4-dioxygenase (2,4QD) is an enzyme that catalyzes the cleavage of the phenyl ring of quercetin to form the phenolic carboxylic acid ester and carbon monoxide in quercetin metabolic process as shown in Scheme 1.4−18 Scheme 1. Dioxygenation Reaction of Quercetin Catalyzed by 2,4-QD

Here, we call the enzyme not 2,3-dioxygenase but 2,4dioxygenase on the basis of its function,6 which incorporates two oxygen atoms at C2 and C4 positions. 2,4-QD from Aspergillus japonicus is a homodimer in the crystal. The electron paramagnetic resonance (EPR) spectroscopy confirms that each monomer has a type 2 copper protein with a Cu2+ complex, the active site of which consists of three histidine residues (His66, His68, and His112) and one glutamic acid © 2015 American Chemical Society

Received: April 13, 2015 Revised: May 17, 2015 Published: May 19, 2015 6952

DOI: 10.1021/acs.jpcb.5b03564 J. Phys. Chem. B 2015, 119, 6952−6962

Article

The Journal of Physical Chemistry B Scheme 2. Proposed Reaction Mechanisms for Ring-Cleaving Reaction of 2,4-QDa

a

In the first step, O2 can bind to a Cu+ Que· tautomer either at the copper ion (i) or at the C2 atom of Que (ii).

In the present study, we have investigated the dioxygenation reaction using a combined quantum mechanics/molecular mechanics (QM/MM) approach.32−37 The QM calculations with a cluster model, namely, QM-only calculations, were also employed for comparison. We would like to propose a plausible reaction pathway by addressing unanswered questions in particular; whether O2 first binds to Cu or Que and why 2,4QD gives an endoperoxide selectively. We have also considered the effect of surrounding protein residues on geometric and electronic structures. As far as we know, this is the first QM/ MM study of the dioxygenation reaction catalyzed by 2,4-QD.

radical (Que·) pair (1taut), associated with one electron transfer from Que− to the Cu2+ ion. According to some similarity with iron(III)-containing intradiol catechol dioxygenases,20−29 O2 can bind to 1taut either at the copper ion (i) or the C2 atom of Que (ii). In the Cu attack (oxygen activation) mechanism (i), dioxygen first binds to Cu+ at the position trans to His68, followed by formation of a [Cu2+-superoxo (Que·)] complex (2A). The subsequent attack on the C2 atom gives rise to a Cu2+-alkylperoxo species (3A) which may undergo nucleophilic attack of a copper-bound peroxo species on the C4 atom, generating an endoperoxide intermediate (EP). The Que attack (substrate activation) mechanism (ii) starts with a direct attack of O2 at the C2 atom of Que·, which is slightly pyramidalized and has weak sp3 character, providing a [Cu2+ Que-peroxo] complex (2B). Formation of EP from 2B can proceed with nucleophilic attack of a substrate-bound peroxo species on C4 atom. In EP, two C−C bonds would be cleaved upon the O−O bond scissions, leading to a product complex (P). There is no consensus on which of these two mechanisms is plausible. Steiner et al. support the substrate activation mechanism (ii) on the basis of structural features.8 They suggested that formation of EP from 3A may not be feasible. Given the position of Que is constrained by the H-bond with Glu73 and coordination bond with Cu2+, even the proximal oxygen atom in 3A is ca. 3.5−4.0 Å apart from the C4 atom. Siegbahn examined the stability for 2A and 2B by using DFT calculations with an active site model.30 2A was calculated to be more stable than 2B by 14 kcal/mol, but transition states toward 2A and 2B and corresponding activation barriers were not computed. Along with the mechanism (i), he explored a reaction mechanism, and concluded that formation of 3A from 2A would be the rate-limiting step. The preceding nucleophilic attack on C4 atom can be a single-step reaction during which the H-bond with Glu73 switches from O3 to O4 spontaneously, upon formation of the O−C4 bond. Antonczak et al.19,31 tested the mechanism (ii) using a model consisting of only O2 and Que (without the copper complex). However, DFT calculations with such a small model may prefer a 1,2dioxetane to an endoperoxide as mentioned above.

2. COMPUTATIONAL DETAILS 2.1. System Setup. The initial X-ray structure of an ES complex was taken from one monomer, chain A in the Protein Data Base (PDB code: 1H1I, resolution 1.75 Å).8 N-acetyl-Dglucosamines and (4S)-2-methyl-2,4-pentanediols were not included in the model. The missing residues 1−2 and 155− 158 were generated by MODELLER (ver. 9.11).38 The protonation states of Asp, Glu, and His residues were assigned by means of estimated pKa values using the PROPKA program39 together with visual inspection. The protonation states were defined as follows: Glu91 is protonated. His71 as well as the metal-ligating His66, His68, and His112 are protonated in the δ position, while His21 and His201 are protonated in the ε position. The rest His13, His24, His57 and His148 are doubly protonated. The CHELPG potential fit40 calculated at the B3LYP/def2-SVP level was used to determine the atomic charges on Que. The Lennard-Jones parameters proposed by Ungar et al.41 were used for the Cu2+ ion. The positions of hydrogen atoms were determined by the HBUILD module in the CHARMM program (ver. c38b2).42,43 The classical molecular dynamics (MD) simulations were performed with the CHARMM program. The model was solvated in a water sphere of 35 Å radius centered at the center of mass with stochastic boundary potentials.44 The solvent water molecules were relaxed by a 50 ps MD simulation, after which the system was resolvated. We repeated this procedure five times. The initial system had a total charge of −12e. The 6953

DOI: 10.1021/acs.jpcb.5b03564 J. Phys. Chem. B 2015, 119, 6952−6962

Article

The Journal of Physical Chemistry B

frequencies were calculated analytically to check whether all stationary points are a local minimum with no imaginary frequencies or a transition state with one imaginary frequency. It should be noted that small imaginary frequencies were obtained in some cases due to freezing atoms. All QM-only computations were performed using the Gaussian 09 program.

net charge of the system was then neutralized, replacing randomly selected solvent water molecules (at least 5.5 Å away from any protein atoms) by Ca2+ and Cl− ions. As depicted in Figure 1, the system consisted of 18,907 atoms, including 4148

3. RESULTS AND DISCUSSION 3.1. ES Complex with Dioxygen. The geometry for the ES complex first was optimized to check whether the O3 oxygen atom or the carboxylate group of Glu73 is protonated. In both QM-only and QM/MM calculations, the proton of the substrate hydroxyl group was transferred to Glu73 spontaneously during the geometry optimization. It suggests that Glu73 energetically favors a protonated state as in the case of the previous QM study.30 The QM/MM-optimized O···O hydrogen bond length of 2.51 Å was in good agreement with the distance in the X-ray structure (2.43 Å), whereas the distance obtained with the QM-only method (2.62 Å) overestimates it. The Cu−N(His68) distances of 2.23 Å (QM/MM) and 2.17 Å (QM) are elongated as compared to the initial value of 2.11 Å, while the Cu−O(Glu73) distances of 2.11 Å (QM/MM) and 2.21 Å (QM) is rather shorter than the distance in the X-ray structure (2.30 Å). The spin populations on the Cu atom are 0.61 and 0.56 calculated at the QM/MM and QM-only levels, supporting a Cu2+ ion. The rest of spin density are distributed mainly on the O3 atom of Que. The sum of the QM/MM and QM-only energy of the ES complex and that of the isolated triplet dioxygen is defined as the reference energy (1 and 1QM). We will discuss potential energies relative to 1 and 1QM. In the QM/MM calculation, the O2 molecule was added to the model around the ES complex (ES-O2 complex), and then it was reoptimized for the doublet and quartet states. There is little interaction between O2 and ES complex judging from the calculated spin populations (see Figure 2). The optimized

Figure 1. QM/MM model with bulk water molecules (left) and QM region (right). The QM region consists of 69 atoms (including link atoms) and has a total charge of +1e. Atoms marked with the asterisks were kept frozen in QM-only calculations.

TIP3P water molecules. Then, a 0.5 ns MD simulation run at 300 K with a time step of 1 fs. During the MD simulations, the QM atoms (see below) were kept fixed. 2.2. QM-Only and QM/MM Calculations. QM/MM geometry optimizations of 2,4-QD start from the final snapshot from the MD trajectory. The QM subsystem consists of Cu, Que, the acetic acid as a model for Glu73, and the imidazole groups of ligating histidine residues (His66, His68, His112). The spin unrestricted B3LYP45−47 (UB3LYP) method was used as with the previous QM studies.19,30,31 Geometry optimizations were performed for both the doublet and quartet states with the def2-SVP basis set.48,49 Single point energy calculations were additionally carried out with the def2-TZVP basis set48,49 based on the QM/MM-optimized geometries. We applied Yamaguchi’s method50,51 to correct the energies of the doublet state. The QM/MM computations were carried out using the ChemShell program.52,53 We used the electronic embedding scheme, in which the fixed MM atomic partial charges are included in the QM Hamiltonian.54 No cutoff for nonbonded interactions was introduced. The link atom method was used at the QM/MM boundary, and the charge shift model was applied to avoid overpolarization.52 We used the Gaussian 0955 and the DL_POLY56 programs to describe the QM and MM regions, respectively. The CHARMM22 force field57 was applied to the MM computations. The forces on the MM point charges were computed from the electric field at the MM atoms58,59 Geometry optimizations were performed using DLFIND,60 which is implemented in the ChemShell program. The positions of all MM atoms more than 12 Å away from the copper ion were kept fixed during the optimizations. The quasiNewton limited-memory Broyden−Fletcher−Goldfarb−Shannon (L-BFGS) method61 was used to locate minima. Transition states were optimized using the climbing image nudged elastic band (CI-NEB) method.62,63 The numerical frequency calculations for core atoms were performed to confirm whether obtained TS structures have one negative eigenvalue. QM-only calculations were also carried out for comparison at the same level of theory (UB3LYP/def2-TZVP//UB3LYP/ def2-SVP). The number of the atoms was the same as the QM region used for the QM/MM computations. Several atoms of ligating residues marked with the asterisks were kept frozen64,65 during the geometry optimizations (see Figure 1). Vibrational

Figure 2. QM/MM-optimized structure of the ES−O2 complex (1′). The selected optimized distances are presented in Å. The spin populations on Cu, O2, and Que for doublet and quartet (in parentheses) states are also presented.

geometries for both spin states are calculated to be identical to each other, but the total QM/MM energy for 21′ is 0.2 kcal/ mol lower than that for 41′. The QM/MM-optimized ES-O2 complex 21′ (41′) lies below 1 by 1.6 (1.4) kcal/mol. As can be seen in Figure 2, the distances between Cu and proximal oxygen (Op) and distal oxygen (Od) atom are 3.85 and 4.11 Å, respectively. The C2 atom of Que is located 3.36 Å away from 6954

DOI: 10.1021/acs.jpcb.5b03564 J. Phys. Chem. B 2015, 119, 6952−6962

Article

The Journal of Physical Chemistry B

Figure 3. Energy profiles for formation of a Cu-alkylperoxo intermediate (4a) via the oxygen activation (Cu attack) mechanism (solid line) and for formation of a 1,2-dioxetane intermediate (2b) via the substrate activation (Que attack) mechanism (dashed line) obtained by the QM/MM calculations. All stationary points were located at the (UB3LYP/def2-SVP)/CHARMM22 level, and relative energies (in kcal/mol) for the doublet and quartet (in parentheses) spin states computed at the (UB3LYP/def2-TZVP)/CHARMM22 level are also presented.

occurs, resulting from an elongated Cu···O3 distance. Formation of 2,42a is endothermic by 7.1 and 7.9 kcal/mol relative to 1 for the doublet and quartet states, respectively. The Cu−O3 bond is completely broken (ca. 2.6 Å) after the Cu−Op bond formation. From the ES-O2 complex, no transition state was located that directly generates another end-on (η1) complex [Cu2+−Od−Op−· Que·] (3a). Nevertheless, the conversion from 22a to 23a via 2TS2 is found to be feasible with a small activation barrier of 0.4 kcal/ mol. The resultant 23a is more stable than 22a by 3.8 kcal/mol. During the step, the Od···C2 distance decreases from 3.87 to 2.67 Å. The difference in spin density between 22a and 23a suggests partial electron transfer from Que to dioxygen. No intermediate corresponding to 3a was located for the quartet state in spite of many efforts. The O2 molecule started dissociating from the Cu ion, once in a position where the Od atom faced the C2 atom. In the QM-only model, meanwhile, the direct formation of 3a (denoted as 2,43aQM) was observed through 2TS2QM

the Op atom. The O2 molecule may bind to any of these sites and will give rise to a Cu−Op, Cu−Od, or Op−C2 bond. The former two corresponds to the oxygen activation (Cu attack) mechanism (i) and the latter one is the substrate activation (Que attack) mechanism (ii). 3.2. O2 Activation (Cu Attack) Mechanism. The mechanism (i) starts from the Cu−Op adduct. Formation of an end-on (η1) [Cu2+−Op−Od−· Que·] complex goes differently between the QM-only and QM/MM models. In the QM/ MM model, an end-on (η1) [Cu2+−Op−Od−· Que·] complex (2,42a) is generated via 2,4TS1 for both spin states (see Figure 3). The quartet state is the ground state in TS1, but potential energy surfaces for two spin states are close in energy (within 0.2 kcal/mol). These transition states are 9.1 (2TS1) and 8.9 (4TS1) kcal/mol higher than 1. As listed in Table 1, in 2TS1 (4TS1), the Cu···O3 and Cu···Op distances are 2.19 (2.38) and 2.23 (2.88) Å, and the spin population on Cu, O2, and Que are 0.20 (0.43), 1.60 (1.60), and −0.85 (0.83), respectively. This reflects that partial electron transfer from Que− to the Cu2+ ion 6955

DOI: 10.1021/acs.jpcb.5b03564 J. Phys. Chem. B 2015, 119, 6952−6962

Article

The Journal of Physical Chemistry B

Table 1. Optimized Structural Parameters and Spin Populations for Stationary Points Involved in Formation of a CuAlkylperoxo Intermediate via the Oxygen Activation (Cu Attack) Mechanism (i)a parametersb

spin population

species

Cu−Op

Op−Od

Od−C2

Cu−O3

Cu···O4

Cu

O2

Que

TS1 4 TS1 2 2a 4 2a 2 TS2 2 3a 2 TS3 2 4a 2 TS2QM 4 TS2QM 2 3aQM 4 3aQM 2 TS3 QM 2 4aQM

2.23 2.18 2.05 2.05 2.04 2.02 1.98 1.97 2.16 2.18 2.02 2.03 1.97 1.93

1.23 1.24 1.26 1.26 1.26 1.27 1.29 1.41 1.23 1.24 1.25 1.25 1.30 1.41

3.89 3.93 3.87 3.87 3.67 2.67 2.16 1.48 3.04 3.12 3.13 3.62 2.08 1.46

2.19 2.38 2.57 2.60 2.62 2.37 2.49 2.86 2.16 2.26 2.55 2.85 2.54 2.85

3.56 3.56 3.61 3.65 3.63 3.46 3.50 3.54 3.71 3.69 3.91 4.21 3.77 3.98

0.20 0.43 0.40 0.43 0.42 0.46 0.55 0.61 −0.18 0.48 0.37 0.39 0.52 0.59

1.60 1.60 1.47 1.48 1.47 1.28 0.87 0.22 1.76 1.65 1.49 1.53 0.81 0.22

−0.85 0.83 −0.96 0.98 −0.98 −0.84 −0.55 0.01 −0.55 0.72 −0.94 0.98 −0.47 0.00

2

a

These stationary points were located using the (UB3LYP/def2-SVP)/CHARMM22 and UB3LYP/def2-SVP computations for the QM/MM and QM-only models, respectively. bIn Å.

complex (2B) but a [Cu+ Que-superoxo] complex. As in the case of the mechanism (i), the QM-only and QM/MM models differ in a reaction mechanism. The Op−C2 bond formation in the QM/MM model occurs via one electron transfer from Que− to Cu2+ together with the complete Cu−O3 bond cleavage. The Cu···O3 distance is 2.00 Å in the ES-O2 complex, which elongates to 3.45 Å in the transition state (2TS2b) as listed in Table 2. In the doublet state, this gives rise to a closedshell Cu+ ion and a substrate radical Que·(↓). The O2 (↑↑) molecule attacks the C2 atom of Que·(↓) to form a Que−O−O (↑) species 22b (see Figure S1 in the Supporting Information), while the Cu+ ion is not involved. The spin populations on Cu (0.00), O2 (0.99), and Que (0.01) confirm its electronic structure. The Op−C2 adduct cannot happen in the quartet state because a triplet coupling Od(↑)···C2(↑) is energetically unfavorable. Formation of 22b requires a substantially high barrier of 21.8 kcal/mol, and 22b is still too unstable as compared to 22a and 23a. No intermediate structure was found, where the Cu−O3 bond remains intact like 2B. All attempts to locate a direct pathway from 22b to an endoperoxide intermediate (2EP) ended up forming only 24a. It seems like that formation of a Cu−Od bond takes place with a small activation barrier of 1.0 kcal/mol, prior to the Od−C4 bond formation. A 1,2-dioxetane intermediate through the C2−Op− Od−C3 bond formation was not located, due to the carbonyl C3O3 group lying out-of-plane of the C-ring in 22b (see Figure S2 in the Supporting Information). The QM-only model provided a substantially different intermediate structure (22bQM), in which the Cu−O3 bond still remains (2.11 Å) as 2B. It is obvious, as shown in Figure 5b, that the B-ring has resulted in a significant out-of-plane distortion with a RMSD value of 0.68 Å. On the other hand, the QM/MM-optimized structure (22b) with a RMSD of 0.36 Å does not differ from the X-ray structure significantly, except for the dissociation of the Cu−O3 bond. Starting from 22bQM, we located a pathway neither to 2EP nor a Cu2+-alkylperoxo complex (24a). Instead, a substrate-bound superoxo species is likely to attacks on the other polarized carbonyl group (C3 O3) in plane with the C2 atom, leading to a 1,2-dioxetane intermediate (2DioQM). The reaction may yield another

(4TS2QM) with an activation barrier of 10.4 (10.3) kcal/mol as illustrated in Figure 4. In 2TS2QM (4TS2QM), the Cu···O3 and Cu···Op distances are 2.16 (2.26) and 2.16 (2.18) Å, and the spin population on Cu, O2, and Que are −0.18 (0.48), 1.76 (1.65), and −0.55 (0.72), respectively (see Table 1). Comparison between 23a and 23aQM reveals that the interval between the Cu ion and Que are significantly different, as indicated by the Od···C2 and Cu···O4 separations. To examine the difference between two approaches, the root-mean-squared deviation (RMSD)66 of the position of the substrate from the crystal structure was calculated. The RMSD values are 0.16 and 0.31 Å for the 23a and 23aQM, respectively. Figure 5a shows the overlay of the computationally optimized and the crystal structures together with Met51 and Phe75. It can be seen that the QM-only model reduces the steric repulsion by separating the substrate from O2 artificially, whereas in the QM/MM model, the nonbonding interactions between Que and the surrounding residues (eg, π−π interaction with Phe75 and CHπ interaction with Met51) keep contributing to impose the orientation of the B-ring of Que. After 23a is formed, the Od atom of the copper-bound superoxo species attacks the C2 atom of Que to yield a Cu2+alkylperoxo complex (24a) at the QM/MM level. The corresponding transition state (2TS3) takes only 1.5 kcal/mol from 23a, which leads to a singlet coupling67−69 Od(↑)···C2(↓) (see Figure S1 in the Supporting Information). The Op-Od bond length of 1.41 Å and decreased spin populations on both O2 (0.22) and Que (0.00) moiety are evidence for the electronic structure of 24a. In the quartet state, however, geometry optimizations fail to locate a minimum because a triplet coupling Od(↑)···C2(↑) is energetically unfavorable. Therefore, results of geometry optimizations for the doublet state will be discussed hereafter. The same tendency is found in the QM-only calculations, except that formation of 24aQM requires a higher activation barrier than TS2QM (12.2 vs 10.4 kcal/mol). The RMSD values are 0.22 and 0.33 Å for 24a and 2 4aQM, respectively. 3.3. Substrate Activation (Que Attack) Mechanism. The substrate activation (Que attack) mechanism (ii) initiates the direct Op-C2 adduct, leading to not a [Cu2+ Que-peroxo] 6956

DOI: 10.1021/acs.jpcb.5b03564 J. Phys. Chem. B 2015, 119, 6952−6962

Article

The Journal of Physical Chemistry B

Figure 4. Energy profiles for formation of a Cu-alkylperoxo intermediate (4aQM) via the oxygen activation (Cu attack) mechanism (solid line) and for formation of a 1,2-dioxetane intermediate (DioQM) via the substrate activation (Que attack) mechanism (dashed line) obtained by the QM-only (UB3LYP/def2-SVP) calculations. Relative energies (in kcal/mol) calculated at the UB3LYP/def2-TZVP for the doublet and quartet (in parentheses) spin states are also presented.

iron-dependent dioxygenases.20,21,29 During the conversion from 24a to 25a, the Od···C4 distance decreased from 3.67 to 2.59 Å as summarized in Table 3. The O(Glu73)···O4 separation decreases from 4.89 to 3.71 Å. This rearrangement facilitates the following H-bond switching. The Cu···O4 separation also reduces to 2.89 Å. The subsequent H-bond switching assisted by the movement of the flexible Glu73 (2TS5) with a small activation barrier of 2.1 kcal/mol. The O(Glu73)···O4 distance has changed from 4.09 to 2.76 Å. The energetics calculated with the QM-only method seems similar to those with the QM/MM method. However, the geometric feature of the corresponding 25aQM is significantly different. In contrast to the QM/MM results, the Cu···C4 distance does not shorten in the conversion process. The O(Glu73)···O4 distance of 5.06 Å is still long in 25aQM. Comparison between 25aQM and 2 6aQM suggests that the structural difference originates from the movement of Que rather than that of Glu73 (see Figure S3 in the Supporting Information), although we could not find a structure that qualifies for a transition state. The switched H-bond in 26a can make the C4O4 carbonyl group more polarized, and electrophilicity of the C4 atom is

product because 2DioQM is stable by 1.9 kcal/mol than 4aQM. Comparison between QM/MM and QM-only results show that noncovalent interactions are important to not only stabilize the substrate dissociated from the copper ion but also impede a nonenzymatic reaction. The present QM/MM calculations strongly indicate that the mechanism (ii) would have to be ruled out as the first step of the reaction. 3.4. Reaction Mechanism for Endoperoxide (EP) Formation. Unlike the previous QM study, no single-step transformation from a Cu2+-alkylperoxo complex (24a or 2 4aQM) to an endoperoxide intermediate was observed for both QM-only and QM/MM models. Discussions about the discrepancy are beyond the scope of this study because a different basis set and a cluster model were used for geometry optimizations.70 As summarized in Figure 6, we would like to propose a multistep mechanism. From the QM/MM calculations, a large geometric rearrangement stemmed from an O−O bond rotation (2TS4) that requires a high activation barrier of 15.1 kcal/mol relative to 24a, providing a rearranged Cu2+-alkylperoxo intermediate (25a). It was also found that a Criegee-type rearrangement is unlikely to occur in contrast to 6957

DOI: 10.1021/acs.jpcb.5b03564 J. Phys. Chem. B 2015, 119, 6952−6962

Article

The Journal of Physical Chemistry B

Figure 5. Overlay of the substrate in (a) an end-on (η1) [Cu2+−Op−Od−· Que·] complex and in (b) a [Cu+ Que-superoxo] complex with ligands and nearby residues (Met51 and Phe75): the X-ray crystal structure, QM-only optimized structure (yellow), and QM/MM-optimized structure (red).

Table 2. Optimized Structural Parameters and Spin Populations for Stationary Points Involved in the Substrate Activation (Que Attack) Mechanism (ii)a parametersb

spin population

species

Cu−Od

Op−Od

Op−C2

Cu···O3

Cu···O4

C3−C2−Op−Od

Cu

O2

Que

TS2b 2 TS2bQM 2 2b 2 2bQM 2 TS2b4a 2 TSdioQM 2 DioQM

4.03 5.01 4.21 4.42 3.79 3.96 3.27

1.25 1.26 1.30 1.31 1.30 1.40 1.47

1.94 1.77 1.53 1.51 1.55 1.41 1.44

3.45 2.09 3.63 2.11 3.52 2.11 1.94

3.44 3.94 3.41 4.10 3.54 4.22 4.07

129.1 −104.5 151.2 −65.0 124.8 −31.8 4.1

0.00 −0.22 0.00 −0.16 0.00 0.43 0.69

1.20 1.13 0.99 0.99 1.20 0.50 0.00

−0.21 0.13 0.01 0.19 0.21 −0.06 0.08

2

a

These stationary points were located using the (UB3LYP/def2-SVP)/CHARMM22 and UB3LYP/def2-SVP computations for the QM/MM and QM-only models, respectively. bIn Å for lengths and deg for angles.

4. CONCLUSION In the present study, we have studied the dioxygenation reaction of quercetin catalyzed by 2,4-QD. Both the O2 activation (Cu attack) and substrate activation (Que attack) mechanisms have been demonstrated using QM/MM and QMonly methods. We would like to emphasize that the substrate must dissociate from the metal center to react with O2, unlike the suggestion based on the geometric feature. The present QM/MM results have shown that steric effects of the protein environment are able to maintain the orientation of the substrate dissociated from the copper ion. The first step of the dioxgenation reaction favors the oxygen activation mechanism (i) with a barrier of ca. 9 kcal/mol relative to the reactant. Then, the Cu2+-superoxo binds to Que, leading to a Cu2+alkylperoxo complex. Note that this alkylperoxo complex is going to be generated even if an unstable intermediate happens to be formed through the substrate activation mechanism (ii). The fact that a metal−alkylperoxo species may be generated as a intermediate in 2,4-QD is consistent with the recent report that revealed the presence of an alkylperoxo intermediate in an intradiol dioxygenase.29 For the following step, we propose a stepwise reaction mechanism in which the Od−C4 bond formation can occur reasonably. It involves a rearrangement of the Cu−Op−Od−C2 moiety, which would be the rate-

promoted. The preceding nucleophilic attack is allowed to occur, leading to 2EP via 2TS6. 2TS6 is described by the Cu− Op bond scission (2.24 Å) and Op−C4 bond formation (1.88 Å). The barrier height is 12.2 kcal/mol. It also involves a Cu− O4 bond formation, and the Cu−O4 bond length is 1.98 Å in 2 EP. The simultaneous cleavage of O−O and two C−C bonds through 2TS7 gives rise to a phenolic carboxylic acid ester and carbon monoxide as a product complex (2P). The QM-only calculations provide a similar reaction pathway concerning 2 EPQM and 2PQM formation from 26aQM. Overall, the ratelimiting step is found to be the rearrangement of peroxo group with a barrier of 16.8 kcal/mol relative to 1. This is in good agreement with the activation energy of 14.5−17.3 kcal/mol observed for a synthetic model complex [Cu2+(fla)(idpa)]ClO4 (fla = falvonolate; idpa =3,3′-iminobis(N,N-dimethylpropylamine)),15 although the reaction does not necessarily proceed with the same way. Since there have not been experimentally observed activation barrier for 2,4-QD so far, direct comparison between theory and experiment will be a future challenge. We should also address the effects of dispersion interactions for QM-QM atom pairs71,72 and performance of different DFT exchange-correlation functionals59,73 on the energies and structures of stationary points. 6958

DOI: 10.1021/acs.jpcb.5b03564 J. Phys. Chem. B 2015, 119, 6952−6962

Article

The Journal of Physical Chemistry B

Figure 6. Energy profile for formation of a product complex (2P or 2PQM) from a Cu2+-alkylperoxo intermediate (4a or 4aQM) via an endoperoxide intermediate (2EP or 2EPQM). All stationary points were located using the (UB3LYP/def2-SVP)/CHARMM22 and UB3LYP/def2-SVP computations for the QM/MM and QM-only models, respectively. Relative energies (in kcal/mol) for the doublet spin states computed at the (UB3LYP/def2-TZVP)/CHARMM22 and UB3LYP/def2-TZVP (in italics, blue) levels are also presented.

Table 3. Optimized Structural Parameters and Spin Populations for Formation of a Product Complex (2P or 2PQM) from a Cu2+−Alkylperoxo Complex (24a or 24aQM) via an Endoperoxide Precursor (2EP or 2EPQM)a parametersb species

Cu−Op

Cu−O4

Op−Od

Op−C4

O(Glu73)···O3

O(Glu)···O4

4a 2 4aQM 2 TS4 2 TS4QM 2 5a 2 5aQM 2 TS5 2 TS5QMc 2 6a 2 6aQM 2 TS6 2 TS6QM 2 EP 2 EPQM 2 TS7 2 TS7QM

1.97 1.93 2.02 1.96 1.98 1.95 1.99 na 1.98 1.93 2.24 2.03 2.54 2.58 2.53 2.54

3.54 3.98 2.86 3.81 2.89 3.72 2.83 na 2.81 3.18 2.41 2.99 1.98 1.98 1.98 1.96

1.41 1.41 1.46 1.45 1.44 1.43 1.44 na 1.43 1.43 1.47 1.47 1.50 1.50 1.87 1.85

3.67 3.91 3.05 3.39 2.59 2.86 2.59 na 2.53 2.66 1.88 1.79 1.49 1.49 1.36 1.36

2.60 2.70 2.56 2.68 2.56 2.65 2.74 na 2.79 2.96 2.89 3.07 2.90 3.01 2.86 3.04

4.89 5.08 4.10 5.29 3.71 5.06 3.00 na 2.76 2.63 2.52 2.47 2.57 2.58 2.57 2.59

2

a

These stationary points were located using the (UB3LYP/def2-SVP)/CHARMM22 and UB3LYP/def2-SVP computations for the QM/MM and QM-only models, respectively. bIn Å. cWe could not locate a saddle point relevant to the hydrogen bond switching.

6959

DOI: 10.1021/acs.jpcb.5b03564 J. Phys. Chem. B 2015, 119, 6952−6962

Article

The Journal of Physical Chemistry B

(13) Matuz, A.; Giorgi, M.; Speier, G.; Kaizer, J. Structural and Functional Comparison of Manganese-, Iron-, Cobalt-, Nickel-, and Copper-Containing Biomimic Quercetinase Models. Polyhedron 2013, 63, 41−49. (14) Sun, Y.-J.; Huang, Q.-Q.; Tano, T.; Itoh, S. Flavonolate Complexes of MII (M = Mn, Fe, Co, Ni, Cu, and Zn). Structural and Functional Models for the ES (Enzyme−Substrate) Complex of Quercetin 2,3-Dioxygenase. Inorg. Chem. 2013, 52, 10936−10948. (15) Barhács, L.; Kaizer, J.; Pap, J.; Speier, G. Kinetics and Mechanism of the Stoichiometric Oxygenation of [CuII(fla)(idpa)]ClO4 [fla = flavonolate, idpa = 3,3′-imino-bis(N,N-dimethylpropylamine)] and the [CuII(fla)(idpa)]ClO4-Catalysed Oxygenation of Flavonol. Inorg. Chim. Acta 2001, 320, 83−91. (16) Balogh-Hergovich, É.; Kaizer, J.; Speier, G.; Huttner, G.; Zsolnai, L. Copper-Mediated Oxygenolysis of Flavonols via Endoperoxide and Dioxetan Intermediates; Synthesis and Oxygenation of [CuII(Phen)2(Fla)]ClO4 and [CuII(L)(Fla)2] [FlaH = Flavonol; L = 1,10-Phenanthroline(Phen), 2,2′-Bipyridine (Bpy), N,N,N′,N′-Tetramethylethylenediamine (TMEDA)] Complexes. Eur. J. Inorg. Chem. 2002, 2287−2295. (17) Czaun, M.; Speier, G.; Kaizer, J.; Bakkali-Taheri, N. El; Farkas, E. Kinetics and Mechanism of the Base-Catalyzed Oxygenation of 1H2Phenyl-3-Hydroxy-4-Oxoquinolines in DMSO/H2O. Tetrahedron 2013, 69, 6666−6672. (18) Hernandez-Ortega, A.; Quesne, M. G.; Bul, S.; Heuts, D. P. H. M.; Steiner, R. A.; Heyes, D. J.; de Visser, S. P.; Scrutton, N. S. Origin of the Proton-Transfter Step in the Cofactor-Free (1H)-3-Hydroxy-4oxoquinaldine 2,4-Dioxygenase. J. Biol. Chem. 2014, 289, 8620−8632. (19) Fiorucci, S.; Golebiowski, J.; Cabrol-Bass, D.; Antonczak, S. Oxygenolysis of Flavonoid Compounds: DFT Description of the Mechanism for Quercetin. ChemPhysChem 2004, 5, 1726−1733. (20) Costas, M.; Mehn, M. P.; Jensen, M. P.; Que, L., Jr. Dioxgen Activation at Mononuclear Nonheme Iron Active Sites: Enzymes, Models, and Intermediates. Chem. Rev. 2004, 104, 939−986. (21) Bugg, T. D. H.; Ramaswamy, S. Non-heme Iron-Dependent Dioxygenase: Unravelling Catalytic Mechanisms for Complex Enzymatic Oxidations. Curr. Top. Chem. Biol. 2008, 12, 134−140. (22) Dong, G.; Shaik, S.; Li, W. Oxygen Activation by Homoprotocatechuate 2,3-Dioxygenase: A QM/MM Study Reveals the Key Intermediate in the Activation Cycle. Chem. Sci. 2013, 4, 3624−3635. (23) Dong, G.; Lai, W. Reaction Mechanism of Homoprotocatecuate 2,3-Dioxygenase with 4-Nitrocatechol: Implication for the Role of Substrate. J. Phys. Chem. B 2014, 118, 1791−1798. (24) Nakatani, N.; Sakaki, S. Multistate CASPT2 Study of Native Iron(III)-Dependent Catachol Dioxygenase and Its Functional Models: Electronic Structure and Ligand-to-Metal Charge Transfer Excitation. J. Phys. Chem. B 2011, 115, 4781−4789. (25) Georgiev, V.; Noack, H.; Borowski, T.; Blomberg, M. R. A.; Siegbahn, P. E. M. DFT Study on the Catalytic Reactivity of a Functional Model Complex for Intradiol-Cleaving Dioxygenase. J. Phys. Chem. B 2010, 114, 5786−5885. (26) Jastrzebsiki, R.; Quesne, M. G.; Weckhuysen, B. M.; de Visser, S. P.; Bruijnincx, P. C. A. Experimental and Computational Evidence for the Mechanism of Intradiol Catechol Dioxygenation by Non-Heme Iron(III). Chem.Eur. J. 2014, 20, 15686−15691. (27) Pau, M. Y. M.; Davis; Pau, M. I.; Orville, A. M.; Lipscomb, J. D.; Solomon, E. I. Spectroscopic and Electronic Structure Study of the Enzyme−Substrate Complex of Intradiol Dioxygenases: Substrate Activation by a High-Spin Ferric Non-heme Iron Site. J. Am. Chem. Soc. 2007, 129, 1944−1958. (28) Pau, M. Y. M.; Lipscomb, J. D.; Solomon, E. I. Substrate Activation for O2 Reactions by Oxidized Metal Centers in Biology. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 18355−18362. (29) Knoot, C. J.; Purpero, V. M.; Lipscomb, J. D. Crystal Structures of Alkylperoxo and Anhydride Intermediates in an Intradiol Ringcleaving Dioxygenase. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 388− 393.

determining step with a barrier of 16.8 kcal/mol. This rearrangement facilitates the following H-bond switching caused by the slight movement of Glu73. Our results also provide better understanding of the roles of flexible Glu73 in the catalytic reaction. The switched H-bond can enhance nucleophilic attack on the target C4 atom to form an endoperoxide intermediate selectively.



ASSOCIATED CONTENT

S Supporting Information *

Natural orbitals relevant to a singlet coupling; comparison between 22b and 22bQM; overlay of 25aQM and 26aQM; full citations for references 43, 52, 55, and 57; optimized Cartesian coordinates of the QM region (in Å). The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b03564.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (T.S.). Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work is partly supported by a grant from Japan Science and Technology Agency (JST) with a Core Research for Evolutional Science and Technology (CREST).



REFERENCES

(1) Oka, T.; Simpson, F. J. Quercetinase, a Dioxygenase Containing Copper. Biochem. Biophys. Res. Commun. 1971, 43, 1−5. (2) Pietta, P.-G. Flavonoids as Antioxidants. J. Nat. Prod. 2000, 63, 1035−1042. (3) Cushnie, T. P. T.; Lamb, A. J. Antimicrobial Activity of Flavonoids. Int. J. Antimicrob. Agents 2005, 26, 343−256. (4) Merkens, H.; Kappl, R.; Jakob, R. P.; Schmid, F. S.; Fentzner, S. Quercetinase QueD of Streptomyces sp. FLA, a Monocupin Dioxygenase with a Preference for Nickel and Cobalt. Biochemistry 2008, 47, 12185−12196. (5) Schaab, M. R.; Barney, B. M.; Francisco, W. A. Kinetic and Spectroscopic Studies on the Quercetin 2, 3-Dioxygenase f rom Bacillus subtilis. Biochemistry 2006, 45, 1009−1016. (6) Solomon, E. I.; Heppner, D. E.; Johnston, E. M.; Ginsbach, J. W.; Cirera, J.; Qayyum, M.; Kieber-Emmons, M. T.; Kjaergaard, C. H.; Hadt, R. G.; Tian, L. Copper Active Sites in Biology. Chem. Rev. 2014, 114, 3659−3853. (7) Fusetti, F.; Schröter, K. H.; Steiner, R. A.; van Noort, P. I.; Pijning, T.; Rozeboom, H. J.; Kalk, K.; Egmond, M.; Dijkstra, B. W. Crystal Structure of the Copper-Containing Quercetin 2,3-Dioxygenase from Aspergillus japonicas. Structure 2002, 10, 259−268. (8) Steiner, R. A.; Kalk, K. H.; Dijkstra, B. W. Anaerobic EnzymeSubstrate Structures Provide Insight into the Reaction Mechanism of the Copper-Dependent Quercetin 2,3-Dioxygenase. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 16225−16230. (9) Steiner, R. A. In Encyclopedia of Inorganic and Biorganic Chemistry; Scott, R. A., Ed.; John Wiley & Sons, Ltd.: Chichester, 2013; pp 1−12. (10) Balogh-Hergovich, É.; Kaizer, J.; Speier, G. Kinetics and Mechanism of the Cu(I) and Cu(II) Flavonolate-Catalyzed Oxygenation of Flavonol. Functional Quercetin 2,3-Dioxygenase Models. J. Mol. Catal. A 2000, 159, 215−224. (11) Pap, J. S.; Kaizer, J.; Speier, G. Model Systems for the COReleasing Flavonol 2, 4-Dioxygenase Enzyme. Coord. Chem. Rev. 2010, 254, 781−793. (12) Kaizer, J.; Pap, J. S.; Speier, G. In Copper-Oxygen Chemistry; Karlin, K. D. Itoh, S., Eds.; John Wiley & Sons, Inc.: Hoboken, NJ, 2011; pp 23−52. 6960

DOI: 10.1021/acs.jpcb.5b03564 J. Phys. Chem. B 2015, 119, 6952−6962

Article

The Journal of Physical Chemistry B (30) Siegbahn, P. E. M. Hybrid DFT Study of the Mechanism of Quercetin 2,3-Dioxygenase. Inog. Chem. 2004, 43, 5944−5953. (31) Antonczak, S.; Fiorucci, S.; Golebiowski, J.; Cabrol-Bass, D. Theoretical Investigations of the Role Played by Quercetinase Enzymes upon the Flavonoids Oxygenolysis Mechanism. Phys. Chem. Chem. Phys. 2009, 11, 1491−1501. (32) Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme. J. Mol. Biol. 1976, 103, 227−249. (33) Field, M. J.; Bash, P. A.; Karplus, M. A Combined Quantum Mechanical and Molecular Mechanical Potential for Molecular Dynamics Simulations. J. Comput. Chem. 1990, 11, 700−733. (34) Lin, H.; Truhlar, D. G. QM/MM: What Have We Learned, Where Are We, and Where Do We Go from Here? Theor. Chem. Acc. 2007, 117, 185−199. (35) Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem., Int. Ed. 2009, 48, 1198−1229. (36) Lonsdale, R.; Harvey, J. N.; Mulholland, A. J. A Practical Guide to Modelling Enzyme-Catalysed Reactions. Chem. Soc. Rev. 2012, 41, 3025−3038. (37) de Visser, S. P.; Quesne, M. G.; Martin, B.; Comba, P.; Ryde, U. Computational Modeling of Oxygenation Processes in Enzymes and Biomimetic Model Complexes. Chem. Commun. 2014, 50, 262−282. (38) Sali, A.; Blundell, T. L. Comparative Protein Modelling by Satisfaction of Spatial Restraints. J. Mol. Biol. 1993, 234, 779−815. (39) Olsson, M. H. M.; Søndergaard, C. R.; Rostkowski, M.; Jensen, J. H. PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa Predictions. J. Chem. Theory Comput. 2011, 7, 525−537. (40) Cox, S. R.; Williams, D. E. Representation of the Molecular Electrostatic Potential by a Net Atomic Charge Model. J. Comput. Chem. 1981, 2, 304−323. (41) Ungar, L. W.; Scherer, N. F.; Voth, G. A. Classical Molecular Dynamics Simulation of the Photoinduced Electron Transfer Dynamics of Plastocyanin. Biophys. J. 1997, 72, 5−17. (42) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4, 187−217. (43) Brooks, B. R.; Brooks, C. L., III; Mackerell, A. D., Jr.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545−1614. (44) Brooks, B. R.; Karplus, M. Deformable Stochastic Boundaries in Molecular Dynamics. J. Chem. Phys. 1983, 79, 6312−6325. (45) Becke, A. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098−3100. (46) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785−789. (47) Becke, A. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 93, 5648−5652. (48) Schäfer, A.; Horn, H.; Ahlrichs, R. Fully Optimized Contracted Gaussian Basis Sets for Atoms Li to Kr. J. Chem. Phys. 1992, 97, 2571− 2577. (49) Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (50) Yamaguchi, K.; Takahara, Y.; Fueno, T. In Applied Quantum Chemistry; Smith, V. H.; Schaefer, H. F., III, Morokuma, K., Eds.; D. Reidel: Boston, MA, 1986; p 15. (51) Yamaguchi, K.; Jensen, F.; Dorigo, A.; Houk, K. N. A Spin Correction Procedure for Unrestricted Hartree-Fock and MøllerPlesset Wavefunctions for Singlet Diradicals and Polyradicals. Chem. Phys. Lett. 1988, 149, 537−542. (52) Sherwood, P.; de Vries, A. H.; Guest, M. F.; Schreckenbach, G.; Catlow, C. R. C.; French, S. A.; Sokol, A. A.; Bromley, S. T.; Thiel, W.;

Turner, A. J.; et al. QUASI: A General Purpose Implementation of the QM/MM Approach and its Application to Problems in Catalysis. J. Mol. Struct. (THEOCHEM) 2003, 632, 1−28. (53) Chemshell, a Computational Chemistry Shell, version 3.5; Science & Technology Facilities Council: Swindon, UK; www.chemshell.org (accessed on April 13, 2015). (54) Bakowies, D.; Thiel, W. Hybrid Models for Combined Quantum Mechanical and Molecular Mechanical Approaches. J. Phys. Chem. 1996, 100, 10580−10594. (55) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian 09, Revision C.01; Gaussian, Inc.: Wallingford CT, 2009. (56) Smith, W.; Forester, T. R. DL_POLY_2.0: A General-Purpose Parallel Molecular Dynamics Simulation Package. J. Mol. Graphics 1996, 14, 136−141. (57) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (58) Okamoto, T.; Yamada, K.; Koyano, Y.; Asada, T.; Koga, N.; Nagaoka, M. A Minimal Implementation of the AMBER−GAUSSIAN Interface for Ab Initio QM/MM-MD Simulation. J. Comput. Chem. 2011, 32, 932−942. (59) Saito, T.; Thiel, W. Quantum Mechanics/Molecular Mechanics Study of Oxygen Binding in Hemocyanin. J. Phys. Chem. B 2014, 118, 5034−5043. (60) Kästner, J.; Carr, J. M.; Keal, T. W.; Thiel, W. Wander, A.; Sherwood, P. DL-FIND: An Open-Source Geometry Optimizer for Atomistic Simulations. J. Phys. Chem. A 2009, 113, 11856−11865. (61) Nocedal, J. Updating Quasi-Newton Matrices with Limited Storage. Math. Comput. 1980, 35, 773−782. (62) Henkelman, G.; Jónsson, H. Improved Tangent Estimate in the Nudged Elastic Band Method for Finding Minimum Energy Paths and Saddle Points. J. Chem. Phys. 1999, 111, 7010−7022. (63) Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113, 9901−9904. (64) Siegbahn, P. E. M. The Performance of Hybrid DFT for Mechanisms Involving Transition Metal Complexes in Enzymes. J. Biol. Inorg. Chem. 2006, 11, 695−701. (65) Himo, F. Quantum Chemical Modeling of Enzyme Active Sites and Reaction Mechanisms. Theor. Chem. Acc. 2006, 116, 232−240. (66) Kearsley, S. K. On the Orthogonal Transformation Used for Structural Comparisons. Acta Crystallogr. 1989, A45, 208−210. (67) Yoshioka, Y.; Yamaguchi, K.; Fueno, T. Heisenberg Models of Radical Reactions II. Conservation of the Local Spin-Permutation Symmetry in Reactions of Biradical Species. Theor. Chim. Acta 1977, 45, 1−20. (68) Shaik, S.; Cohen, S.; Wang, Y.; Chen, H.; Kumar, D.; Thiel, W. P450 Enzymes: Their Structure, Reactivity, and Selectivity-Modeled by QM/MM Calculations. Chem. Rev. 2010, 110, 949−1017. (69) Shaik, S.; Chen, H. Lessons on O2 and NO Bonding to Heme from Ab Inito Multireference/Multiconfiguration and DFT Calculations. J. Biol. Inorg. Chem. 2011, 16, 841−855. (70) In the previous QM model, the ligating histidines and glutamic acid were modeled by imidazoles and a formic acid. The hydroxyl groups of A- and C-rings were also replaced by hydrogen atoms. The LACVP (LANL2DZ for copper and 6-31G for other atoms) for geometry optimizations. (71) Lonsdale, R.; Harvey, J. N.; Mulholland, A. J. Inclusion of Dispersion Effects Significantly Improves Accuracy of Calculated Reaction Barriers for Cytochrome P450 Catalyzed Reactions. J. Phys. Chem. Lett. 2010, 1, 3232−3237. (72) Lonsdale, R.; Harvey, J. N.; Mulholland, A. J. Effects of Dispersion in Density Functional Based Quantum Mechanical/ Molecular Mechanical Calculations on Cytochrome P450 Catalyzed Reactions. J. Chem. Theory Comput. 2012, 8, 4637−4645. 6961

DOI: 10.1021/acs.jpcb.5b03564 J. Phys. Chem. B 2015, 119, 6952−6962

Article

The Journal of Physical Chemistry B (73) Saito, T.; Kataoka, Y.; Nakanishi, Y.; Matsui, T.; Kitagawa, Y.; Kawakami, T.; Okumura, M.; Yamaguchi, K. Which Hybrid GGA DFT is Suitable for Cu2O2 Systems If the Spin Contamination Error is Removed? Chem. Phys. 2010, 368, 1−6.

6962

DOI: 10.1021/acs.jpcb.5b03564 J. Phys. Chem. B 2015, 119, 6952−6962