MM Study of the Conversion of Oxophlorin into Verdoheme by

Nov 1, 2017 - The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b08332. Relative energi...
0 downloads 0 Views 1MB Size
Subscriber access provided by READING UNIV

Article

QM/MM Study of the Conversion of Oxophlorin Into Verdoheme by Heme Oxygenase Fatemeh Sadat Alavi, Mansour Zahedi, Nasser Safari, and Ulf Ryde J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b08332 • Publication Date (Web): 01 Nov 2017 Downloaded from http://pubs.acs.org on November 7, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

215x279mm (144 x 144 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

QM/MM Study of the Conversion of Oxophlorin into Verdoheme by Heme Oxygenase Fatemeh Sadat Alavia, Mansour Zahedia*, Nasser Safaria, Ulf Rydeb a

Department of Chemistry, Faculty of Sciences, Shahid Beheshti University, G.C., Evin,

19839-6313, Tehran, Iran b

Department of Theoretical Chemistry, Lund University, Chemical Centre, P.O. Box 124,

SE-221 00 Lund, Sweden

Abstract Heme oxygenase is an enzyme that degrades heme, thereby recycling iron in most organisms, including human. Pervious density functional theory (DFT) calculations have suggested that iron(ІІI) hydroxyheme, an intermediate generated in the first step of heme degradation by heme oxygenase, is converted to iron(ІІІ) superoxo oxophlorin in the presence of dioxygen. In this paper, we have studied the detailed mechanism of conversion of iron(ІІІ) superoxo oxophlorin to verdoheme by using combined quantum mechanics and molecular mechanics (QM/MM) calculations. The calculations employed the B3LYP method and the def2-QZVP basis set, considering dispersion effects with the DFT-D3 approach, obtaining accurate energies with large QM regions of almost 1000 atoms. The reaction was found to be exothermic by –35 kcal/mol, with a rate-determining barrier of 19 kcal/mol in the doublet state. The protein environment and especially water in the enzyme pocket significantly affects the reaction by decreasing the reaction activation energies and changing the structures by providing strategic hydrogen bonds.

*

Corresponding author

ACS Paragon Plus Environment

Page 2 of 23

Page 3 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Introduction Heme oxygenase (HO) is an enzyme found in many organisms, including bacteria, plants, and vertebrates. It is important in human physiology, degrading heme and recovering iron.1 In fact, it has been shown that only 3% of the human iron requirement is fulfilled via food intake, whereas the remainder is obtained by the action of HO.2 Consequently, much effort has been spent to study the function of this enzyme. Three different isoforms of HO are known.3, 4 The first type, HO-1, is induced by oxidative stress, whereas the second, HO-2 is a constitutive isozyme.5,

6

Both isoforms degrade heme,

whereas HO-3, the third type, has no catalytic activity and seems to act as an oxygen sensor.7 HO-1 has a molecular weight of 32 kDa and is a member of the stress protein superfamily (HSP32). It has a broad spectrum of inducers and is found in spleen, liver, and bone marrow.8 On the other hand, the 36 kDa HO-2 does not respond to factors that induce HO-19 but is observed in high levels in testis, brain, and endothelial cells of cerebral vessels.10 Expression of HO-1 provides antioxidant protection in a variety of cells and tissues. In fact, it has been suggested that modulation of HO-1 activity could have a potential therapeutic value.11 Researchers have paid much attention to HO-1 for a long time.12 In contrast, HO-2 has been less studied, owing to its apparent constitutive role. Both HO-1 and HO-2 oxidize heme to biliverdin ІΧα in a multistep reaction mechanism, requiring O2 and reducing equivalents from NADPH cytochrome P450 reductase.13-16 The first step of heme degradation catalyzed by HO is the oxidation of heme to α-meso-hydroxyheme.17 The second step is the formation of verdoheme with concomitant release of the hydroxylated αmeso carbon as carbon monoxide.18 In the third step, biliverdin is formed from verdoheme, concomitantly with release of Fe2+19 (see Scheme 1). Interestingly, all the products of the HO reaction are biologically active species: iron is a gene inducer, biliverdin is an antioxidant, and carbon monoxide is a neurotransmitter, anti-inflammatory, vasodilator, and a promoter of neovascular growth.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Scheme 1. The three-step degradation of heme catalyzed by HO.

During the latest decades, computational chemistry has been established as a useful complement to experimental work for the study of the structure and function of many important biological processes.20, 21 Fe-containing proteins, being one of the most important enzymes, have been the subject of many molecular modeling studies.22-25 In particular, quantum mechanical (QM) calculations have been shown to be very powerful to identify intermediates and transition states in enzyme reactions.26, 27 Unfortunately, QM calculations are too demanding to study the reactions and dynamics of an entire protein; due to restrict to small models of the active site. By the QM/MM approach, the reactive system (20–200 atoms) is treated by accurate QM methods, whereas the influence of the surroundings is taken into account with MM calculations. Using such methods, it is possible to simulate the entire protein and the surrounding solvent with a good accuracy. In this work, we have studied the heme cleavage mechanism of HO with QM/MM methods. Previous computational and experimental studies, have concentrated on the first and the third steps of the HO mechanism in Scheme 1. However, the second step has attracted less attention and several features of this step are still controversial. In this step, Fe(III) hydroxyheme is converted to Fe(II) verdoheme and CO, consuming an O2 molecule and one electron (Scheme 1). The CO molecule comes from the α-meso carbon and the attached hydroxyl group. The ring is re-sealed with an oxygen atom from O2.28 This step is special in that all four meso-hydroxyheme isomers (α, β, γ, and δ, cf. Scheme 1) can be converted to verdoheme, in contrast to the other two steps, in which only the α isomer is active.29 In fact, it has been observed that hydroxyheme reacts with oxygen and produces verdoheme even in the absence of HO.30 This indicates that the second step of the HO mechanism can be considered as a spontaneous auto-oxidation of hydroxyheme. Moreover, only this step of the HO reaction is not inhibited by

ACS Paragon Plus Environment

Page 4 of 23

Page 5 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

CO ligation to the iron and it cannot be supported with H2O2 instead of O2.11, 31 Several investigations of the conversion of α-meso hydroxyheme to verdoheme have been performed by employing chemical model systems such as α-meso hydroxyheme in HO-132 or myoglobin.33 These studies have indicated that radical species are formed. However, the reduction equivalent necessary for the reaction have been the subject of a controversy. 31 In a previous QM-cluster study, our group studied this part of the HO reaction, suggesting a threestep mechanism, the reaction is exothermic by –15 kcal/mol and the rate-limiting barrier is 12 kcal/mol. Moreover, the nature of the intermediates and the reaction mechanism for these transformations are still not clear. Here the conversion of Fe(III) hydroxyheme to verdoheme in HO has been studied. For such multistep process which includes carbon-oxygen bond formation, O-O bond cleavage and further rearrangement and loss of CO moiety; no information is available and its mechanism is totally unknown. Due to the above considerations and in order to provide an answer for the fact that what influences is imparted to oxophlorin reactivity by the surrounding environment and in continuation of previous studies

34-36

, study of the enzymatic

reaction environment path has been considered necessary. In this work, we have studied the heme cleavage mechanism of HO with QM/MM methods.

Methods The QM/MM optimizations were performed with the COMQUM software.37-39 It uses the Turbomole software40 for the QM part and Amber41 for MM part. In this approach, the protein and solvent were divided into two regions: System 1 was optimized at the QM level, system 2 was kept fixed at the starting structure and was treated by MM. In the QM calculations, system 1 was represented by a wave function, whereas all of the other atoms were modeled by partial point charges, one for each atom, taken from MM libraries.42 Thereby, the polarization of the QM system by the surrounding is included in a self-consistent manner (electrostatic embedding). Covalent bonds between systems 1 and 2 were treated by the hydrogen link-atom approach:43 The QM system was truncated by hydrogen atoms, the positions of which are related to the corresponding carbon atoms in the protein.37

The total QM/MM energy was calculated as

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

EQM/MM = EQM + EMM12 – EMM1

(1)

where EQM is the QM energy of the quantum system truncated by hydrogen atoms, including a point-charge model of the surrounding systems 2. EMM1 is the MM energy of the quantum system, still truncated by hydrogen atoms, but without any electrostatic interactions. Finally, EMM12 is the classical energy of all atoms with carbon atoms at the junctions and with the charges of the QM region zeroed, to avoid double-counting of the electrostatic interactions. This approach is similar to the one used in ONIOM method.44 The QM calculations were carried out using the TPSS functional45 with the Turbomole package.40 The def2-SV(P) basis set46, 47 was used for all atoms in the geometry optimizations. Dispersion effects were treated by the DFT-D3 approach and Becke Johnson damping.48 More accurate energies were obtained with single-point energy calculations with the def2-TZVP and def2-QZVP basis sets49 (still with TPSS-D3), and with the B3LYP-D3 method50, 51 and the def2TZVP basis set. All single-point calculations were done with a point-charge model of the surroundings. All QM calculations were sped up by expanding the Coulomb interactions in an auxiliary basis set, the resolution-of-identity approximation.49, 52 It is well-known that QM calculations are insensitive to the method used for the geometry calculations.53 Therefore, most QM studies employ a cheap method (e.g. a medium-sized basis set) for geometry optimizations and then a better QM method for single-point energy calculations. We have followed this approach in the present study, optimizing geometries with the TPSS-D3/def2-SV(P) approach and then calculating energies with both TPSS-D3 and B3LYP-D3 with the much larger def2-TZVPD basis set. Owing to the resolution-ofidentity approach, a pure DFT functional is ~10 times faster than a hybrid functional and they often give more accurate metal-ligand bonds than B3LYP (which tends to overestimate such bonds).54 The single-point energy calculations with large basis sets give much improved energies and they also allow us to judge whether the energies are sensitive to the amount of exchange in the functional. This approach has successfully been used for many different bioinorganic systems;55, 56 and has sometimes pointed out unexpectedly large errors in previous calculations.57 The MM calculations were run with the sander module of the Amber software41 using the Amber ff14SB force field58 for the protein and the general AMBER force field (GAFF)59 for heme, O2. The active site was represented by charges fitted to the electrostatic potential (ESP), calculated in the following way: The structure was optimized at DFT level with TPSS functional

ACS Paragon Plus Environment

Page 6 of 23

Page 7 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

and def2-SV(P) basis set. Then, the ESP was calculated with single-point calculation, at positions selected according to the Merz–Kollman scheme.60 Finally, charges were fitted to the ESP at these points.

Modeled system The calculation in this study were based on the X-ray structure of human HO-1 in complex with heme at 1.5 Å resolutions (PDB entry 1N45).61 The high resolution of this structure revealed important details of the active site, such as a network of hydrogen-bonded water molecules. The protein is dimeric in the crystal, but we included only chain A in the calculations. All crystal-water molecules around this subunit were kept in the calculations. The heme substrate was manually converted to α-meso hydroxyheme and O2 was added as a sixth ligand of the Fe ion. All residues were assumed to be in their most stable protonation state at neutral pH. Thus, all Arg and Lys residues were protonated and positively charged, whereas all Asp and Glu residues, as well as the two propionate groups of hydroxyheme, were deprotonated and negatively charged. After a detailed study of the surroundings, solvent accessibility and possible hydrogenbond networks around the His residues, it was assumed that His25 is protonated on Nδ1 atom, His56, 132 and 223 are protonated on the Nɛ2 atom, and the remaining His residues (all solventexposed) are doubly protonated. The protonation states of all residues were checked with the PROPKA software62. This assignment gave a total charge of –2. According to many experimental investigations,63 the distal acid–base pair Arg127–Asp131 and several crystal water molecules form a cluster in the distal pocket, which is crucial for the catalytic function of HO. Therefore, our QM region comprised iron porphyrin (including all side chains) with its proximal ligand His25 (modeled as methyl imidazole), the distal residues Arg127 (modeled as CH3NHC(NH2)2+) and Asp131 (modeled as CH3COO–) and six crystal waters. This QM system is shown in Figure 1. The remaining parts of protein together with the surrounding water molecules were included in system 2. Consequently, the simulated system contained 24157 atoms with 126 atoms in the QM region.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 1. The QM region used in the calculations. Scheme 2. Atom numbering used in this study

ACS Paragon Plus Environment

Page 8 of 23

Page 9 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Big-QM calculations QM/MM calculations show a fairly slow convergence with respect to the size of the QM system.64 The big-QM approach is a method to post-process QM/MM calculations to get more stable energies.65,

66

. It involves calculations with very large QM systems, including atoms

according to the following three rules: 

All atoms within 6 Å of a minimal QM system involving the porphyrin ring, His25 and the distal ligand.



All buried charged groups in the protein, which for HO are Arg127 and Asp131.



Junctions are moved two residues away from the minimal QM system.

This gave a QM region of almost 1000 atoms. The calculations employed the TPSS method with the def2-SV(P) basis set. This energy was enhanced with a DFT-D3 dispersion correction, calculated for the same big QM region with Becke Johnson damping,48 third-order terms, and default parameters for the TPSS functional using dftd3 program.67 We also included a standard MM correction for this big-QM system:

Ebq= Ebig-QM + EMM12 – EMM1 + Edisp (2) Finally, the energies were extrapolated to the B3LYP-D3/def2-QZVP level of theory by four single-point energy calculations for the small QM system, including a point-charge model of the surroundings: Etot = Ebq + ETPSS/QZ – ETPSS/SV + EB3LYP/TZ – ETPSS/TZ

(3)

The B3LYP method was selected because it has been shown to give reliable energies of heme systems. The calculations on the big-QM systems employed the multipole-accelerated resolution-of-identity J approach.68

Result and Discussion Previous studies34 have shown that hydroxyheme (FeIII(POH)), which is the product of the first step of heme degradation by HO, can be oxidized in the presence of O2 to give Fe(ІІІ)–

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 23

superoxo–oxophlorin (FeIII(PO•)(O•2)). This deprotonated radical species is essential for the subsequent reaction with O2 to yield verdoheme and to cleave the porphyrin ring because it is highly reactive toward oxygen.35 In this paper, we study the conversion of this intermediate to verdoheme. We have used QM/MM calculations with the entire oxophlorin molecule (including all side chains), His25 as an axial ligand, the distal residues Arg127 and Asp131, as well as six waters in the QM region. Starting from Fe(ІІІ)–superoxo–oxophlorin (FeIII(PO•)(O•2)) (1), the terminal oxygen atom in O2 attacks the carbon atom adjacent to the keto group of oxophlorin C7 (Scheme 3; atom numbers are shown in Scheme 2). In transition state (TS1), the distance between OB and C7 is shortened to 1.80 Å, while the OA–OB bond is extended to 1.38 Å. In the resulting intermediate 2, O2 bridges between the iron ion and the ring. Next, the O–O bond is cleaved homolytically, which leads to intermediate 3. In this, the Fe–OA bond length is 1.64 Å, indicating that it contains a FeIV=O group. The bond between C7 and OB is reduced to 1.36 Å. Subsequently, the OB atom attacks the other carbon atom (C8) adjacent to the keto group, giving rise to intermediate 4, in which the C6– C8 bond is cleaved. In the next step, the C6 and C7 bond is homolytically broken, which leads to the release of a CO molecule and the formation of a bond between OB and C8, giving rise to Fe(IV)-oxo verdoheme (5). Finally, an internal electron transfer from the ring to Fe(IV) produces the Fe(III)-oxo oxidation state and a closed-shell verdoheme macrocycle that after double protonation of the oxo group and dissociation of water molecule gives rise to Fe(II) in verdoheme. Optimized QM/MM structures of all intermediates in this reaction mechanism are shown in Figure 2. It can be seen that the water cluster forms strong hydrogen bonds for all intermediates (shown as green dotted lines in Figure 2). In particular, there is a strong hydrogen bond between OB and water in the intermediate 3, which facilitates the formation of the bond between OB and C8. Initially, we saw some reorientation of the hydrogen bonds in the water cluster during the reaction. However, by running the reaction forth and back, this turned out to represent only different local minima. In the structures shown in Figure 2, the water cluster shows the same hydrogen-bond pattern throughout the whole reaction mechanism. The only notable difference is that there are hydrogen bonds only to OB in reactant 1, but this shifts to OA in the other intermediates. However, new hydrogen bonds are established to OB in the crucial intermediate 3.

ACS Paragon Plus Environment

Page 11 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1

2

3

4

5

6

Fig. 2. QM/MM optimized structures of the various intermediates in the suggested mechanism.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 23

Scheme 3. The suggested mechanism for the heme degradation from Fe(ІІІ)superoxo oxophlorin to verdoheme.

The reactant 1 has a complicated electronic structure, as has been discussed previously. 36,35 Superoxide has one unpaired electron, the doubly anionic oxophlorin radical a second unpaired electron and Fe(III) one, three or five unpaired electrons, depending on the spin state. All these unpaired electrons can be either parallel or partly antiparallel. We have studied the three lowest spin states and the calculated QM/MM energies are reported in Table 1. The ground state was found to be the doublet with the electronic structure

(the first three are Fe 3d

orbitals, the other two being unoccupied, and the last two are orbitals on O2 and on oxophlorin, respectively), with a low-spin state on Fe(III) and with the unpaired electrons antiferromagnetically coupled. The corresponding quartet state with the three spins parallel was found to be only 0.5 kcal/mol less stable and it is expected to give a similar reactivity. The sextet state (

), with an intermediate-spin state on Fe(III) and with the unpaired

electrons ferromagnetically coupled, is 31 kcal/mol less stable than the doublet state.

ACS Paragon Plus Environment

Page 13 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Key geometric parameters of the optimized structure of 1 in the three spin states are also collected in Table 1. The Fe–O distance is much longer in sextet state (2.50 Å) than in the doublet and quartet states (1.84–1.90 Å), in line with the Fe–O and Fe–NIm antibonding character of the Fe

orbital, which is singly occupied in the sextet state, but vacant in the doublet and

quartet states. As an effect, the O–O bond in the sextet state is considerably shorter than in the other two spin states (1.23 Å, compared to 1.27–1.28 Å; cf. Table 1). Consequently, it will be harder to cleave. Therefore, the sextet surface is unlikely to play any role in the reaction mechanism. On the other hand, we have optimized both the doublet and quartet spin states for all species with the TPSS/def2-SV(P) method. The energies are shown in Table S2 in the Supporting Information. For some species, the quartet state is only 0.1–0.6 kcal/mol higher in energy than the doublet state. However, for the first species like, TS1, 2 and TS2, as well as the product (6), the energy difference is much higher, 11–16 kcal/mol, showing that the quartet spin state does not play a main role in mechanism. Interestingly, the rate-limiting step is different for the two spin states. In the quartet spin state, the OA–OB distance is shorter than in the doublet state. Therefore, the O–O bond cleavage (second reaction step) requires more energy than the OB–C7 bond formation (first step) in the quartet state. Consequently, the second step of the mechanism (TS2) is rate-limiting for the quartet spin state, whereas the first step (TS1) is rate-limiting for the doublet state. Table 1. Key Bond Distances (in Å) of the Three Spin States of Iron (ІІІ)–superoxo–oxophlorin Fe–NIm Fe–OA OA–OB Fe–NPO (av) Energy 2

1

2.00

1.84

1.27

2.01

0.0

4

1

2.05

1.90

1.28

2.02

0.5

6

1

2.10

2.50

1.23

2.02

30.9

Key bond distances of the various intermediates and transition states in the doublet state are listed in Table 2. It can be seen that the Fe–OA bond length is 1.84–1.90 Å in the first three states of the reaction, when it represents a Fe–O single bond, but it then decreases to around 1.64 Å in the other states, reflecting the change to a FeIV=O double bond.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 23

Table 2. Key Bond Distance of the various intermediates and transition states in the reaction mechanism, calculated for the doublet state Distance OA–OB Fe–OA OB–C8 OB–C7 C6–C7 C6–C8 OA–HWa OB–HWa

a

1

1.28

1.84

3.27

3.50

1.46

1.46

2.52

1.77

TS1

1.38

1.90

2.70

1.80

1.51

1.45

2.11

2.60

2

1.52

1.85

2.79

1.49

1.53

1.45

1.87

2.50

TS2

1.60

1.77

2.81

1.46

1.54

1.56

1.95

2.72

3

3.34

1.64

2.48

1.36

1.54

1.49

1.84

2.20

TS3

3.35

1.64

1.90

1.40

1.54

1.51

1.80

2.80

4

3.33

1.65

1.46

1.47

1.56

1.56

1.77

3.00

TS4

3.40

1.65

1.44

1.43

1.80

1.57

1.76

2.80

5

3.77

1.64

1.35

1.35

3.60

2.74

1.82

3.30

Shortest O–H distance to a water molecule of the cluster in the distal pocket

Fig. 3. Relative Energies throughout the reaction mechanism, calculated with QM/MM (Etot energies, cf. Eqn. 3).

ACS Paragon Plus Environment

Page 15 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Energies were calculated with three different basis sets, def2-SV(P), def2-QZVP and def2TZVP and two different methods, TPSS and B3LYP. The various results are listed in the SI as Table S1 and depicted as Figure S1 and they are combined to Etot according to Eqn. 3. The relative Etot energies for various states in the reaction mechanism are shown in Figure 3. It can be seen that the first step is rate limiting with an activation energy of 19 kcal/mol. In a previous QM-cluster study,34 our group studied this part of the HO reaction, suggesting a three-step mechanism. The first step was similar to the one in the current mechanism and it was also rate limiting, by 12 kcal/mol. However, in the later stages, the two mechanisms differ in that the formation of the bond between OB and C8 and the release of carbon monoxide, took place in separate steps in the QM/MM calculations, whereas these changes occurred simultaneously in the QM-cluster calculations. Moreover, the final product is notably more stable in the QM/MM calculations (ΔEtot= –35 kcal/mol) than in the QM-cluster calculations (ΔEtot = –15 kcal/mol). It is interesting to compare the structures obtained in the two studies. In the QM-cluster calculations, the Fe–O–O moiety and the imidazole group of His25 were parallel (the O A–Fe– N3–C4 dihedral angle was 1º) for the reactant state 1. The calculations in gas phase showed that this conformation was 1 kcal/mol more stable than the orthogonal conformation. The orientation of the histidine ligand causes only minor differences in the relative energies but according to pervious article,34 it has a great impact on the geometry of the macrocycle. However, during the reaction, the conformation changed so that for the product state, the orthogonal orientation was preferred. In the current protein environment, the orientation of His25 is constrained to the orthogonal orientation by the protein environment, in particular by a hydrogen bond to Glu29. In the QM-cluster calculations, only the bare oxophlorin ring without any side chains, and with no protein residues, except the axial His25 ligand, was considered. Instead, the surroundings were implicitly treated by continuum-solvation methods. In contrast, the QM/MM calculations included iron-oxophlorin with all side chains, O2, His25, as well as the distal residues Arg127 and Asp131 and six crystal waters in QM region and the remainder of the protein and some solvent in the MM region. In the big-QM calculations, all atoms within 6 Å of the minimal QM-cluster system were considered by QM (almost 1000 atoms) and the remainder of the protein (almost 23000 atoms), as well as the surrounding solvent were included as point charges. Undoubtedly, this should give more realistic results.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 23

The Etot energies contain several contributions as can be seen in Eqns. 2–3. They are based on the big-QM (TPSS/def2-SV(P)) energies. To these, dispersion and MM corrections are added (Eqn. 2). The latter term is minimal (less than 0.7 kcal/mol), whereas the former is somewhat larger (up to 3 kcal/mol), as can be seen from Table S1 in the supporting information. Then, the energy is corrected to the B3LYP-D3/def2-QZVP level of theory by a number of single-point calculations. From the raw data in Table S1, it can be seen that the basis-set correction is large, 1–7 kcal/mol, typically favoring the intermediates and transition states. In fact, already the def2TZVP basis set is almost saturated – relative energies calculated with this state differs by less than 0.6 kcal/mol from those with the def2-QZVP basis, Calculations with the B3LYP functional also change the relative energies substantially (1–7 kcal), but normally with an opposite sign, compared to the basis set (i.e. increasing the energies of the intermediates and transition states). To have a better insight of how the surrounding protein affect the reaction, we have performed a series of single-point calculations on the QM/MM structures for the various states in the reaction mechanism (all at the TPSS-D3/def2-SV(P) level of theory). First, we considered the various terms in Eqn. 1. The EMM = EMM12 – EMM1 energy is a MM energy correction to the QM energy, representing primarily the van der Waals interactions of the QM region with the surroundings for the various states in the reaction. From Table 3, it can be seen that EMM is rather small for all states, 0–4 kcal/mol. It is always positive, indicating that this term favors the starting state and increases all reaction and activation energies. Next, we recalculated all QM energies without the point-charge model. The difference between calculations with and without the point charges show the electrostatic effect from the surroundings. From Table 3 (column ptch), it can be seen that this effect is large, up to 13 kcal/mol. It decreases all reaction and activation energies, except for the intermediate 3. In particular, it decreases the rate-limiting barrier by 6 kcal/mol. However, the effect is even larger for the three last states of the reaction mechanism, which are stabilized by 11–13 kcal/mol. We have used a quite large QM region in the QM/MM calculations, with 126 atoms, including the porphyrin side chains, Arg127, Asp131 and six water molecules. Therefore, the effects of these groups are not included in the electrostatic stabilization, discussed in the previous paragraph. To estimate these effects, QM energies (without the point-charge model) were recalculated with additional QM regions in which we first replaced the porphyrin side chains with hydrogen atoms, then the models of Arg127 and Asp131 were deleted, and finally the six

ACS Paragon Plus Environment

Page 17 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

water molecules were omitted. These calculations gave three new reaction profiles that directly show the effect of these groups on the reaction energies. These results are also included in Table 3. It can be seen that the porphyrin side chains have a strong effect (up to 14 kcal/mol) on the reaction energies, essentially counteracting the electrostatic effect of the surroundings (within 3 kcal/mol, except for intermediate 3), reflecting the large effect of two net charges, which are compensated by the surrounding protein. On the other hand, Arg127 and Asp131, had a smaller effect, up to 3 kcal/mol, again except for intermediate 3 (8 kcal/mol). They increase all energies, except for the rate-limiting step (decreased by 1 kcal/mol). The water cluster has a somewhat larger effect, up to 7 kcal/mol. It stabilizes all intermediates and transition states. When all these groups are removed from the QM region, we have arrived at the QM model used in the previous QM-cluster calculations.34 The remaining energy difference between the two sets of calculations represent the effect of geometric constraints in the protein (i.e. the effect of optimizing the various states in the protein or in a vacuum), and the effect of the continuumsolvation model. From the last column of Table 3, it can be seen that the geometric constraints of the protein have a rather small effect on several states (1, TS1,and 3). However, states 2 and TS2 are actually more stable in the vacuum geometry than in the protein by 10–12 kcal/mol. On the other hand, the two last states are stabilized in the protein by 5–10 kcal/mol. If we sum the first five effects in Table 3, it can be seen that the total effect of the surroundings is rather small (up to 2 kcal/mol) for all states except intermediate 3 and the following TS3, which are stabilized by 10 and 5 kcal/mol, respectively. However, these conclusions are obtained at the QM/MM level. The more accurate big-QM calculations indicate that all reaction and activation barriers are further decreased by up to 5 kcal/mol when the larger QM system is employed (difference by the QM/MM and total big-QM energies in Table S1). In particular, the rate-limiting barrier is decreased by 4.9 kcal/mol. Unfortunately, it is harder to obtain reliable contributions to these energies. The previous QM-cluster calculation reported that the rate-limiting activation barrier decreased by 1.6 kcal/mol when recalculated in a continuum solvent with a dielectric constant of 5.7. The present QM/MM calculations indicate that the effect can be appreciably larger and that a detailed account of the surroundings is important to get accurate reaction energies of protein reactions.

Table 3. Effects of the surrounding protein for the various states in the reaction mechanism. The

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 23

first column shows the effect of the EMM = EMM12 – EMM1 energy term in Eqn. 1. The second column (ptch) is the difference between the QM energy with and without point charges. The next three terms represent the effects when first the porphyrin side chains are replaced by hydrogen atoms, then the model of the Arg127 and Asp131 models are removed from the QM region, and finally the effect when also the six water molecules are removed from the QM region. The last column shows the effect of the geometry, i.e. the difference in energy of the minimal QM system (Fe, oxophlorin without side chains, O2, and His25), when optimized with QM/MM or in the previous QM-cluster calculations in a vacuum.34 All energies were obtained by single-point TPSS-D3/def2-SV(P) calculations. MM

ptch sidech Arg+Asp Wat

Geo

1

0.0

0.0

0.0

0.0

0.0

0.0

TS1

2.2

-5.7

7.8

-0.7

-1.6

1.3

2

3.3

-7.5

6.5

2.5

-5.5

12.4

TS2

3.6

-7.4

6.6

0.8

-4.0

10.0

3

0.2

3.8

-14.3

7.6

-6.9

-1.7

TS3

0.0

-5.7

3.2

2.1

-4.7 -10.3

4

0.8 -10.9

11.4

2.5

-5.4

TS4

1.4 -12.7

13.6

2.7

-5.5

5

2.3 -11.4

13.5

0.9

-3.4

-4.9

Conclusions In this article, we have performed a QM/MM and big-QM study of the second part of heme degradation mechanism by heme oxygenase, i.e. the conversion of Fe(ІІІ) superoxo-oxophlorin to verdoheme. We suggest that the terminal oxygen atom attacks the carbon atom adjacent to keto group of oxophlorin, followed by a hemolytic cleavage of O2 and an attack of the same oxygen atom also on the other carbon adjacent to keto group of oxophlorin, which leads to the release of carbon monoxide and the formation of verdoheme (cf. Scheme 3 and Figure 2). This mechanism is similar to that suggested by QM-cluster calculations, but there are important differences, e.g. in the orientation of the His residue. Moreover, the formation of the second C–O bond in the verdoheme ring and the release of CO take place in separate steps in the QM/MM mechanism. In the QM/MM calculations, the reaction is exothermic by –35 kcal/mol and the

ACS Paragon Plus Environment

Page 19 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

rate-limiting barrier is 19 kcal/mol, whereas the corresponding data for the QM-cluster calculations were -15 and 12 kcal/mol, respectively. However, both calculations suggested that the first step of the mechanism is rate limiting. The QM/MM calculations are expected to give a better description of the energetic of mechanism through explicit inclusion of the steric and electrostatic effects of the protein environment. The water cluster in the distal pocket (which was not considered in the previous study) was found to play an important role, decreasing the reaction and activation energies by 2–7 kcal/mol. Furthermore, it is responsible for the orientation of the O2 group, by providing proper hydrogen bonds, directing the attack exclusively towards the α-meso carbon atom.

Acknowledgment This investigation has been supported by grants from the research council of Shahid Beheshti University and the Swedish research council (project 2014-5540). The computations were performed on computer resources provided by the Swedish National Infrastructure for Computing (SNIC) at Lunarc at Lund University and HPC2N at Umeå University.

ASSOCIATED CONTENT Supporting Information The supporting information is available free of charge via the Internet at http://pubs.ac.org. Email: [email protected] ; m-zahedi@sbu.,ac.ir Phone number: 98-21-29902889 ; 09127905750 Abbreviations DFT, density functional theory; HO, heme oxygenase; MD, molecular dynamics; QM, quantum mechanics; QM/MM, quantum mechanics/molecular mechanics.

Notes The authors declare no competing financial interest.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References 1. Unno, M.; Matsui, T.; Ikeda-Saito, M., Structure and catalytic mechanism of heme oxygenase. Nat. Prod. Rep. 2007, 24 (3), 553-570. 2. Uzel, C.; Conrad, M. E. In Absorption of heme iron, Semin. Hematol., 1998; pp 27-34. 3. Maines, M. D., Heme oxygenase: Function, multiplicity, regulatory mechanisms, and clinical applications. FASEB J. 1988, 2 (10), 2557-2568. 4. Maines, M. D., The heme oxygenase system: A regulator of second messenger gases. Annu. Rev. Pharmacol. 1997, 37 (1), 517-554. 5. Shibahara, S.; Müller, R.; Taguchi, H.; Yoshida, T., Cloning and expression of cdna for rat heme oxygenase. P. Natl. A. Sci. 1985, 82 (23), 7865-7869. 6. Maines, M. D.; Trakshel, G.; Kutty, R. K., Characterization of two constitutive forms of rat liver microsomal heme oxygenase. Only one molecular species of the enzyme is inducible. J. Biol. Chem. 1986, 261 (1), 411-419. 7. Mccoubrey, W. K.; Huang, T.; Maines, M. D., Isolation and characterization of a cdna from the rat brain that encodes hemoprotein heme oxygenase‐3. FEBS J. 1997, 247 (2), 725-732. 8. Ryter, S. W.; Alam, J.; Choi, A. M., Heme oxygenase-1/carbon monoxide: From basic science to therapeutic applications. Physiol. Rev. 2006, 86 (2), 583-650. 9. Ishikawa, K.; Matera, K. M.; Zhou, H.; Fujii, H.; Sato, M.; Yoshimura, T.; Ikeda-Saito, M.; Yoshida, T., Identification of histidine 45 as the axial heme iron ligand of heme oxygenase-2. J. Biol. Chem. 1998, 273 (8), 4317-4322. 10. McCoubrey, W. K.; Ewing, J. F.; Maines, M. D., Human heme oxygenase-2: Characterization and expression of a full-length cdna and evidence suggesting that the two ho-2 transcripts may differ by choice of polyadenylation signal. Arch. Biochem. Biophys. 1992, 295 (1), 13-20. 11. Liu, Y.; Moënne-Loccoz, P.; Loehr, T. M.; de Montellano, P. R. O., Heme oxygenase-1, intermediates in verdoheme formation and the requirement for reduction equivalents. J. Biol. Chem. 1997, 272 (11), 6909-6917. 12. Schuller, D. J.; Wilks, A.; Ortiz de Montellano, P. R.; Poulos, T. L., Crystal structure of human heme oxygenase-1. Nat Struct Mol Biol 1999, 6 (9), 860-867. 13. Ortiz de Montellano, P. R., Heme oxygenase mechanism: Evidence for an electrophilic, ferric peroxide species. Acc. Chem. Res. 1998, 31 (9), 543-549. 14. Ishikawa, K.; Takeuchi, N.; Takahashi, S.; Matera, K. M.; Sato, M.; Shibahara, S.; Rousseau, D. L.; Ikeda-Saito, M.; Yoshida, T., Heme oxygenase-2 properties of the heme complex of the purified tryptic fragment of recombinant human heme oxygenase-2. J. Biol. Chem. 1995, 270 (11), 6345-6350. 15. Migita, C. T.; Matera, K. M.; Ikeda-Saito, M.; Olson, J. S.; Fujii, H.; Yoshimura, T.; Zhou, H.; Yoshida, T., The oxygen and carbon monoxide reactions of heme oxygenase. J. Biol. Chem. 1998, 273 (2), 945-949. 16. Matsui, T.; Iwasaki, M.; Sugiyama, R.; Unno, M.; Ikeda-Saito, M., Dioxygen activation for the self-degradation of heme: Reaction mechanism and regulation of heme oxygenase. Inorg. Chem. 2010, 49 (8), 3602-3609. 17. Kumar, D.; de Visser, S. P.; Shaik, S., Theory favors a stepwise mechanism of porphyrin degradation by a ferric hydroperoxide model of the active species of heme oxygenase. J. Am. Chem. Soc. 2005, 127 (22), 8204-8213. 18. Balch, A. L., Coordination chemistry with meso-hydroxylated porphyrins (oxophlorins), intermediates in heme degradation. Coord. Chem. Rev. 2000, 200, 349-377. 19. Lai, W.; Chen, H.; Matsui, T.; Omori, K.; Unno, M.; Ikeda-Saito, M.; Shaik, S., Enzymatic ringopening mechanism of verdoheme by the heme oxygenase: A combined x-ray crystallography and qm/mm study. J. Am. Chem. Soc. 2010, 132 (37), 12960-12970. 20. Karplus, M.; McCammon, J. A., Molecular dynamics simulations of biomolecules. Nat Struct Mol Biol 2002, 9 (9), 646-652.

ACS Paragon Plus Environment

Page 20 of 23

Page 21 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

21. Jensen, F., Introduction to computational chemistry. John Wiley & Sons: 2013. 22. de Visser, S. P., Elucidating enzyme mechanism and intrinsic chemical properties of short-lived intermediates in the catalytic cycles of cysteine dioxygenase and taurine/α-ketoglutarate dioxygenase. Coord. Chem. Rev. 2009, 253 (5–6), 754-768. 23. Siegbahn, P. E. M.; Blomberg, M. R. A., Transition-metal systems in biochemistry studied by high-accuracy quantum chemical methods. Chem. Rev. 2000, 100 (2), 421-438. 24. Loew, G. H.; Harris, D. L., Role of the heme active site and protein environment in structure, spectra, and function of the cytochrome p450s. Chem. Rev. 2000, 100 (2), 407-420. 25. 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 (2), 949-1017. 26. Siegbahn, P. E.; Blomberg, M. R., Density functional theory of biologically relevant metal centers. Annu. Rev. Phys. Chem. 1999, 50 (1), 221-249. 27. Neese, F., A critical evaluation of dft, including time-dependent dft, applied to bioinorganic chemistry. JBIC, J. Biol. Inorg. Chem. 2006, 11 (6), 702-711. 28. Liao, M. S.; Watts, J. D.; Huang, M. J., Assessment of the performance of density‐functional methods for calculations on iron porphyrins and related compounds. J. Comput. Chem. 2006, 27 (13), 1577-1592. 29. Zhang, X.; Fujii, H.; Matera, K. M.; Migita, C. T.; Sun, D.; Sato, M.; Ikeda-Saito, M.; Yoshida, T., Stereoselectivity of each of the three steps of the heme oxygenase reaction: Hemin to mesohydroxyhemin, meso-hydroxyhemin to verdoheme, and verdoheme to biliverdin. Biochemistry 2003, 42 (24), 7418-7426. 30. Morishima, I.; Fujii, H.; Shiro, Y.; Sano, S., Studies on the iron (ii) meso-oxyporphyrin pi-neutral radical as a reaction intermediate in heme catabolism. Inorg. Chem. 1995, 34 (6), 1528-1535. 31. Matera, K. M.; Takahashi, S.; Fujii, H.; Zhou, H.; Ishikawa, K.; Yoshimura, T.; Rousseau, D. L.; Yoshida, T.; Ikeda-Saito, M., Oxygen and one reducing equivalent are both required for the conversion of-hydroxyhemin to verdoheme in heme oxygenase. J. Biol. Chem. 1996, 271 (12), 6618-6624. 32. Yoshinaga, T.; Sudo, Y.; Sano, S., Enzymic conversion of α-oxyprotohaem ix into biliverdin ixα by haem oxygenase. Biochem. J. 1990, 270 (3), 659-664. 33. Sano, S.; Sano, T.; Morishima, I.; Shiro, Y.; Maeda, Y., On the mechanism of the chemical and enzymic oxygenations of alpha-oxyprotohemin ix to fe. Biliverdin ix alpha. P. Natl. A. Sci. 1986, 83 (3), 531-535. 34. Gheidi, M.; Safari, N.; Zahedi, M., Density functional theory studies on the conversion of hydroxyheme to iron-verdoheme in the presence of dioxygen. Dalton Trans. 2017. 35. Gheidi, M.; Safari, N.; Zahedi, M., Chameleonic nature of hydroxyheme in heme oxygenase and its reactivity: A density functional theory study. Inorg. Chem. 2014, 53 (6), 2766-2775. 36. Gheidi, M.; Safari, N.; Zahedi, M., Effect of axial ligand on the electronic configuration, spin states, and reactivity of iron oxophlorin. Inorg. Chem. 2012, 51 (13), 7094-7102. 37. Ryde, U., The coordination of the catalytic zinc ion in alcohol dehydrogenase studied by combined quantum-chemical and molecular mechanics calculations. J. Comput.-Aided Mol. Des. 1996, 10 (2), 153-164. 38. Ryde, U.; Olsson, M. H., Structure, strain, and reorganization energy of blue copper models in the protein. Int. J. Quantum Chem. 2001, 81 (5), 335-347. 39. Ryde, U.; Olsen, L.; Nilsson, K., Quantum chemical geometry optimizations in proteins using crystallographic raw data. J. Comput. Chem. 2002, 23 (11), 1058-1070. 40. Treutler, O.; Ahlrichs, R., Efficient molecular numerical integration schemes. J. Chem. Phys. 1995, 102 (1), 346-354. 41. Case, D.; Darden, T.; Cheatham III, T.; Simmerling, C.; Wang, J.; Duke, R.; Luo, R.; Merz, K.; Pearlman, D.; Crowley, M., Amber 9 university of california. San Francisco 2006. 42. Claeyssens, F.; Harvey, J. N.; Manby, F. R.; Mata, R. A.; Mulholland, A. J.; Ranaghan, K. E.; Schütz, M.; Thiel, S.; Thiel, W.; Werner, H. J., High‐accuracy computation of reaction barriers in enzymes. Angew. Chem. 2006, 118 (41), 7010-7013.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

43. Lecerof, D.; Fodje, M.; Hansson, A.; Hansson, M.; Al-Karadaghi, S., Structural and mechanistic basis of porphyrin metallation by ferrochelatase. J. Mol. Biol. 2000, 297 (1), 221-232. 44. Svensson, M.; Humbel, S.; Froese, R. D.; Matsubara, T.; Sieber, S.; Morokuma, K., Oniom: A multilayered integrated mo+ mm method for geometry optimizations and single point energy predictions. A test for diels-alder reactions and pt (p (t-bu) 3) 2+ h2 oxidative addition. J. Phys. Chem. 1996, 100 (50), 19357-19363. 45. Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E., Climbing the density functional ladder: Nonempirical meta–generalized gradient approximation designed for molecules and solids. Phys. Rev. Lett. 2003, 91 (14), 146401. 46. Schäfer, A.; Horn, H.; Ahlrichs, R., Fully optimized contracted gaussian basis sets for atoms li to kr. J. Chem. Phys. 1992, 97 (4), 2571-2577. 47. 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 (18), 3297-3305. 48. Grimme, S.; Ehrlich, S.; Goerigk, L., Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32 (7), 1456-1465. 49. Eichkorn, K.; Treutler, O.; Öhm, H.; Häser, M.; Ahlrichs, R., Auxiliary basis sets to approximate coulomb potentials. Chem. Phys. Lett. 1995, 240 (4), 283-290. 50. Becke, A. D., A new mixing of hartree–fock and local density‐functional theories. J. Chem. Phys. 1993, 98 (2), 1372-1377. 51. Hertwig, R. H.; Koch, W., On the parameterization of the local correlation functional. What is becke-3-lyp? Chem. Phys. Lett. 1997, 268 (5), 345-351. 52. Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R., Auxiliary basis sets for main row atoms and transition metals and their use to approximate coulomb potentials. Theor. Chem. Acc. 1997, 97 (1), 119-124. 53. Blomberg, M. R.; Borowski, T.; Himo, F.; Liao, R.-Z.; Siegbahn, P. E., Quantum chemical studies of mechanisms for metalloenzymes. Chem. Rev. 2014, 114 (7), 3601-3658. 54. Ryde, U., Accurate metal-site structures in proteins obtained by combining experimental data and quantum chemistry. Dalton Trans. 2007, (6), 607-625. 55. Delcey, M. G.; Pierloot, K.; Phung, Q. M.; Vancoillie, S.; Lindh, R.; Ryde, U., Accurate calculations of geometries and singlet–triplet energy differences for active-site models of [nife] hydrogenase. Phys. Chem. Chem. Phys. 2014, 16 (17), 7927-7938. 56. Dong, G.; Ryde, U., O2 activation in salicylate 1, 2-dioxygenase: A qm/mm study reveals the role of his162. Inorg. Chem. 2016, 55 (22), 11727-11735. 57. Li, J.-L.; Mata, R. A.; Ryde, U., Large density-functional and basis-set effects for the dmso reductase catalyzed oxo-transfer reaction. J. Chem. Theory Comput. 2013, 9 (3), 1799-1807. 58. Case, D.; VB JTB, B. R.; Cai, Q.; Cerutti, D.; Cheatham III, T.; Darden, T.; Duke, R.; Gohlke, H.; Goetz, A.; Gusarov, S., The ff14sb force field. Amber 2014, 14, 29-31. 59. Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A., Development and testing of a general amber force field. J. Comput. Chem. 2004, 25 (9), 1157-1174. 60. Besler, B. H.; Merz, K. M.; Kollman, P. A., Atomic charges derived from semiempirical methods. J. Comput. Chem. 1990, 11 (4), 431-439. 61. Schuller, D. J., Wilks, A., Ortiz de Montellano, P.R., Poulos, T.L., Comparison of the heme-free and -bound crystal structures of human heme oxygenase-1. http://www.rcsb.org/pdb/explore/explore.do?structureId=1N45. 62. Li, H.; Robertson, A. D.; Jensen, J. H., Very fast empirical prediction and rationalization of protein pka values. Proteins: Struct., Funct., Bioinf. 2005, 61 (4), 704-721. 63. Bartlett, G. J.; Porter, C. T.; Borkakoti, N.; Thornton, J. M., Analysis of catalytic residues in enzyme active sites. J. Mol. Biol. 2002, 324 (1), 105-121. 64. Hu, .; S derh elm, P. r.; Ryde, U., On the convergence of qm/mm energies. J. Chem. Theory Comput. 2011, 7 (3), 761-777.

ACS Paragon Plus Environment

Page 22 of 23

Page 23 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

65. Sumner, S.; S derh elm, P. r.; Ryde, U., Effect of geometry optimizations on qm-cluster and qm/mm studies of reaction energies in proteins. J. Chem. Theory Comput. 2013, 9 (9), 4205-4214. 66. Hu, L.; Soderhjelm, P.; Ryde, U., Accurate reaction energies in proteins obtained by combining qm/mm and large qm calculations. J. Chem. Theory Comput. 2012, 9 (1), 640-649. 67. Becke, A. D., Density‐functional thermochemistry. Iii. The role of exact exchange. J. Chem. Phys. 1993, 98 (7), 5648-5652. 68. Sierka, M.; Hogekamp, A.; Ahlrichs, R., Fast evaluation of the coulomb potential for electron densities using multipole accelerated resolution of identity approximation. J. Chem. Phys. 2003, 118 (20), 9136-9148.

ACS Paragon Plus Environment