Iron Hydroperoxide Intermediate in Superoxide ... - ACS Publications

May 16, 2017 - ... CEA/DRF/BIG/CBM/MCT, CNRS UMR 5249, Université Grenoble Alpes, ... Superoxide reductase is a mononuclear iron enzyme involved in ...
0 downloads 0 Views 5MB Size
Subscriber access provided by UB + Fachbibliothek Chemie | (FU-Bibliothekssystem)

Article

Iron hydroperoxide intermediate in Superoxide Reductase: Protonation or Dissociation First? MM Dynamics and QM/MM Metadynamics study. Rolf David, Helene Jamet, Vincent Niviere, Yohann Moreau, and Anne Milet J. Chem. Theory Comput., Just Accepted Manuscript • Publication Date (Web): 16 May 2017 Downloaded from http://pubs.acs.org on May 20, 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.

Journal of Chemical Theory and Computation 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 73

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

Journal of Chemical Theory and Computation

Iron hydroperoxide intermediate in Superoxide Reductase: Protonation or Dissociation First? MM Dynamics and QM/MM Metadynamics study. Rolf David, †,‡,§ Hélène Jamet, †,‡ Vincent Nivière, ǁ Yohann Moreau, § Anne Milet *,†,‡ † Univ. Grenoble Alpes, DCM, F-38000 Grenoble, France ‡ CNRS, DCM, F-38000, Grenoble, France § Laboratoire de Chimie et Biologie des Métaux, CEA/DRF/BIG/CBM/MCT, CNRS UMR 5249, Université Grenoble Alpes, Grenoble, France

ǁ Laboratoire de Chimie et Biologie des Métaux, CEA/DRF/BIG/CBM/BioCat, CNRS UMR 5249, Université Grenoble Alpes, Grenoble, France

Abstract Superoxide reductase is a mononuclear iron enzyme involved in superoxide radical detoxification in some bacteria. Its catalytic mechanism is associated with the remarkable formation of a ferric hydroperoxide Fe3+-OOH intermediate, which is specifically protonated on its proximal oxygen to generate the reaction product H2O2. Here we present a computational study of the protonation mechanism of the Fe3+-OOH intermediate, at different levels of theory.

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation

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 2 of 73

This was performed on the whole system (solvated protein) using well-tempered metadynamics at the QM/MM (B3LYP/AmberFF99SB) level. Enabled by the development of a new set of forcefield parameters for the active site, a conformational MM study of the Fe3+-OOH species gave insights into its solvation pattern, in addition of generating the two starting conformations for the ab initio metadynamics setup. Two different protonation mechanisms for the Fe3+-OOH intermediate have been found depending on the starting structure. Whereas a possible mechanism involves at first the protonation of the hydroperoxide ligand and then dissociation of H2O2, the most probable one starts with an unexpected dissociation of the HOO- ligand from the iron, followed by its protonation. This favored reactivity was specifically linked to the influence of both the nearby conserved lysine 48 residue and the micro-solvatation on the charge distribution of the oxygens of the HOO- ligand. These data highlight the crucial role of the whole environment, solvent and protein, to describe accurately this second protonation step in superoxide reductase. This is clearly not possible with smaller models unable to reproduce correctly the mechanistically determinant charge distribution.

Introduction Oxidative stress is a major problem encountered by living organisms. It originates from an abnormal formation of the so-called reactive oxygen species (ROS) in the cells, such as the superoxide radical (O2●-), the hydrogen peroxide (H2O2) or the hydroxyl radical (HO●)1. Although ROS also play crucial functions in cell physiology, it is now well documented that their overproduction and/or mismanagement are capable of inflicting serious biological damages, which could be implicated in the genesis of human pathologies such as cancers, cardiovascular, inflammatory and degenerative diseases1,2. To counteract these toxic effects of ROS, living organisms have developed whole enzymatic machineries aimed at catalyzing very efficiently

ACS Paragon Plus Environment

2

Page 3 of 73

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

Journal of Chemical Theory and Computation

their elimination3. Concerning the superoxide radical O2●-, for years, superoxide dismutase enzymes (SOD) were reported to be the only enzymatic system known to detoxify it, catalyzing O2●- dismutation into hydrogen peroxide and molecular oxygen3. Nevertheless, more recently, it was discovered that some bacteria use an alternative system to eliminate O2●-, involving superoxide reductase (SOR). This iron-containing enzyme catalyzes the one electron reduction of superoxide into hydrogen peroxide without production of molecular oxygen, and thus involves a different reaction from SOD4–9. The electron involved in this reaction is provided by cellular reductases, acting as physiological partners for SOR activity: O2●- + 1 e - + 2 H+  H2O2 (SOR) Although SOD and SOR are unrelated enzymes9, they both scavenge superoxide with an extremely fast reaction rate close to the diffusion limit, exhibiting in vivo a similar efficiency in superoxide detoxification10. Due to the critical importance of ROS detoxification processes in the cells, gaining information about these systems is a great challenge. Whereas the reaction mechanisms of superoxide scavenging by SODs have been intensively documented, that of the more recently characterized SOR still require substantial clarifications9. The active site of SOR consists in an atypical non-heme mononuclear iron center coordinated in a square pyramidal geometry by four equatorial nitrogens from histidines and a sulfur from a cysteine in the axial position ([FeN4S1])6–9. The iron active site is very solvent accessible since it is positioned at the surface of the enzyme (i.e. not buried as in SOD, e.g.). Although the active site is identical in all the SORs characterized so far, some enzymes like those from Desulfoarculus baarsii11 and Desulfovibrio vulgaris6, possess a second mononuclear iron center localized in an additional N-terminal domain6,11. This center is of desulforedoxin type (DX), with an iron chelated by four sulfurs from cysteines in a rubredoxin-distorted manner ([FeS4])12. It is

ACS Paragon Plus Environment

3

Journal of Chemical Theory and Computation

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 4 of 73

located at 22 Å from the active site, and up to now when present in SOR, its function in catalysis has not been yet definitively established, although it has been proposed to act as an electron relay toward the active site13–15. The catalytic cycle of SOR has been extensively studied both experimentally16–20 and theoretically20–24 and it is usually described in four main steps, as reported for the well-characterized enzyme from D. baarsii20 (Scheme 1). Along the catalytic cycle, the iron center remains in a high spin electronic configuration20.

Scheme 1. Native catalytic cycle of Superoxide Reductase

Starting from the reduced resting state (Fe2+, top structure in Scheme 1), the superoxide substrate binds to the sixth position of iron to form a first intermediate, which has been proposed to be a ferrous iron-superoxide (IUPAC: superoxido) intermediate Fe2+-OO●-20,23. This first

ACS Paragon Plus Environment

4

Page 5 of 73

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

Journal of Chemical Theory and Computation

intermediate is then protonated on its distal oxygen (Od), allowing a one electron transfer from the ferrous iron to the superoxo adduct to form a ferric iron-hydroperoxide (IUPAC: hydrogenperoxido) Fe3+-OOH species (second intermediate, at the bottom of Scheme 1). In this first protonation step, there are no clear experimental evidences for the involvement of second coordination sphere amino-acid residues, the proton could be directly provided by H3O+ from the bulk solvent16–18,20. The Fe3+-OOH intermediate, where the hydroperoxide binds in an end-on fashion to the iron center25, has been experimentally trapped and well characterized in some mutants of the SOR from D. baarsii by the mean of high resolution X-Ray structures25, resonance Raman26–28 and Mössbauer spectroscopies29. A second protonation step, specifically on the proximal oxygen (Op) of the hydroperoxide moiety, was then proposed to lead to the formation and subsequent release of the H2O2 product20. However, this second protonation step could not be directly observed and studied experimentally, and it was proposed that it could correspond to a very fast, non-rate limiting step20. Nevertheless, X-ray structures of the Fe3+-OOH intermediate suggested that Lys48 (numbered accordingly to D. baarsii SOR sequence), a quasi conserved residue located close to the iron site*, could play a critical role in this second protonation step25. These structures showed that Lys48 is involved in a H-bond network with the Op of the Fe3+-OOH, via one molecule of water. Such a conformation may favor a specific protonation of the proximal oxygen of the intermediate to form the final product H2O2. In line with this hypothesis, mutation of Lys48 into isoleucine (K48I mutant) in the SOR from D. baarsii was shown to modify its reactivity with superoxide20. In the K48I mutant, the first reaction intermediate Fe2+-OO●- is still formed and its decay still remains associated with a protonation step on its distal oxygen.

*

In a very few examples, such as in the SORs from Ignicoccus hospitalis, this lysine residue is not conserved106.

ACS Paragon Plus Environment

5

Journal of Chemical Theory and Computation

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 6 of 73

However, the H2O2 product is not formed any more. It was further shown that in the K48I SOR mutant, the Fe3+-OOH intermediate evolves into a high valent iron-oxo species (Fe=O), through a O-O bond cleavage of the hydroperoxide adduct30. Formation of such iron-oxo species in the K48I mutant was associated with a lack of a specific protonation on the proximal oxygen of the Fe3+-OOH species, most-likely favoring a double protonation of its distal oxygen and a subsequent O-O bond cleavage30. In cytochrome P450 oxygenase, such a double protonation process of the distal oxygen of the Fe3+-OOH intermediate (compound 0 in P450) which leads to the formation of a high valent iron-oxo species (compound I in P450) was described as one of the key steps in its reaction mechanism, allowing the catalysis of oxidation of very stable chemical bonds of substrates31. Thus, the mechanisms by which these enzymes control the protonation of the Fe3+-OOH intermediate appears determining for their activity, detoxification for SOR, oxygenation for cytochrome P450. Finally, after its release from the iron, H2O2 is replaced by an OH-/ H2O ligand, itself further replaced by the binding of the side chain of the loop-well conserved glutamate residue (Glu47, in D. baarsii)18,20. The catalytic cycle can restart after the 1 e- reduction by cellular reductases of the ferric resting state into the ferrous one10 (Scheme 1). Lys48 is not the only neighboring residue that helps tuning the reactivity of SOR. The X-ray structures of SOR have shown that Ile118 makes H-bond interaction through its NH main chain with the sulfur atom of the cysteinate ligand11,25. Site-directed mutagenesis experiments on this residue associated with DFT calculations24 showed that this H-bond network modulates the ”push effect” of the cysteinate ligand32, namely the electro-donor properties of the thiolate to the iron. In particular, the modulation of this “push effect” was shown to control the pKa of the first reaction intermediate and its rate of protonation24.

ACS Paragon Plus Environment

6

Page 7 of 73

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

Journal of Chemical Theory and Computation

Thus, the whole SOR machinery is optimized to favor the efficient reduction of superoxide into hydrogen peroxide, and the involvement of the two protons to form H2O2 represents a central point in its catalytic mechanism. However, as mentioned above, unlike the properties of the active site and its surrounding, these protonation steps required to convert O2●- into H2O2, and particularly the second one, were much less well addressed experimentally20. Theoretical calculations, in order to overcome the lack of experimental information, are particularly well suited to gain into detailed information about these critical protonation processes. Previous work based on theoretical calculations23 presented the investigation of the first protonation step, i.e. protonation of the Od of the Fe2+-OO●- intermediate to generate the Fe3+-OOH species. It was proposed that Lys48 might act as a proton relay during this first protonation step, either with a direct interaction with the substrate or via a H-bond network involving one or several water molecule(s). To assess the first proton transfer, the authors performed QM/MM molecular dynamics simulations as well as QM calculations. In term of methodology, this work showed the importance of a proper description of the spin state of the iron and assessed the validity of several levels of calculations including the use of the B3LYP33– 35

functional. In the present work, we present a series of theoretical calculations focused on the mechanism

of the second protonation step of SOR, i.e. protonation of the key Fe3+-OOH intermediate to generate the H2O2 reaction product. To get a clear view of the whole mechanism as well as of the impact of the environment around the active site and the H-bonding network with solvent molecules, we opted for a QM/MM dynamic study of the solvated protein along with an exploratory MM study of the H-bonding network dynamic on the Fe3+-OOH intermediate.

ACS Paragon Plus Environment

7

Journal of Chemical Theory and Computation

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 8 of 73

To accelerate the dynamic and describe the reactivity, we opted for metadynamics36,37 simulations, which have shown to play a major role in understanding the mechanism and the energetics of reactions38–42 especially for complex systems such as biochemical compounds43–46. Indeed, the bias introduced to accelerate the rare events usually allows enough flexibility to take into account several possible mechanisms. These calculations allowed the description of different protonation pathways involving directly or indirectly the lysine 48 residue. We showed that Lys48 has a major impact on the charges of the proximal oxygen of the Fe3+-OOH intermediate. Coupled with the micro-solvation of water around the active site, our data showed that Lys48 triggers a totally unexpected protonation mechanism: the dissociation of the hydroperoxo moiety from the iron atom occurring at first, followed in a second step by its protonation and thus forming the hydrogen peroxide reaction product.

Computational details and methods Params setup To model the two iron centers of the SOR from Desulfoarculus baarsii (active site and DX centers) we used the MCPB47 program in the AmberTools1248 package. We made the choice to describe the metal-ligand interaction with a bonded model. Bond-stretching & angle-bending force constants were determined from submatrices of the Cartesian Hessian following a methodology described by Seminario49. Hoops et al50 showed that dihedral parameters were not needed to describe metallic center. Dihedral angle force constants involving iron were thus set to zero. Charges on atoms were calculated using the RESP51 approach using the electrostatic potential sampled in the Merz-Kollman52 scheme. Following the MCPB protocol, Gaussian0353

ACS Paragon Plus Environment

8

Page 9 of 73

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

Journal of Chemical Theory and Computation

was used to determine the forces subsystems (for other versions of Gaussian frequencies calculation output was not supported by MCPB at the time) whereas Gaussian0954 was used for the charges subsystems. Subsystems used to define the force constants (sidechain model) includes the high spin irons and the first sphere of amino acids and are shown in Figure 1 (See SI for Cartesian coordinates).

Figure 1. Side chain models for the D. baarsii SOR. Active site in the Fe3+-OOH state (bottom) and DX site (top). Highlighted carbon atoms are kept fixed during geometries optimization. Colors for iron, sulfur, oxygen, nitrogen, carbon and hydrogen atoms are respectively ochre, yellow, red, blue, gray and white.

ACS Paragon Plus Environment

9

Journal of Chemical Theory and Computation

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 73

Structures of active sites were extracted from the chain B of the structure with PDB entry 2JI325, and an optimization of the geometry has been performed at the B3LYP33–35/6-31+G(d)55– 62

level of theory with the alpha carbons frozen at their crystallographic position. The Hessian

matrices were determined via harmonic frequency calculations, which also allowed confirming the structures to be being actual minima on the potential energy surface. To determine atomic charges, larger model were used (see Figure 2 and SI for the Cartesian coordinates) with the addition of the N-methyl amide capping group or the capping acetyl group in order to ensure the consistency of the AMBER charges fitting procedure used in the FF9X63 force fields. An optimization of geometry was also performed at the B3LYP/6-31G(d) level with the C-N-O-C of the backbone frozen followed by frequency calculation. Mulliken charges were calculated at the B3LYP/6-311G(d,p)58,60,64–67 level and were fitted in two steps. A first charge fitting without any restriction was performed. The second fitting with the inclusion of the charges on the backbones atoms (including hydrogens) of the amino-acids ligands constrained to their Amber FF9463 values. Also, in this step, the charge on the iron was fixed to its unconstrained fitting values (from the first step) to be consistent with the oxidation number of the iron. The van der Waals non-bonded parameters were taken from reference 68 for Fe2+ and from reference 69 for Fe3+. Following this procedure, three sets of parameters were obtained for the DX center (the iron center is in the formal oxidation state Fe3+ and in a high spin configuration) and the two states for the SOR active site (SOR-Fe2+-OO●- & SOR-Fe3+-OOH). (Complete sets of parameters available in the Supporting Information, section I)

ACS Paragon Plus Environment

10

Page 11 of 73

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

Journal of Chemical Theory and Computation

Figure 2. Large model for the DX (top) and SOR-Fe3+-OOH (bottom) sites. Highlighted atoms are kept fixed during geometries optimization. Colors for iron, sulfur, oxygen, nitrogen, carbon and hydrogen atoms are respectively ochre, yellow, red, blue, gray and white.

MM Setup To check the stability of these three sets of parameters (SOR-Fe2+-OO●-, SOR-Fe3+-OOH & DX), and to get insight on the structures of the SOR centers, classical molecular dynamics simulations in explicit solvatation (water) were run for the dimeric system including chains A and B, extracted from the 2JI3 X-Ray structure25, but replaced back to the wild type. (Ala114 –>

ACS Paragon Plus Environment

11

Journal of Chemical Theory and Computation

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 73

Glu114). The two systems (SOR-Fe2+-OO●- & SOR-Fe3+-OOH) were neutralized by adding the needed number of Na+ ions and then solvated with a TIP3P70 cubic box of water with a solvation layer of the enzyme of 8.0 Å. The molecular dynamics trajectories were produced using the Amber1248 package with the Amber force field FF99SB71. Periodic boundary conditions were applied and a smooth particle mesh Ewald72,73 method was used to calculate non-bonded interactions with a cutoff of 10 Å. The whole system was minimized in two steps: a first one with the protein fixed via a harmonic constraint with a global force constant of 500 kcal mol-1 Å-2 in order to minimize the solvent and suppress close contacts. The second minimization was performed with no constraint. For all classical molecular dynamics simulations a timestep of 0.5 fs was used, for consistency with further QM/MM calculations. After both minimizations, the system was heated 25 ps to reach an equilibrium temperature of 300 K, using a Langevin thermostat74,75 and a friction constant of 5 fs-1 in the NVT ensemble. This step was followed by 975 ps of equilibration in the NPT ensemble, to ensure a correct density of 1. Production dynamics in the NVT ensemble during 100 ns were then performed.

QM/MM Setup To explore the reaction mechanisms along the protonation step of the SOR-Fe3+-OOH intermediates, a series of QM/MM76 dynamic simulations were run with the CP2K77,78 program. In all the QM/MM simulations, the QM part (only one SOR- Fe3+-OOH site from only one chain, A or B) was constituted of the ferric iron, the hydroperoxide substrate, the side chains of residues His49, His69, His75, His119 and Cys116 (i.e. the first sphere of coordination of the metal) and the one of Lys48. All the water molecules at a distance of 5 Å of the hydroperoxide moiety as well as the ones within 3 Å of the amino group of Lys48 were also considered. The QM/MM

ACS Paragon Plus Environment

12

Page 13 of 73

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

Journal of Chemical Theory and Computation

boundaries, between the Cα and Cβ of the five amino acids binding Fe and the Lys48 were treated using the link atom scheme79. The GEEP80,81 scheme was used to calculate the QM/MM interactions. The QM cell has the size of the QM subsystems plus a distance of 5 Å on every direction. The all-electron GAPW82,83 scheme was employed at the B3LYP level using a truncated coulomb operator84 of 6 Å. B3LYP, a hybrid density functional, was preferred over pure functionals due to its ability to provide the right spin state for the system. Dispersion was added by the empirical dispersion correction D385,86. The 6-31G(d,p) basis set was used for all atoms except for the oxygen atoms of the substrate, the four nitrogens and the sulfur binding the iron and the nitrogen of the lysine side, for which 6-31+G(d,p) basis set was used. A cutoff of 400 Ry was retained to project the electronic structure. The MM part has the same parameters as the full MM setup. Periodic boundary conditions were applied to the MM and the QM system87. Points of interest from the MM runs were extracted; geometries and velocities were kept to start Born-Oppenheimer MD simulations for 5 ps in the NVT ensemble with a Nose-Hoover thermostat88–90 with a 100 fs-1 constant. The time step for QM/MM simulations was 0.5 fs.

Metadynamics QM/MM Setup Well-tempered metadynamics36,37 was used to explore and later reconstruct the free-energy surface along the proton transfer process. The level of theory is the same as the one used for the QM/MM dynamics calculations, B3LYP/6-31G(d,p) and B3LYP/6-31+G(d,p) for the QM part (vide supra). Different sets of two collective variables (CV, hereafter) were used to describe the reaction depending on the initial conformation and they were described in the results part. The Gaussian bias potential parameters used in this study were the following: the height of the Gaussian was set to 1.0 kcal mol-1 and the width was set to 0.05 Å. This corresponds to an

ACS Paragon Plus Environment

13

Journal of Chemical Theory and Computation

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 73

oscillation of the CV in unbiased dynamics, with a deposition step of 10 fs (every 20 MD steps). Quadratic function walls are applied at various distances with a K constant of 500 kcal mol-1 to restrict the exploration of the subspace. The ∆T parameter of the well-tempered metadynamic (for restricting the exploration to an energy range of T + ∆T) was set to 6000 K (vide infra).

Results and discussions In order to obtain a global and accurate view of the mechanism of the second protonation step in the SOR reaction cycle (Scheme 1), as well as the significant factors governing this reactivity, we have first gained insights into the H-bond network around the active site. This has been carried out through a classical molecular dynamics study enhanced by the set-up of new parameters for the iron atoms, whose validity was assessed through comparisons with experimental and theoretical works. Then, starting from selected configuration issued from the classical simulations, the protonation of the SOR Fe3+-OOH intermediate was simulated at the QM/MM level, using metadynamics at the B3LYP/MM level.

H-Bonding network around the active site Validity of the models The SOR from the Desulfoarculus baarsii contains two non-heme mononuclear irons, the rubredoxin-like desulforedoxin (DX) center [FeS4] and the active site center [FeN4S1]. Although the DX iron center is located quite far from the SOR active site (22 Å), with a role in catalysis not clear yet, this center was nevertheless taken into account for the construction of our SOR models. Some parameterizations already exist in the literature for the DX center in the Desulfovibrio gigas protein91, but to ensure consistency between the two SOR metallic centers,

ACS Paragon Plus Environment

14

Page 15 of 73

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

Journal of Chemical Theory and Computation

we chose to perform a complete parameterization for both centers as no parameterization with similar protocol was available. Upon reaction with superoxide, the SOR active site becomes hexacoordinated, with the binding of the superoxide substrate on the sixth coordination position. This leads to the formation of at least two SOR reaction intermediates, the ferrous iron superoxide species (Fe2+-OO●-) and the ferric iron hydroperoxide species (Fe3+-OOH) species20 (Scheme 1). The sets of parameters which have been developed here for the iron centers, with the active site in different conformations and bound intermediates, were used in a series of classical molecular dynamics simulations to assess their validity by comparing structural results to experimental data (X-Ray structures) and to previous calculations. For the DX site, the force constants provided by our models are in good agreement with the parameterization work of Carvalho et al91. For example, the bond equilibrium distance for the Fe-S distance is found to be 2.3 Å in both our work and the Carvalho’s one, with force constant of ≈ 80 kcal mol-1 Å-2 for each of the four cysteine ligand. Considering the SOR active site, two sets of parameters had to be derived, depending on the nature of the reaction intermediates. A first set was thus determined for the Fe2+-OO●- species and a second one for the Fe3+-OOH species. Two noticeable differences between the two sets of parameters can be pointed out. First, the force constant for the Fe-Op bond is greater in the Fe3+-OOH structure than in the Fe2+-OO●- one, with, respectively, 76.4 kcal mol-1 Å-2 and 25.5 kcal mol-1 Å-2. Consistently, the Fe-Op equilibrium distances are 2.0 and 2.2 Å, in the Fe3+-OOH and Fe2+-OO●- species, respectively. Second, the Op-Od bond distances differ between the two states. In Fe2+-OO●-, the equilibrium distance is 1.33 Å with a force constant of 483.1 kcal mol-1 Å-2, while the corresponding values in Fe3+-OOH are 1.46 Å and 265.5 kcal

ACS Paragon Plus Environment

15

Journal of Chemical Theory and Computation

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 73

mol-1 Å-2, respectively. These differences are consistent with the evolution from a Fe2+-OO●- species towards a Fe3+-OOH one (Scheme 1), as already pointed out by previous QM calculations24. Using these sets of parameters, classical MD simulations were carried out, considering the active site of SOR in either Fe2+-OO●- or Fe3+-OOH states. Because the structure of the D. baarsii SOR is homodimeric†,25 (i.e. chains A and B in PDB 2JI3), the sampling for analysis was doubled for each SOR active site (i.e. presence of two chains per simulation) and fourfold for the DX site (two chains per simulation and two simulations : one for each SOR active site state). This means there are four DX samples, two Fe2+-OO●- samples and two Fe3+-OOH samples. Root Mean Squared Deviations (RMSD) of the relevant distances were computed considering the 100 ns of the MD production run and are presented in the Figure 3. The RMSD for all the atoms of the protein were found to be stable for the simulation time within 1.8 Å (Fig. 3 d, in red), and within 0.5 Å for each active site as an average of both chains (i.e. the metallic center and the amino acids directly bound to it, Fig. 3 a, b & c, in black).



The dynamics simulation of only one chain led to the unfolding of the protein meaning that two chains are essentials to keep the whole structure of the D. baarsii SOR during simulation.

ACS Paragon Plus Environment

16

Page 17 of 73

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

Journal of Chemical Theory and Computation

Figure 3. (a) RMSD of the SOR-Fe2+-OO●- active site. (b) RMSD of the SOR-Fe3+-OOH active site. (c) RMSD of the DX site. (d) RMSD of the entire protein. SOR DX site. To assess the validity of the DX parameters, we compared our simulations with previous experimental results and calculations. Concerning the experimental data, we used two sets of X-Ray structures: the DX center present in the structure of the SOR from D. baarsii with the Fe3+-OOH intermediate (PDB entry 2JI3) and the one in the desulforedoxin from D. gigas (PDB entry 1DXG) obtained by Archer et al12. We also compared our results to the previous theoretical work of Carvalho et al91, which developed parameterization for the DX iron center in 1DXG. Note that the DX center is the only metallic center present in 1DXG. Our set of parameters allows a good description of the DX center within the SOR from D. baarsii and is in good

ACS Paragon Plus Environment

17

Journal of Chemical Theory and Computation

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 73

agreement with the one in the desulforedoxin from D. gigas. Our data also reproduce well the H-bonding network around the DX site (More details can be found in the SI, section II.1).

Fe2+-OO●- site. Due to the lack of direct crystallographic data, our simulations of the SOR active site in its ferrous-superoxide state could be solely compared to previous QM models described in the literature23,24. Our calculations were in very good agreement with these two models, confirming the validity of the protocol for the parameterization of the Fe2+-OO●- state. The details of the models used, geometry and active site H-bonding can be found in the SI.

Fe3+-OOH site. Our values calculated for the iron-hydroperoxide site were compared to the crystallographic data from 2JI325, to the small model of Sit et al23 (made of 45 atoms, imidazoles & methyl thiolate, gas phase) and to the large one by Tremey et al24 (made of 88 atoms, Cα histidines, cysteine & backbone atoms between Cys116 and His119 (C model in Ref

24

) – implicit water).

Table 1 gathers characteristics distances for the active site with the hydroperoxide moiety bound to the iron, determined from our calculations, the 2JI3 crystallographic structure and both QM models. For the Fe-S bond, we obtained a value of 2.42 ± 0.1 Å. The averaged distance in the crystal 2JI3 is 2.53 ± 0.1 Å25. The large QM model reported a value of 2.53 Å24. No value was reported for the model of Sit et al23. The Fe-Op distance is 2.04 ± 0.1 Å in our simulations. In the crystal structure, the Fe-Op distance is 2.00 ± 0.1 Å whereas the small QM model reported a value of 2.00 Å and the large one a value of 1.99 Å. For the Fe-Op-Od angle, we found a value of 126 ± 5°, which is close to the X-Ray value of 125 ± 3°25. Values of 123.9° for the small

ACS Paragon Plus Environment

18

Page 19 of 73

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

Journal of Chemical Theory and Computation

model23 and 119.1° for the large model24 were reported. Altogether, these data show that, our parameters preserve very well the structure of the SOR active site, which is found to be equivalent for chain A or chain B. The value of the Cβ-S-Op-Od dihedral angle will be discussed later in this work. Table 1. Relevant distances and angles (bending and dihedral) values for the Fe3+-OOH species from our simulations (averaged), QM models and the X-Ray structure (averaged over the three (hydro)peroxide-containing chains). Distance (Å)

Fe-Op

S-Fe

Op-Od

This Work (averages)

2.04 ± 0.1

2.42 ± 0.1

1.47 ± 0.0

Tremey et al Work24

1.99

2.52

1.44

Sit et al Work23

2.00

X-Ray 2JI325 (averages)

2.00 ± 0.1

2.53 ± 0.1

1.47 ± 0.0

Angle (°)

Fe-Op-Od

Op-Od-Hd

Cβ-S-Op-Od

Fe-Op-Od-Hd

This Work (averages)

126 ± 5

107 ± 5

-50 ± 21

-126 ± 56

Tremey et al Work24

119.1

100.7

105.1

-171.1

Sit et al Work23

123.9

X-Ray 2JI325 (averages)

125 ± 3

1.47

128 ± 15

The Hydrogen-bond network‡ around the cysteine sulfur ligand was also checked for our model. Similarly to the Fe2+-OO●- intermediate, we found Ile118 and His119 residues to be engaged with the sulfur ligand in a hydrogen bond network via their amide backbones. In chain A, the involvements of Ile118 and His119 in H-bonding with the sulfur ligand are 26% and 65%, respectively, while for the B chain site, the corresponding values are 30% and 50%. Relevant



Parameters pf H-bond analysis are detailled in the SI, section I.1.

ACS Paragon Plus Environment

19

Journal of Chemical Theory and Computation

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 20 of 73

bonds and angles are reported in Table 2 and our simulations showed no difference between chains A and B of SOR. Table 2. Hydrogen bond distances between the sulfur ligand of the SOR-Fe3+-OOH species and the main chain amide groups of Ile118 and His119, calculated in this work and in comparison with the values reported in the X-Ray 2JI3 structure25. Acceptor

Chain A

Chain B

2JI325 Averages over the three (hydro)peroxidecontaining chains

Donor

Distance H-A (Å)

Distance N-A (Å)

Angle D-H-A (°)

S Cys116

H Ile118

N Ile118

2.34 ± 0.2

3.27 ± 0.1

153 ± 10

S Cys116

H His119

N His119

2.20 ± 0.1

3.17 ± 0.1

162 ± 11

S Cys245

H Ile247

N Ile247

2.46 ± 0.2

3.39 ± 0.2

154 ± 11

S Cys245

H His248

N His248

2.30 ± 0.2

3.24 ± 0.2

157 ± 15

S Cys116

H Ile118

N Ile118

3.41 ± 0.1

S Cys116

H His119

N His119

3.20 ± 0.1

Description of the H-Bonding network around the Fe3+–OOH intermediate. One point of interest was the role of the water molecules in interaction with Lys48, whose importance has been foreseen due to its closeness to the OOH- moiety in the Fe3+-OOH species. Along our MD simulations, Lys48 was found to interact with the superoxide or hydroperoxide intermediates either directly or indirectly via one or more molecules of water with the superoxide or hydroperoxide moiety. When considering the interaction and H-bond network between lysine and the –-O-O● or the OOH- moiety, three different patterns were found to be relevant. The first one involved a direct lysine-to-intermediates interaction and was noted ‘direct’. In the second one, one molecule of water was found between the lysine and the intermediates and was noted ‘one water bridge’. In

ACS Paragon Plus Environment

20

Page 21 of 73

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

Journal of Chemical Theory and Computation

the third one, two interstitial molecules of water were inserted between the lysine and the intermediates and was noted ‘two water bridge’. Figure 4 depicts examples of these configurations for the hydroperoxide§.

Figure 4. Three different H-bond networks between lysine 48 and the hydroperoxide intermediate. (a) Lys48 directly H-bonded to the hydroperoxide noted as ‘direct’). (b) One water molecule involed in the H-bond network between Lys48 and the hydroperoxide (noted as ‘one water bridge’). (c) Two interstitial water molecules between Lys48 and the hydroperoxide (noted as ‘two water bridge’). Colors for iron, sulfur, oxygen, nitrogen, carbon and hydrogen atoms are respectively ochre, yellow, red, blue, gray and white. Over the 100 ns simulation run, we extracted geometries every 12.5 ps generating 8000 structures. Each pattern was found in comparable statistical occurrence in chains A and B. Percentages of occurrence of the five different patterns** (Direct, one water bridge, two water

§

The corresponding figure for the SOR-superoxide system is given as SI, section II.1. The geometrical characteristic and cutoffs used to define the different patterns are described in the SI, section II.2.

**

ACS Paragon Plus Environment

21

Journal of Chemical Theory and Computation

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 22 of 73

bridge, water H-bonded to the hydroperoxide without Lys48 involvement, none) are shown in Table 3, including statistics for the two chains. These data clearly differentiate the interactions concerning the proximal and the distal oxygen atoms of the hydroperoxide intermediate. Table 3. Hydrogen bonds between the Fe2+-OO●- and Fe3+-OOH species and the lysine-water system. Op stands for the proximal oxygen and Od for the distal oxygen. SOR-Fe2+-OO●-

SOR-Fe3+-OOH

Op

Od

Op

Od

Direct

0.11%

1.97%

0.08%

0.21%

One water bridge

6.43%

18.97%

1.72%

1.18%

Two water bridge

11.12%

15.72%

18.11%

11.60%

Total involving Lys

18%

37%

20%

13%

Only H2O

38%

25%

56%

17%

None

44%

38%

24%

70%

SOR-Fe2+-OO●- intermediate. In the Fe2+-superoxide species, the distal oxygen is more prone to interact with the water molecule or the lysine than the proximal one. This is consistent with the fact that the protonation of this first reaction intermediate takes place on the distal oxygen (Scheme 1), as previously shown20. Since in this work we mainly focus in the protonation step of the second Fe3+-OOH intermediate, details for the Fe2+-OO●- species are reported in supporting information, section II.2.

ACS Paragon Plus Environment

22

Page 23 of 73

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

Journal of Chemical Theory and Computation

SOR-Fe3+-OOH intermediate. For the Fe3+-hydroperoxide species, the direct interaction between Lys48 with the OOH- moiety is very rare compared to the interactions involving one or two bridging molecules of water. Being protonated, the distal oxygen is much less involved as an acceptor of H-bonds than in the iron-superoxide state (70% of the conformations without H-bonds). On the other hand, the proximal oxygen is an acceptor of H-bonds with water molecules without the involvement of Lys48 for 56% of the conformations. When Lys48 is involved with or without water bridges, the picture is similar to the one obtained for the superoxide system. However, the proportion of ‘two water bridge’ is higher for the proximal oxygen of the hydroperoxide species (18%) when compared to the superoxide one (11%). The hydrogen carried by the distal oxygen is always H-bonded to a water molecule (97%), making it an effective donor. Overall, even if lysine 48 has a clear direct or indirect interaction with the OOH- intermediate, configurations not involving this residue clearly dominate. Beside the study of the H-bond network itself, an interesting correlation between the value of the dihedral angle Cβ-S-Op-Od and the occurrence of the ‘direct’ (lysine only) conformation was noted. While the dihedral angle Cβ-S-Op-Od can adopt two specific values in our simulations (Figure 5), ≈ 120° (b) and ≈ -50° (c), only one value around 128° (a) is found in X-Ray structures (2JI3)25. In the MM-MD simulations, this 120° conformation occurs only when lysine 48 is in direct H-bonding with the hydroperoxide intermediate, which is a very rare event (less than 0.2%). We can also note that the 120° conformation is more flexible (the amplitude is wider) than the -50° conformation.

ACS Paragon Plus Environment

23

Journal of Chemical Theory and Computation

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 24 of 73

Figure 5. Conformations of the Fe3+-OOH with different values of the Cβ-S-Op-Od dihedral angle. (a) From the X-Ray 2JI3 structure, with a dihedral Cβ-S-Op-Od angle equals to 128° (b) Simulation ‘direct’ lysine with a 120° Cβ-S-Op-Od dihedral angle. (c) Simulation ‘one water bridge’ lysine and with a -50° Cβ-S-Op-Od dihedral angle. The conformations are shown with a representaion of solvent accessible surface. In the 120° conformation (Fig. 5 b), the OOH- moiety is ‘outside’, making it easily accessible to the lysine for direct H-bonding. The arrangement of the environment around the lysine gives a solvent accessible surface with a narrower cavity. In the -50° conformation (Fig. 5 c), the OOH- moiety points towards the lysine making the proximal oxygen more easily accessible to the water molecules, with the lysine nearby to form a bridge. In this case, the cavity is more exposed to the solvent. Whereas the -50° conformation is the more preponderant one in our simulations, the 120° conformation observed in the X-ray structure might be stabilized by crystal packing forces92,93, favoring a more compact stacking. The 120° confirmation is still present in our work and with a conformation of the lysine in agreement with the X-Ray data.

ACS Paragon Plus Environment

24

Page 25 of 73

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

Journal of Chemical Theory and Computation

Spin state of the iron As reported in the literature, the iron present in the active site of SOR has a high spin configuration25,29. Hybrid functionals like B3LYP33–35 or PBE094,95 predict the correct spin state, whereas pure GGA functionals (BLYP33,34 or PBE96) favor low spin configuration22,23. One solution is to add a Hubbard correction potential97–99 to the iron to avoid overhybridization and so, obtain the high spin as the ground state and resulting correct spin densities. The cost of this approach is equivalent to the standard GGA whereas hybrid functionals are more costly because of the calculations of the Hartree-Fock exchange. With evolution of supercomputing power, QM/MM calculations with hybrid functionals are now more readily accessible100,101 and in this work we choose this method to avoid the need of a parameter for each different state102 of the system allowing us to ensure continuity. From previous work20,24, we knew that B3LYP was shown to predict the correct spin state and geometries for this system (for the QM model), so it was our choice. To assess the validity of this choice, we undergo for the Fe3+-OOH system geometries optimization at a QM/MM level, starting from geometry extracted from classical molecular dynamics. In this calculation, only the QM part was optimized, while the MM part was kept fixed. Consistently with results obtained on small QM models23,24 the hybrid B3LYP (and PBE0) functional is able to describe the system in the correct spin state (high spin, S=5/2) whereas the non-hybrid BLYP functional (as well as PBE) favors the low spin configuration (see Table 4). The geometrical properties of the Fe3+-OOH species (the only species considered in the following part) are reported hereafter. In the high spin state, the Fe-S bond distance was found to be 2.42 Å using B3LYP (B3LYP-HS) and 2.44 Å with BLYP (BLYP-HS). In the same species, the Fe-Op distances, has a value of 1.98 Å at the B3LYP-HS level and 2.01 Å at the BLYP-HS

ACS Paragon Plus Environment

25

Journal of Chemical Theory and Computation

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 26 of 73

level while the Op-Od distances were found to be 1.45 Å and 1.48 Å at the B3LYP-HS and BLYP-HS level, respectively. Unlike B3LYP (and PBE0), BLYP (and PBE) predicts the lowspin Fe3+-OOH species to be the most stable configuration. For the low spin configuration, BLYP (BLYP-LS) provides distances of 2.34 Å, 1.84 Å and 1.53 Å for the Fe-S, Fe-Op and Op-Od bonds respectively. This trend is consistent with the overhybridization of the d-orbitals of the iron, leading to a strengthening of the Fe-Op bond and a weakening of the Op-Od bond. As pointed out, PBE0 gave comparable results to B3LYP and PBE to BLYP. Table 4. Energies of the lower spin values systems relative to the high spin ones for the Fe3+OOH species. Energies are expressed in kcal mol-1. BLYP

PBE

B3LYP

PBE0

High Spin (Fe3+ 5/2)

0.0

0.0

0.0

0.0

Intermediate Spin (Fe3+ 3/2)

-4.5

-3.3

3.3

7.4

Low Spin (Fe3+ 1/2)

-13.9

-10.9

3.0

Conv. I.

To bridge the gap between purely classical simulations and metadynamics simulation, we ran a series of QM/MM MD simulations (5 ps each), starting from three MM configuration points selected to be representative of two H-bonding patterns as described earlier. (Snapshots ‘5050A’, chain A site, ‘one water bridge’; ‘5741B’, chain B, ‘direct’; ‘5497B’, chain B, ‘one water bridge’). The QM/MM simulations showed a rather good geometric stability of the active site, in continuity with classical MM dynamics described in the previous section††. These data confirm the validity of the force field parameters we have determined as well as the fact that the structure is equilibrated. We assumed that the ‘two water bridge’ is close to the ‘one water bridge’ in terms of mechanism.

††

Details can be found in the supporting information, section II.3.

ACS Paragon Plus Environment

26

Page 27 of 73

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

Journal of Chemical Theory and Computation

Reactivity of Fe3+-OOH via metadynamics As underlined in the Introduction, the protonation process of the Fe3+-OOH intermediate is thought to be the key step in the SOR catalytic cycle to form the H2O2 reaction product. Indeed, the SOR active site allows the formation of H2O2 by the specific protonation of the oxygen, initially bound to iron in the Fe3+-OOH species. It should be stressed out that the protonation of the distal oxygen of a Fe3+-OOH species might lead to a O-O bond cleavage and to the formation of high-valent iron-oxo Fe=O species, as it occurs during the catalytic cycle of cytochrome P450 oxygenases31. In this latter case, the formation of the high-valent iron-oxo species allows the catalysis of substrate oxidation and oxygenation. Thus, it is of critical relevance to provide a clear description of the protonation process of the Fe3+-OOH intermediate in order to understand the specificity of SOR as not being an oxygenase-type enzyme. Classical molecular dynamics simulations have shown that Lys48 interacts strongly, either directly, or via a H-bond network, with the hydroperoxide intermediate. In conjunction with experimental data20,25,30, our results thus have lead us to consider that this residue is the final proton source/relay in the protonation mechanism. To study the formation of the SOR reaction product H2O2, we have focused on the mechanism of proton transfer from the -NH3+ part of Lys48 to the OOH- moiety by the means of metadynamics simulations. Two starting conformations have been considered. A first one with the ‘direct’ conformation with the -NH3+ part of lysine 48 directly interacting with the hydroperoxide (i.e. lysine 48 is the direct source of proton) and a second one with one water molecule intercalated between the two (‘one water bridge’). This second pattern is believed to be representative of all those involving one or more molecules of water.

ACS Paragon Plus Environment

27

Journal of Chemical Theory and Computation

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 28 of 73

Between the two starting structures, apart from two variables, no major discrepancies were found between the structures: for example, the iron-sulfur bond distance was 2.49 ± 0.10 Å for the ‘direct’ conformation and 2.48 ± 0.08 Å for the ‘one water bridge’ conformation. The only noticeable differences were found for the iron-oxygen bond distance, which is slightly longer in the ‘direct’ conformation at 1.99 ± 0.05 Å versus 1.96 ± 0.05 Å for the ‘one water bridge’ conformation. Also, as discussed before, the dihedral angle Cβ-S-Op-Od is positive in the ‘direct’ conformation (70 ± 16°) and negative in the ‘one water bridge’ conformation (-40 ± 17°). For each simulation, two collective variables (CV) were used: one to describe the proton transfer and the other to describe the dissociation of the hydrogen peroxide product from the iron center (Fe-Op bond cleavage). In the system where Lys48 interacts directly with the OOH- moiety, the CV describing the proton transfer is the N-H distance with the closest H atom to the hydroperoxide (i.e. H-bonded). In the system with one molecule of water bridging the lysine and OOH-, the corresponding CV is the Owater-H distance. In both the cases, the second variable is the distance between the iron atom and the proximal oxygen bound to it. A first series of simulations were carried out using untempered metadynamics (UTMetaD) with coarse parameters in order to get an estimated barrier for the two starting conformations. These preliminary simulations have first confirmed the feasibility of the guessed process and the choice of the CV. They also provided a first estimate of the activation free energy of the protonation process, which was no more than 15 kcal mol-1. These results also enabled us to set the ∆T parameter for the well-tempered metadynamics (WTMetaD, presented hereafter) to a value of 6000 K.

ACS Paragon Plus Environment

28

Page 29 of 73

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

Journal of Chemical Theory and Computation

Direct protonation pathway. In this paragraph, the ‘direct’ protonation pathway of the Fe3+-OOH intermediate is detailed and the corresponding free energy surface (FES) reconstructed from the well-tempered metadynamics (WTMetaD) is shown in Figure 6. The metadynamics was started from the end of the QM/MM dynamics run (after 5 ps, position and velocities were conserved). It was stopped after 22 ps, corresponding to 2200 Gaussians added.

Figure 6. Free energy surface of the direct protonation pathway of the Fe3+-OOH intermediate with respects to the chosen CVs. The preferred reaction path is depicted as a continuous black line. An alternative path is shown in dashed line. The free energy reaction path found from the metadynamics simulations could be described in several steps from (AD) to the product (ED). Starting from the initial configuration (AD), taken from previous QM/MM molecular dynamics in which the Fe-Op distance is 1.9 Å, the system

ACS Paragon Plus Environment

29

Journal of Chemical Theory and Computation

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 30 of 73

evolves via a small activation barrier of about 1 kcal mol-1 towards a second minimum lying 1 kcal mol-1 below the initial one. This intermediate, labeled (A’D), is characterized by an increase of the Fe-O distance to 2.2 Å and might be stabilized by the nearby charged side-chain of the lysine. The second step consists in the further elongation of the Fe-Op distance. From (A’D), the system evolves towards (BD), an intermediate in which the Fe-Op bond distance is now 2.6 Å. Along these two first steps, the proton remains bound to the lysine (1.2 Å in intermediate BD). The second step occurs via a barrier of 7 kcal mol-1, which is the highest value found for this pathway. Using a small static DFT model, Attia et al103 have reported an energy cost of 17 kcal mol-1(internal energy) to separate the hydroperoxide from the Fe3+ center. The charged lysine side chain thus greatly helps the process of dissociation of the OOH- ion from the iron active site. Starting from (BD), two pathways are now possible, either via intermediate (CD) or intermediate (DD). In the path involving (CD), the protonation of Op occurs first and is followed by the complete release of H2O2. The other path, involving (DD) starts with the complete dissociation of the HOO- moiety followed by the protonation process (see Figure 6). Thus for the path through (CD), the protonation of the Op by the lysine occurs first. The N-H distance goes from 1.1 Å to 1.5 Å, while the Op-H distance decreases simultaneously from 1.9 Å to 1.1 Å. The associated transition barrier from (BD) to (CD) is 3.5 kcal mol-1. The system then evolves from (CD) with a slight barrier of 2 kcal mol-1, downhill towards the product well (ED), with the separated H2O2 moiety, the uncoordinated Fe3+ and the deprotonated lysine. The second path evolves from (BD) towards (DD) with a barrier of 4 kcal mol-1, so 1 kcal mol-1 higher than the previous path. It leads to a dissociation of the HOO- from the iron (Fe-Op is now 3.2 Å). The protonation occurs only after the Fe-Op dissociation has occurred, with a barrier of 1.4 kcal mol-1 and leads to the product (ED), similar to the one described above.

ACS Paragon Plus Environment

30

Page 31 of 73

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

Journal of Chemical Theory and Computation

The two paths lead to the same product-like state (ED), which is 5 kcal mol-1 above the initial reactants (A’D). The endergonicity of this reaction process is only in appearance since product (ED) is a reaction intermediate along the reaction path leading to the release of H2O2. Indeed, hydrogen peroxide will be dislodged and solvated while the ‘free’ sixth coordination of Fe3+ will eventually be further occupied by the Glu47 carboxylate side chain (Scheme 1). These processes then induce a stabilization and a re-organization of the surrounding of the active site, especially of the lysine loop. However these long-time scale processes are out of the scope of this work. Furthermore, the constraints used in these simulations restrict the explored subspace‡‡ and prevent the evolution toward the resting oxidized state, in which the sixth coordination position of the Fe3+ is occupied by the Glu47 carboxylate side chain (Scheme 1). Moreover, the CVs chosen are not fitted to describe the glutamate coordination state.

‡‡

The amplitude of the CVs was adequately restricted to ensure the convergence of the free energy through back and forth simulations between (AD) and (ED).

ACS Paragon Plus Environment

31

Journal of Chemical Theory and Computation

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 32 of 73

Figure 7. Structures of the minimum points of interests of the free energy surface along the ‘direct’ protonation pathway of the Fe3+-OOH species. The QM level is represented in ball and stick model (H atoms of amino acids residues were removed for clarity except for the ammonium group of the lysine 48), MM protein in the New Cartoon, MM solvent molecules are removed: (a) A’D, (b) BD, (c) CD, (d) DD, (e) ED. Figure 7 shows the evolution of the whole system from a geometric point of view instead of just the CVs. As the hydroperoxide uncoordinates from the iron and forms H2O2, the sulfur from Cys116 makes a tighter bond with the ferric iron going from 2.5 Å up to 2.2 Å. Also the deprotonated Lys48 slightly moves from 5.4 Å to 5.7 Å. However, the Lys48 motion is hindered, as one of the CVs is the N-H distance with a wall at 3 Å.

ACS Paragon Plus Environment

32

Page 33 of 73

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

Journal of Chemical Theory and Computation

One water bridge protonation pathway In this paragraph, the protonation pathway of the Fe3+-OOH intermediate involving one water molecule is studied and the free energy surface (FES) reconstructed from the WTMetaD (Figure 8). The metadynamics was started from the end of the QM/MM dynamics run (after 5 ps, positions and velocities were conserved). It was stopped after 15 ps corresponding to 1500 Gaussians added.

Figure 8. Free energy surface of the ‘one water bridge’ protonation pathway of the Fe3+-OOH intermediate, respective to the chosen CV. From a mechanistic point of view, this pathway is different from the direct pathway (vide supra). This can be observed from the space above the starting point of the simulation (Aw) which is clearly visited, contrary to the ‘direct’ pathway. Thus, the first step of this pathway is clearly the protonation of the hydroperoxide intermediate, with a very small Fe-Op elongation.

ACS Paragon Plus Environment

33

Journal of Chemical Theory and Computation

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 34 of 73

Going from (AW) to (Bw), which displays bond distances for H-OW of 1.4 Å, H-Op of 1.2 Å and a Fe-Op of 2.05 Å, involves a barrier around 11 kcal mol-1. The (Bw) intermediate evolves with a barrier of 2 kcal mol-1 into (Cw), which displays bond distances for H-Ow of 1.6 Å, H-Op of 1.0 Å and Fe-Op of 2.2 Å1. Also the lysine protonates the hydroxide ion to form back H2O. Then (Cw) evolves towards (Dw) with a transition barrier of 1 kcal mol-1. The Fe-Op distance elongates further to 2.4 Å and the H2O2 is starting to get dissociated. Several intermediates with different H-Ow and Fe-Op distances are found and we were able to carry out back and forth simulations from (Dw) to (Aw). The transition (Aw) to (Dw) describes partly the proton hopping between the lysine and the Fe3+-hydroperoxide via a water molecule. However, the proton transfer from this water molecule to the surrounding water molecules prevents us to have a clearly defined product by the used CV and gives this scattered surface around (Ew). In spite of all this, the estimation of the barrier is still valid since the subsequent proton transfer does not influence it. We report the geometries in Figure 9. The ‘one water bridge’ pathway shows a free energy difference of an extra 4 kcal mol-1 compared to the first rate-limiting step of the direct pathway.

ACS Paragon Plus Environment

34

Page 35 of 73

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

Journal of Chemical Theory and Computation

Figure 9. Structures of the minimum points of interests of the free energy surface of the ‘one water bridge’ protonation pathway of the Fe3+-OOH intermediate. The QM level is represented in ball and stick model (H atoms of amino acids residues were removed for clarity except for the ammonium group of the Lysine), MM protein in the New Cartoon, MM solvent molecules are removed: (a) AW, (b) BW, (c) CW, (d) DW, (e) EW.

The protonation is the first step for the ‘one water bridge’ pathway. The hydroperoxide is still bounded to the iron and get protonated by the water molecule (BW), which in turn gets very quickly protonated by the lysine (CW, 2-5 fs). No minimum was found for a Fe3+-H2O2--OH---NH3+-Lys48 intermediate. In this case, the lysine arm is more mobile as it is not

ACS Paragon Plus Environment

35

Journal of Chemical Theory and Computation

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 36 of 73

hindered by any CVs. It can move§§ from 6.8 Å to 7.3 Å from the iron unlike, the water molecule involved in the protonation for which the Ow-Hw distance is limited by the metadynamic wall on the CV at 3 Å. The sulfur atom follows the same trend as the direct pathway: as soon as H2O2 is formed (EW) and gets dissociated, the Fe-S bond is strengthened.

Reaction mechanism As described above, for the protonation of the Fe3+-OOH intermediate two different pathways have been identified. i) The ‘direct’ mechanism involving at first the dissociation of the OOH- moiety from the iron, followed by the protonation of the proximal oxygen. ii) The ‘one water bridge’ mechanism, with the protonation of the proximal oxygen occurring first, followed by the dissociation of the H2O2 from the iron. We assume that the ‘two water bridge’ pathway will follow the same process than the ‘one water bridge’. In both cases, the first step of the pathway is always the limiting one, whichever it concerns the dissociation (‘direct’) or the protonation (‘one water bridge’) process. The ‘direct’ process has a free energy activation barrier of 7 kcal mol-1 which is lower by 4 kcal mol-1 compared the ‘one water bridge’ step (11 kcal mol-1). Despite the fact that the ‘direct’ conformation is less probable than the ‘water bridged’ ones by a factor of 100, vide supra, this 4 kcal mol-1 energy difference could be enough to make the direct pathway the most favored one. In fact, using the Eyring-Polanyi equation and assuming ∆F‡ ≈ ∆G‡ + C (with C = -PV which is the expansion work and cancel out between the two mechanisms), we can estimate rate constants of 5 107 s-1 for the ‘direct’ mechanism and 6 104 s-1 for ‘one water bridge’ mechanism. These data strongly support the fact that the ‘direct’ mechanism associated with the dissociation of the

§§

The lysine is found to be in a distance from 6.6 Å to 8 Å when the proton hopping has happened.

ACS Paragon Plus Environment

36

Page 37 of 73

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

Journal of Chemical Theory and Computation

OOH- in a first step followed by its protonation is the major, if not the sole, process which occurs during SOR catalysis. This direct pathway highlights an unprecedented dissociation mechanism of the OOH- moiety in a Fe3+OOH intermediate, occurring without being directly associated with the protonation of the proximal oxygen29. This reactivity favoring in a first step the dissociation of the OOH- moiety from the iron was unexpected since previous studies using small quantum models (without Lys48) to describe the same process exhibited an energetical cost of 17 kcal mol-1

103

. In order to explain if this

difference in reactivity could be due to the presence of the lysine in the proximity of the iron-hydroperoxide species, we carried out several calculations using different models and levels of theory. In particular, we used Mulliken charges, directly available from QM and QM/MM calculations to provide a qualitative picture of the effects of environment on the electron structure of the systems. Since Mulliken charges are sensible to the method used (especially basis set), complementary calculations of density concentration, based on AIM partition,104,105 were also carried to assess the influence of the environment on the Fe-O bond strength. Results indicate a weakening of the Fe-O bond in presence of the environment (see details in SI, section II.5). In the last 2 ps of the QM/MM dynamics, the charge analysis showed a significant charge difference on the proximal oxygen of the Fe3+-OOH intermediate: -0.48 ± 0.08 e for the ‘one water bridge’ system and -0.28 ± 0.06 e for the ‘direct’ system. On the other hand, the distal oxygen showed a less pronounced difference in charge between both systems: -0.53 ± 0.05 e for ‘one water bridge’ and -0.45 ± 0.05 e for ‘direct’ process. In other words, the ‘direct’ system showed a marked difference of charge between the distal and the proximal oxygens, whereas the charges of the two oxygens from the ‘one water bridge’ were more similar. It should be noted that the charges on the iron atom is also in pair with the charge of the proximal oxygen of the

ACS Paragon Plus Environment

37

Journal of Chemical Theory and Computation

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 38 of 73

hydroperoxide moiety. The iron atom is more negative by 0.50 e in the ‘direct’ system compared to the ‘one water bridge’. In order to understand if this difference can be explained by the close proximity of the lysine towards the hydroperoxide solute in the ‘direct’ system (less than 3 Å), we underwent calculations on small models. The models consisted of the iron atom, the cysteine and the four histidines ligand, and the hydroperoxide adduct (called hereafter AS for Active Site) at the B3LYP level using the same basis set as the QM/MM on respective atoms. Then, we added to this model the explicit QM solvent water molecules and/or the lysine also at the QM level. Geometries were extracted from various points of the ‘direct’ QM/MM dynamic and single point calculations were performed with both CP2K and Gaussian09 program packages. We observed that the Mulliken charges on this small model without the environment, even if they come from the direct QM/MM simulations, are closer to the charges observed for the one water mechanism (vide supra), i.e. the proximal oxygen is more negative than the distal one (AS, Table 5). Addition of the lysine, which is a positively charge moiety close to active site, does not reverse the charge population (AS with lysine, Table 5). The addition of the QM solvent water molecules only, leads to the same result as the addition of the lysine, namely an equilibrated charge distribution between the two oxygens (AS with explicit water, Table 5). Interestingly, the presence of both the lysine and the solvent gave results close to the full QM/MM level (AS with explicit water and lysine, Table 5). This tends to indicate that the solvent environment amplifies the role of the lysine. This effect can be seen as local. The charges are not much impacted by the bulk of the solvent or by the protein but mostly by the proximity of the charge of the lysine (NH3+) and the hydrogen bond network involving the water molecules, the lysine and the solute.

ACS Paragon Plus Environment

38

Page 39 of 73

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

Journal of Chemical Theory and Computation

Despite being due to a local surrounding (microsolvation, neighboring residues), the polarization effect described above is a consequence of the whole environment and its dynamics. Therefore, small static quantum models can hardly describe it properly.

Table 5. Values of the Mulliken charges for the 2 oxygen atoms of the Fe3+-OOH intermediate in the ‘direct’ protonation pathway for various models in gas phase. Values were extracted from CP2K software. In italic, values were computed with Gaussian09 software. The values in parenthesis are the mean values. AS

AS with lysine

AS with explicit water

AS with explicit water and lysine

Full system QM/MM

qOp (Proximal oxygen) (e)

-0.48 / -0.48

-0.40 / -0.41

-0.40 / -0.41

-0.28 / -0.29

-0.24 (-0.28)

qOd (Distal oxygen) (e)

-0.37 / -0.37

-0.40 / -0.41

-0.40 / -0.41

-0.45 / -0.47

-0.44 (-0.45)

In Table 6, we studied the influence of the use of implicit solvent compared to Table 5 by including in all the calculations the implicit solvation (IEF-PCM, water) as implemented in Gaussian09. Replacing the water by implicit solvent does not contribute to gain on the charge distribution with the distal oxygen as the most negative oxygen atom. This conclusion is even true when the lysine is included in the system. In this case, the charge population between the two oxygen atoms is the same, exactly like in the gas phase. Altogether, these data highlight the importance of hydrogen bonding, and not only a bulk dielectric media, i.e. the micro solvation effects and the way it amplifies the effect of the charge near the active site in enzyme reactivity. Also, our data clearly underline the limitation of implicit models to describe this effect, which once again emphasized the need of a complete

ACS Paragon Plus Environment

39

Journal of Chemical Theory and Computation

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 40 of 73

model to get a proper description of the reactivity both for the intrinsic mechanism as well as the energetic features. Table 6. Values of the Mulliken charges for both hydroperoxide oxygens in the ‘direct’ configuration for various models with implicit solvation in G09. AS

AS with lysine

AS with explicit water

AS with explicit water and lysine

qOp (Proximal oxygen) (e)

-0.51

-0.39

-0.46

-0.26

qOd (Distal oxygen) (e)

-0.38

-0.41

-0.46

-0.46

The same types of calculations were carried out with the ‘one water bridge’ system and details can be found in the SI, section II.5. It is important to note that in the latter case the influence of the lysine is weaker than the direct interaction. The QM/MM charge distribution of both oxygens is already closer to that of the active site model. The charges are fairly equilibrated between both oxygen, less than 0.10 e compared to the 0.20 e difference in the ‘direct’ system.

After studying the charge population on the starting geometry, we focused on the evolution of the charge along the metadynamics. For the ‘direct’ system the charge of the proximal oxygen evolves from -0.28 e in the Fe3+-OOH species to -0.45 e in the dissociated form Fe3+ + OOH- After the protonation has occurred, both oxygens acquire the same charge value of around -0.38 e as found in hydrogen peroxide alone (Calculations were made for the hydrogen peroxide alone at the same level of theory). For the ‘one water bridge’ system, the protonation is occurring first. During this first step, the charge on the proximal oxygen that is already close to the value of the hydrogen peroxide goes from -0.50 e to -0.40 e. Then dissociation of H2O2 slightly modifies the charges towards -0.38 e as for the ‘direct’ system.

ACS Paragon Plus Environment

40

Page 41 of 73

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

Journal of Chemical Theory and Computation

From these results, we can propose that the dissociation occurs first for the ‘direct’ because of the impact of the lysine in close proximity and the solvation by water molecule around it. This makes the proximal oxygen less negatively charged when bounded to the iron and less susceptible to direct protonation. The equilibrium is furthermore displaced towards the Fe3+ + OOH- dissociated form assisted by the lysine close enough (and its NH3+ group) to stabilize the OOH- ion. Upon the dissociation, the proximal oxygen reacquires the negative charge surplus, which easily triggers the subsequent proton transfer on the proximal oxygen (less than 4 kcal mol-1), giving the H2O2 reaction product. A contrario, the ‘one water bridge’ configuration with the lysine further away doesn’t have any impact on the charge of the proximal oxygen and the equilibrium is not displaced towards the dissociated form Fe3+ + OOH-. The protonation occurs first, which in return promotes the dissociation of the H2O2 reaction product from the iron. The free energy values in association with the charges analysis emphasize the role of lysine 48. Indeed, lysine 48 appears essential as a proton donor and also has a major impact on the charges of the proximal oxygen. Coupled with the micro solvation of water around the active site, it favors the mechanism in which the dissociation of the hydroperoxo moiety from the iron atom occurs first, followed by its protonation and the formation of the hydrogen peroxide reaction product. Also without direct hydrogen bonding with the solute, the lysine still has a role in the protonation of the proximal oxygen by a water molecule as a subsequent proton source. This leads to the second mechanism, only in the case of an indirect bonding of the lysine, consisting of the protonation first and then the dissociation of the H2O2 product. In this case, water could also be a proton donor even if the cost in an aqueous environment is expected to be slightly higher.

ACS Paragon Plus Environment

41

Journal of Chemical Theory and Computation

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 42 of 73

Conclusion The ferric iron hydroperoxo (Fe3+-OOH) species is a key reaction intermediate in the catalytic cycle of Superoxide Reductase. Upon the protonation of its proximal oxygen, the Fe-O bond is cleaved and the final product H2O2 is released from the active site. Despite being a determining step along the catalytic cycle, this second protonation is a poorly experimentally documented process. A complete simulation workflow has thus been designed to understand how the protonation of the proximal oxygen of the Fe3+-OOH intermediate is efficiently catalyzed by SOR, hence giving its specificity as a detoxifying enzyme and not as an oxygenase (such as P450, with which SOR shares a similar Fe3+-OOH intermediate). This simulation workflow involves a MM-MD conformational study of the iron-hydroperoxide (Fe3+-OOH) key intermediate, followed by a QM/MM well-tempered metadynamic study of the protonation process. This state-of-the-art approach has allowed us to highlight two possible mechanisms of formation of the final product H2O2. In the first one, the hydroperoxide moiety dissociates from the iron center prior to its protonation on the proximal oxygen. This rather unexpected mechanism has a rate-limiting barrier of only 7 kcal mol-1 and involves an initial geometry with a direct interaction between the hydroperoxo ligand and the lysine 48. The second mechanism found occurs when the starting structure involves interstitial molecules of water between the lysine and the hydroperoxo moiety. In this case, the first step is the protonation, with a rate-limiting barrier of 11 kcal mol-1 followed by the dissociation of H2O2 from the iron. From this study, we have drawn not only mechanistic conclusions with an uncommon first rate-limiting step but we also have highlighted the need for a precise description of the

ACS Paragon Plus Environment

42

Page 43 of 73

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

Journal of Chemical Theory and Computation

environment as well as developed a protocol adapted to SOR metalloenzymes, including the development of parameters for iron active site. The MM parameters developed for this study are very consistent and robust across long simulation times (100 ns) for the reaction intermediate Fe3+ OOH. The subsequently ran MM dynamics have allowed to sample the conformations of the enzyme, especially around the active site. We particularly found that the specific micro-solvation patterns, leading to direct or water-mediated interactions of the substrate with its nearby proteic environment was of critical relevance in the protonation mechanism. Indeed, the hydroperoxide substrate actually interacts with the side chain of lysine 48 in three distinct ways: either directly, or with the involvement of one or two interstitial molecules of water. The MM-MD simulations have also showed the high flexibility of the side chain of Lys48, which supports the hypothesis for its role as a proton source/relay during the second protonation step in the catalytic cycle of SOR. During QM/MM well-tempered metadynamics simulations, the hybrid density functional B3LYP has been used since it accurately accounts for the specific high spin configuration of the metallic active site. To our knowledge, no similar method, coupling hybrid DFT with QM/MM metadynamics, has ever been documented to date for such a metallo-enzymatic system. The QM/MM well-tempered metadynamics approach also allowed accounting for the environment effects on the electronic structure of the active site, especially interactions with lysine 48 and the micro solvation structures. Environment plays a critical role since it determines the protonation mechanism. Indeed, depending on the interaction pattern of lysine 48 with the substrate, two competing mechanisms were found. The lowest rate-limiting barrier (7 kcal mol-1) was found for the ‘direct’ mechanism (dissociation first, followed by protonation), which therefore seems to be the preferred mechanism for the second protonation in the catalytic cycle of SOR. The less

ACS Paragon Plus Environment

43

Journal of Chemical Theory and Computation

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 44 of 73

favorable second mechanism (protonation followed by dissociation) was found when one molecule of water is intercalated between the substrate and lysine 48. Furthermore, by lowering the activation barriers along the protonation mechanism, environmental effects caused by both charged ammonium group of lysine 48 and the micro-solvation structure, directly impose the preferred mechanism, which could not have been described with a QM-only model. Our work, based on a state-of-the-art methodology, has thus lead to the elucidation of the mechanism by which the formation of the H2O2 product is catalyzed by SOR from the key intermediate Fe3+-OOH. Our results have also evidenced the critical role of the environment of the active site (lysine 48 and solvent structure) in enhancing the second protonation reaction, thus making SOR a very efficient enzymatic system. Site directed mutagenesis experiments on SOR have shown that, upon the replacement of lysine 48 by an isoleucine (K48I mutant), the enzyme undergoes a dramatic alteration of its functioning. In this mutant, although upon reaction with O2●- the protonation of the first reaction intermediate (Fe2+-OO●-, Scheme 1) is not affected compared to the wild type, H2O2 is no more produced. Under different experimental conditions, it was shown that the Fe3+-OOH intermediated in the K48I mutant leads to the formation of a high-valent iron-oxo Fe=O species, which could not be observed for the wild-type SOR. This mutation has attracted attention on SOR since it appeared possible to finely tune its properties and even expect turning this reductase into an oxygenase. Thanks to the information our methodological approach is able to provide, we have now started to study the properties and reactivity of K48I and other mutants of SOR to understand the different factors affecting the reactive processes.

ACS Paragon Plus Environment

44

Page 45 of 73

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

Journal of Chemical Theory and Computation

Finally, it is important to underline that the whole protocol presented here, from MM parameters determination to dynamic modeling of the reactive process, is transferable to any enzymatic system.

ASSOCIATED CONTENT Supporting Information Additional details about computational details (FF parameters models and forces) and resultst (validity of the models, more description of the H-bonding network, QM/MM dynamics and metadynamics). FF charges and forces in the Amber suite format are provided in an archive. The Supporting Information is available free of charge on the ACS Publications website.

AUTHOR INFORMATION Corresponding Author *Email: [email protected]

ACKNOWLEDGEMENTS The authors thank the LabEx ARCANE (ANR-11-LABX-0003-01) for funding the RD grant. All the computations presented in this paper were performed using the Froggy platform of the CIMENT infrastructure, which is supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA) and the Equip@Meso project (ANR-10-EQPX-29-01) and the Ceciccluster platform of the PCECIC infrastructure, which is supported by ICMG (FR 2607). We also thank Dr Pierre Girard for his technical expertise and support.

ACS Paragon Plus Environment

45

Journal of Chemical Theory and Computation

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 46 of 73

ABBREVIATIONS SOR, superoxide reductase; MM, Molecular mechanics; QM, Quantum mechanics; QM/MM, Quantum mechanics/Molecular mechanics; ROS, reactive oxygen species; SOD, superoxide dismutase; DFT, Density Functional Theory; MCPB, Metal Center Parameter Builder; FF, Force Field; RESP, Restrained electrostatic potential; DX, desulforedoxin; GEEP, Gaussian expansion of the electrostatic potential; GAPW, Gaussian augmented plane wave; RMSD, Root mean square deviation; MD, Molecular dynamic; GGA, Gradient generalized approximation; CV, Collective variable; FES, Free energy surface; IEFPCM, Integral equation formalism polarizable continuum model.

REFERENCES (1)

Winterbourn, C. C. Reconciling the Chemistry and Biology of Reactive Oxygen Species. Nat. Chem. Biol. 2008, 4 (5), 278–286.

(2)

Sabharwal, S. S.; Schumacker, P. T. Mitochondrial ROS in Cancer: Initiators, Amplifiers or an Achilles’ Heel? Nat. Rev. Cancer 2014, 14 (11), 709–721.

(3)

Imlay, J. A. Cellular Defenses against Superoxide and Hydrogen Peroxide. Annu. Rev. Biochem. 2008, 77 (1), 755–776.

(4)

Jenney Jr., F. E.; Verhagen, M. F. J. M.; Cui, X.; Adams, M. W. W. Anaerobic Microbes: Oxygen Detoxification Without Superoxide Dismutase. Science (80-. ). 1999, 286 (5438), 306–309.

(5)

Lombard, M.; Fontecave, M.; Touati, D.; Niviere, V.; Nivière, V.; Niviere, V. Reaction of the Desulfoferrodoxin from Desulfoarculus Baarsii with Superoxide Anion: EVIDENCE

ACS Paragon Plus Environment

46

Page 47 of 73

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

Journal of Chemical Theory and Computation

FOR A SUPEROXIDE REDUCTASE ACTIVITY. J. Biol. Chem. 2000, 275 (1), 115– 121. (6)

Kurtz, D. M. Microbial Detoxification of Superoxide: The Non-Heme Iron Reductive Paradigm for Combating Oxidative Stress. Acc. Chem. Res. 2004, 37 (11), 902–908.

(7)

Nivière, V.; Fontecave, M. Discovery of Superoxide Reductase: An Historical Perspective. J. Biol. Inorg. Chem. 2004, 9 (2), 119–123.

(8)

Pereira, A. S.; Tavares, P.; Folgosa, F.; Almeida, R. M.; Moura, I.; Moura, J. J. G. Superoxide Reductases. Eur. J. Inorg. Chem. 2007, 2007 (18), 2569–2581.

(9)

Sheng, Y.; Abreu, I. a.; Cabelli, D. E.; Maroney, M. J.; Miller, A.-F.; Teixeira, M.; Valentine, J. S. Superoxide Dismutases and Superoxide Reductases. Chem. Rev. 2014, 114 (7), 3854–3918.

(10)

Emerson, J. P.; Coulter, E. D.; Phillips, R. S.; Kurtz, D. M. Kinetics of the Superoxide Reductase Catalytic Cycle. J. Biol. Chem. 2003, 278 (41), 39662–39668.

(11)

Adam, V.; Royant, A.; Nivière, V.; Molina-Heredia, F. P.; Bourgeois, D. Structure of Superoxide Reductase Bound to Ferrocyanide and Active Site Expansion upon X-RayInduced Photo-Reduction. Structure 2004, 12 (9), 1729–1740.

(12)

Archer, M.; Huber, R.; Tavares, P.; Moura, I.; Moura, J. J. G.; Carrondo, M. A.; Sieker, L. C.; LeGall, J.; Romão, M. J.; Romco, M. J. Crystal Structure of Desulforedoxin from Desulfovibrio Gigas Determined at 1.8 A Resolution: A Novel Non-Heme Iron Protein Structure. J. Mol. Biol. 1995, 251 (5), 690–702.

(13)

Bonnot, F.; Duval, S.; Lombard, M.; Valton, J.; Houée-Levin, C.; Nivière, V.

ACS Paragon Plus Environment

47

Journal of Chemical Theory and Computation

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 48 of 73

Intermolecular Electron Transfer in Two-Iron Superoxide Reductase: A Putative Role for the Desulforedoxin Center as an Electron Donor to the Iron Active Site. JBIC J. Biol. Inorg. Chem. 2011, 16 (6), 889–898. (14)

Folgosa, F.; Cordas, C. M.; Santos, J. A.; Pereira, A. S.; Moura, J. J. G.; Tavares, P.; Moura, I. New Spectroscopic and Electrochemical Insights on a Class I Superoxide Reductase: Evidence for an Intramolecular Electron-Transfer Pathway. Biochem. J. 2011, 438 (3), 485–494.

(15)

Horch, M.; Utesch, T.; Hildebrandt, P.; Mroginski, M. A.; Zebger, I. Domain Motions and Electron Transfer Dynamics in 2Fe-Superoxide Reductase. Phys. Chem. Chem. Phys. 2016, 18 (33), 23053–23066.

(16)

Emerson, J. P.; Coulter, E. D.; Cabelli, D. E.; Phillips, R. S.; Kurtz, D. M. Kinetics and Mechanism of Superoxide Reduction by Two-Iron Superoxide Reductase from Desulfovibrio Vulgaris. Biochemistry 2002, 41 (13), 4348–4357.

(17)

Nivière, V.; Asso, M.; Weill, C. O.; Lombard, M.; Guigliarelli, B.; Favaudon, V.; HouéeLevin, C. Superoxide Reductase from Desulfoarculus Baarsii : Identification of Protonation Steps in the Enzymatic Mechanism. Biochemistry 2004, 43 (3), 808–818.

(18)

Rodrigues, J. V.; Abreu, I. A.; Cabelli, D.; Teixeira, M. Superoxide Reduction Mechanism of Archaeoglobus Fulgidus One-Iron Superoxide Reductase. Biochemistry 2006, 45 (30), 9266–9278.

(19)

Huang, V. W.; Emerson, J. P.; Kurtz, D. M. Reaction of Desulfovibrio Vulgaris Two-Iron Superoxide Reductase with Superoxide: Insights from Stopped-Flow Spectrophotometry.

ACS Paragon Plus Environment

48

Page 49 of 73

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

Journal of Chemical Theory and Computation

Biochemistry 2007, 46 (40), 11342–11351. (20)

Bonnot, F.; Molle, T.; Ménage, S.; Moreau, Y.; Duval, S.; Favaudon, V.; Houée-Levin, C.; Nivière, V. Control of the Evolution of Iron Peroxide Intermediate in Superoxide Reductase from Desulfoarculus Baarsii. Involvement of Lysine 48 in Protonation. J. Am. Chem. Soc. 2012, 134 (11), 5120–5130.

(21)

Silaghi-Dumitrescu, R.; Silaghi-Dumitrescu, I.; Coulter, E. D.; Kurtz, D. M. Computational Study of the Non-Heme Iron Active Site in Superoxide Reductase and Its Reaction with Superoxide. Inorg. Chem. 2003, 42 (2), 446–456.

(22)

Dey, A.; Jenney, F. E.; Adams, M. W. W.; Johnson, M. K.; Hodgson, K. O.; Hedman, B.; Solomon, E. I.; Jenney, F. E.; Adams, M. W. W.; Johnson, M. K.; Hodgson, K. O.; Hedman, B.; Solomon, E. I. Sulfur K-Edge X-Ray Absorption Spectroscopy and Density Functional Theory Calculations on Superoxide Reductase: Role of the Axial Thiolate in Reactivity. J. Am. Chem. Soc. 2007, 129 (41), 12418–12431.

(23)

Sit, P. H.-L.; Migliore, A.; Ho, M.-H.; Klein, M. L. Quantum Mechanical and Quantum Mechanical/Molecular Mechanical Studies of the Iron−Dioxygen Intermediates and Proton Transfer in Superoxide Reductase. J. Chem. Theory Comput. 2010, 6 (9), 2896– 2909.

(24)

Tremey, E.; Bonnot, F.; Moreau, Y.; Berthomieu, C.; Desbois, A.; Favaudon, V.; Blondin, G.; Houée-Levin, C.; Nivière, V. Hydrogen Bonding to the Cysteine Ligand of Superoxide Reductase: Acid-Base Control of the Reaction Intermediates. J. Biol. Inorg. Chem. 2013, 18 (7), 815–830.

ACS Paragon Plus Environment

49

Journal of Chemical Theory and Computation

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

(25)

Page 50 of 73

Katona, G.; Carpentier, P.; Niviere, V.; Amara, P.; Adam, V.; Ohana, J.; Tsanov, N.; Bourgeois, D. Raman-Assisted Crystallography Reveals End-On Peroxide Intermediates in a Nonheme Iron Enzyme. Science (80-. ). 2007, 316 (5823), 449–453.

(26)

Mathé, C.; Mattioli, T. A.; Horner, O.; Lombard, M.; Latour, J.-M.; Fontecave, M.; Nivière, V. Identification of Iron(III) Peroxo Species in the Active Site of the Superoxide Reductase SOR from Desulfoarculus Baarsii. J. Am. Chem. Soc. 2002, 124 (18), 4966– 4967.

(27)

Mathé, C.; Nivière, V.; Houée-Levin, C.; Mattioli, T. A. Fe3+–η2–peroxo Species in Superoxide Reductase from Treponema Pallidum. Comparison with Desulfoarculus Baarsii. Biophys. Chem. 2006, 119 (1), 38–48.

(28)

Mathe, C.; Weill, C. O.; Mattioli, T. A.; Berthomieu, C.; Houee-Levin, C.; Tremey, E.; Niviere, V. Assessing the Role of the Active-Site Cysteine Ligand in the Superoxide Reductase from Desulfoarculus Baarsii. J. Biol. Chem. 2007, 282 (30), 22207–22216.

(29)

Horner, O.; Mouesca, J.-M.; Oddou, J.-L.; Jeandey, C.; Nivière, V.; Mattioli, T. A.; Mathé, C.; Fontecave, M.; Maldivi, P.; Bonville, P.; Halfen, J. A.; Latour, J.-M. Mössbauer Characterization of an Unusual High-Spin Side-On Peroxo−Fe 3+ Species in the Active Site of Superoxide Reductase from Desulfoarculus Baarsii . Density Functional Calculations on Related Models †. Biochemistry 2004, 43 (27), 8815–8825.

(30)

Bonnot, F.; Tremey, E.; von Stetten, D.; Rat, S.; Duval, S.; Carpentier, P.; Clemancey, M.; Desbois, A.; Nivière, V. Formation of High-Valent Iron-Oxo Species in Superoxide Reductase: Characterization by Resonance Raman Spectroscopy. Angew. Chemie Int. Ed. 2014, 53 (23), 5926–5930.

ACS Paragon Plus Environment

50

Page 51 of 73

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

Journal of Chemical Theory and Computation

(31)

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.

(32)

Dawson, J. H.; Holm, R. H.; Trudell, J. R.; Barth, G.; Linder, R. E.; Bunnenberg, E.; Djerassi, C.; Tang, S. C. Magnetic Circular Dichroism Studies. 43. Oxidized Cytochrome P-450. Magnetic Circular Dichroism Evidence for Thiolate Ligation in the SubstrateBound Form. Implications for the Catalytic Mechanism. J. Am. Chem. Soc. 1976, 98 (12), 3707–3709.

(33)

Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38 (6), 3098–3100.

(34)

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 (2), 785–789.

(35)

Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98 (7), 5648.

(36)

Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99 (20), 12562–12566.

(37)

Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100 (2), 1–4.

(38)

Laio, A.; Gervasio, F. L. Metadynamics: A Method to Simulate Rare Events and Reconstruct the Free Energy in Biophysics, Chemistry and Material Science. Reports Prog. Phys. 2008, 71 (12), 126601.

ACS Paragon Plus Environment

51

Journal of Chemical Theory and Computation

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

(39)

Page 52 of 73

Michel, C.; Laio, A.; Mohamed, F.; Krack, M.; Parrinello, M.; Milet, A. Free Energy Ab Initio Metadynamics: A New Tool for the Theoretical Study of Organometallic Reactivity? Example of the C−C and C−H Reductive Eliminations from Platinum(IV) Complexes. Organometallics 2007, 26 (5), 1241–1249.

(40)

Ensing, B.; De Vivo, M.; Liu, Z.; Moore, P.; Klein, M. L. Metadynamics as a Tool for Exploring Free Energy Landscapes of Chemical Reactions. Acc. Chem. Res. 2006, 39 (2), 73–81.

(41)

Blumberger, J.; Ensing, B.; Klein, M. L. Formamide Hydrolysis in Alkaline Aqueous Solution: Insight from Ab Initio Metadynamics Calculations. Angew. Chemie Int. Ed. 2006, 45 (18), 2893–2897.

(42)

Ensing, B.; Klein, M. L. Perspective on the Reactions between F- and CH3CH2F: The Free Energy Landscape of the E2 and SN2 Reaction Channels. Proc. Natl. Acad. Sci. 2005, 102 (19), 6755–6759.

(43)

Sgrignani, J.; Iannuzzi, M.; Magistrato, A. On the Role of Water in The Puzzling Mechanism of Final Aromatization Step Promoted by the Human Aromatase Enzyme. Insights from QM/MM MD Simulations. J. Chem. Inf. Model. 2015, 151002105138002.

(44)

Gouron, A.; Milet, A.; Jamet, H. Conformational Flexibility of Human Casein Kinase Catalytic Subunit Explored by Metadynamics. Biophys. J. 2014, 106 (5), 1134–1141.

(45)

Biarnés, X.; Ardèvol, A.; Iglesias-Fernández, J.; Planas, A.; Rovira, C. Catalytic Itinerary in 1,3-1,4-β-Glucanase Unraveled by QM/MM Metadynamics. Charge Is Not Yet Fully Developed at the Oxocarbenium Ion-like Transition State. J. Am. Chem. Soc. 2011, 133

ACS Paragon Plus Environment

52

Page 53 of 73

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

Journal of Chemical Theory and Computation

(50), 20301–20309. (46)

Alfonso-Prieto, M.; Biarnés, X.; Vidossich, P.; Rovira, C. The Molecular Mechanism of the Catalase Reaction. J. Am. Chem. Soc. 2009, 131 (33), 11751–11761.

(47)

Peters, M. B.; Yang, Y.; Wang, B.; Füsti-Molnár, L.; Weaver, M. N.; Merz, K. M. Structural Survey of Zinc-Containing Proteins and Development of the Zinc AMBER Force Field (ZAFF). J. Chem. Theory Comput. 2010, 6, 2935–2947.

(48)

Case, D. A.; Darden, T. A.; Cheatham, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Swails, J.; Goetz, A. W.; Kolossváry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wolf, R. M.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M. J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; SalomonFerrer, R.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER 12. University of California, San Francisco 2012.

(49)

Seminario, J. M. Calculation of Intramolecular Force Fields from Second-Derivative Tensors. Int. J. Quantum Chem. 1996, 60 (7), 1271–1277.

(50)

Hoops, S. C.; Anderson, K. W.; Merz, K. M. Force Field Design for Metalloproteins. J. Am. Chem. Soc. 1991, 113 (22), 8262–8270.

(51)

Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. a. A Well-Behaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97 (40), 10269–10280.

(52)

Besler, B. H.; Merz, K. M.; Kollman, P. A. Atomic Charges Derived from Semiempirical

ACS Paragon Plus Environment

53

Journal of Chemical Theory and Computation

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 54 of 73

Methods. J. Comput. Chem. 1990, 11 (4), 431–439. (53)

Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, Jr., J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03 Revision C.02. Gaussian Inc. Wallingford CT 2004.

(54)

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.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., J.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.;

ACS Paragon Plus Environment

54

Page 55 of 73

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

Journal of Chemical Theory and Computation

Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09 Revision D.01. Gaussian Inc. Wallingford CT 2009. (55)

Hehre, W. J. Self—Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian—Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. J. Chem. Phys. 1972, 56 (5), 2257.

(56)

Binkley, J. S.; Pople, J. a; Hehre, W. J. Self-Consistent Molecular Orbital Methods. 21. Small Split-Valence Basis Sets for First-Row Elements. J. Am. Chem. Soc. 1980, 102 (3), 939–947.

(57)

Gordon, M. S.; Binkley, J. S.; Pople, J. A.; Pietro, W. J.; Hehre, W. J. Self-Consistent Molecular-Orbital Methods. 22. Small Split-Valence Basis Sets for Second-Row Elements. J. Am. Chem. Soc. 1982, 104 (10), 2797–2803.

(58)

Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. Self-Consistent Molecular Orbital Methods. XXIII. A Polarization-Type Basis Set for Second-Row Elements. J. Chem. Phys. 1982, 77 (7), 3654.

(59)

Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. V. R. Efficient Diffuse Function-Augmented Basis Sets for Anion Calculations. III. The 3-21+G Basis Set for First-Row Elements, Li–F. J. Comput. Chem. 1983, 4 (3), 294–301.

(60)

Frisch, M. J.; Pople, J. a.; Binkley, J. S. Self-Consistent Molecular Orbital Methods 25.

ACS Paragon Plus Environment

55

Journal of Chemical Theory and Computation

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 56 of 73

Supplementary Functions for Gaussian Basis Sets. J. Chem. Phys. 1984, 80 (7), 3265. (61)

Raghavachari, K.; Trucks, G. W. Highly Correlated Systems. Excitation Energies of First Row Transition Metals Sc–Cu. J. Chem. Phys. 1989, 91 (2), 1062.

(62)

Rassolov, V. a.; Pople, J. a.; Ratner, M. a.; Windus, T. L. 6-31G∗ Basis Set for Atoms K through Zn. J. Chem. Phys. 1998, 109 (4), 1223.

(63)

Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. a. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117 (19), 5179–5197.

(64)

Wachters, a. J. H. Gaussian Basis Set for Molecular Wavefunctions Containing ThirdRow Atoms. J. Chem. Phys. 1970, 52 (3), 1033.

(65)

Hay, P. J. Gaussian Basis Sets for Molecular Calculations. The Representation of 3d Orbitals in Transition-Metal Atoms. J. Chem. Phys. 1977, 66 (10), 4377.

(66)

Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. a. Self-Consistent Molecular Orbital Methods. XX. A Basis Set for Correlated Wave Functions. J. Chem. Phys. 1980, 72 (1), 650.

(67)

McLean,

a. D.; Chandler, G. S. Contracted Gaussian Basis Sets for Molecular

Calculations. I. Second Row Atoms, Z=11–18. J. Chem. Phys. 1980, 72 (10), 5639. (68)

Li, P.; Roberts, B. P.; Chakravorty, D. K.; Merz, K. M. Rational Design of Particle Mesh Ewald Compatible Lennard-Jones Parameters for +2 Metal Cations in Explicit Solvent. J. Chem. Theory Comput. 2013, 9 (6), 2733–2748.

ACS Paragon Plus Environment

56

Page 57 of 73

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

Journal of Chemical Theory and Computation

(69)

Giammona, D. A. Ph.D. Thesis, University of California, Davis, 1984.

(70)

Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 (2), 926.

(71)

Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins Struct. Funct. Bioinforma. 2006, 65 (3), 712–725.

(72)

Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98 (12), 10089.

(73)

Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103 (19), 8577.

(74)

Adelman, S. A. Generalized Langevin Equation Approach for Atom/solid-Surface Scattering: General Formulation for Classical Scattering off Harmonic Solids. J. Chem. Phys. 1976, 64 (6), 2375.

(75)

Pastor, R. W.; Brooks, B. R.; Szabo, A. An Analysis of the Accuracy of Langevin and Molecular Dynamics Algorithms. Mol. Phys. 1988, 65 (6), 1409–1419.

(76)

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 (2), 227–249.

(77)

VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian

ACS Paragon Plus Environment

57

Journal of Chemical Theory and Computation

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 58 of 73

and Plane Waves Approach. Comput. Phys. Commun. 2005, 167 (2), 103–128. (78)

Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. cp2k: Atomistic Simulations of Condensed Matter Systems. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4 (1), 15–25.

(79)

Maseras, F.; Morokuma, K. IMOMM: A New Integratedab Initio + Molecular Mechanics Geometry Optimization Scheme of Equilibrium Structures and Transition States. J. Comput. Chem. 1995, 16 (9), 1170–1179.

(80)

Laino, T.; Mohamed, F.; Laio, A.; Parrinello, M. An Efficient Real Space Multigrid QM/MM Electrostatic Coupling. J. Chem. Theory Comput. 2005, 1 (6), 1176–1184.

(81)

Laino, T.; Mohamed, F.; Laio, A.; Parrinello, M. An Efficient Linear-Scaling Electrostatic Coupling for Treating Periodic Boundary Conditions in QM/MM Simulations. J. Chem. Theory Comput. 2006, 2 (5), 1370–1378.

(82)

Hutter, J.; Parrinello, M.; Lippert, G. The Gaussian and Augmented-Plane-Wave Density Functional Method for Ab Initio Molecular Dynamics Simulations. Theor. Chem. Acc. 1999, 103, 124–140.

(83)

Krack, M.; Parrinello, M. All-Electron Ab-Initio Molecular Dynamics. Phys. Chem. Chem. Phys. 2000, 2 (10), 2105–2112.

(84)

Spencer, J.; Alavi, A. Efficient Calculation of the Exact Exchange Energy in Periodic Systems Using a Truncated Coulomb Potential. Phys. Rev. B 2008, 77 (19), 193110.

(85)

Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132 (15), 154104.

ACS Paragon Plus Environment

58

Page 59 of 73

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

Journal of Chemical Theory and Computation

(86)

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.

(87)

Blöchl, P. E. Electrostatic Decoupling of Periodic Images of Plane-Wave-Expanded Densities and Derived Atomic Point Charges. J. Chem. Phys. 1995, 103 (17), 7422.

(88)

Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81 (1), 511.

(89)

Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31 (3), 1695–1697.

(90)

Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nosé–Hoover Chains: The Canonical Ensemble via Continuous Dynamics. J. Chem. Phys. 1992, 97 (4), 2635.

(91)

Carvalho, A. T. P.; Teixeira, A. F. S.; Ramos, M. J. Parameters for Molecular Dynamics Simulations of Iron-Sulfur Proteins. J. Comput. Chem. 2013, 34 (18), 1540–1548.

(92)

Dauber, P.; Hagler, A. T. Crystal Packing, Hydrogen Bonding, and the Effect of Crystal Forces on Molecular Conformation. Acc. Chem. Res. 1980, 13 (4), 105–112.

(93)

Thompson, H. P. G.; Day, G. M. Which Conformations Make Stable Crystal Structures? Mapping Crystalline Molecular Geometries to the Conformational Energy Landscape. Chem. Sci. 2014, 5 (8), 3173.

(94)

Ernzerhof, M.; Scuseria, G. E. Assessment of the Perdew–Burke–Ernzerhof ExchangeCorrelation Functional. J. Chem. Phys. 1999, 110 (11), 5029.

(95)

Adamo, C.; Barone, V. Toward Reliable Density Functional Methods without Adjustable

ACS Paragon Plus Environment

59

Journal of Chemical Theory and Computation

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 60 of 73

Parameters: The PBE0 Model. J. Chem. Phys. 1999, 110 (13), 6158. (96)

Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77 (18), 3865–3868.

(97)

Hubbard, J. Electron Correlations in Narrow Energy Bands. Proc. R. Soc. A Math. Phys. Eng. Sci. 1963, 276 (1365), 238–257.

(98)

Cococcioni, M.; de Gironcoli, S. A Linear Response Approach to the Calculation of the Effective Interaction Parameters in the LDA+U Method. Phys. Rev. B 2004, 71 (3), 35105.

(99)

Kulik, H. J.; Cococcioni, M.; Scherlis, D. A.; Marzari, N. Density Functional Theory in Transition-Metal Chemistry: A Self-Consistent Hubbard U Approach. Phys. Rev. Lett. 2006, 97 (10), 103001.

(100) Guidon, M.; Schiffmann, F.; Hutter, J.; Vandevondele, J. Ab Initio Molecular Dynamics Using Hybrid Density Functionals. J. Chem. Phys. 2008, 128 (2008). (101) Guidon, M.; Hutter, J.; VandeVondele, J. Robust Periodic Hartree-Fock Exchange for Large-Scale Simulations Using Gaussian Basis Sets. J. Chem. Theory Comput. 2009, 5 (Mic), 3010–3021. (102) Sit, P. H.-L.; Car, R.; Cohen, M. H.; Selloni, A. Simple, Unambiguous Theoretical Approach to Oxidation State Determination via First-Principles Calculations. Inorg. Chem. 2011, 50 (20), 10259–10267. (103) Attia, A. A. A. A.; Cioloboc, D.; Lupan, A.; Silaghi-Dumitrescu, R. Fe-O versus O-O Bond Cleavage in Reactive Iron Peroxide Intermediates of Superoxide Reductase. J. Biol.

ACS Paragon Plus Environment

60

Page 61 of 73

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

Journal of Chemical Theory and Computation

Inorg. Chem. 2013, 18 (1), 95–101. (104) Keith, T. A. AIMAll (15.09.27). TK Gristmill Software, Overland Park KS 2015. (105) Cremer, D.; Kraka, E. Chemical Bonds without Bonding Electron Density- Does the Difference Electron-Density Analysis Suffice for a Description of the Chemical Bond? Angew. Chemie Int. Ed. 1984, 41 (8), 627–628. (106) Horch, M.; Pinto, a. F.; Utesch, T.; Mroginski, M. a.; Romão, C. V.; Teixeira, M.; Hildebrandt, P.; Zebger, I. Reductive Activation and Structural Rearrangement in Superoxide Reductase: A Combined Infrared Spectroscopic and Computational Study. Phys. Chem. Chem. Phys. 2014, 16 (27), 14220–14230.

ACS Paragon Plus Environment

61

Journal of Chemical Theory and Computation

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 62 of 73

Insert Table of Contents artwork here

ACS Paragon Plus Environment

62

Page 63 of 73

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

Journal of Chemical Theory and Computation

84x140mm (300 x 300 DPI)

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

84x139mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 64 of 73

Page 65 of 73

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

Journal of Chemical Theory and Computation

177x103mm (300 x 300 DPI)

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

177x81mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 66 of 73

Page 67 of 73

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

Journal of Chemical Theory and Computation

177x58mm (300 x 300 DPI)

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

177x101mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 68 of 73

Page 69 of 73

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

Journal of Chemical Theory and Computation

177x130mm (300 x 300 DPI)

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

177x101mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 70 of 73

Page 71 of 73

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

Journal of Chemical Theory and Computation

177x122mm (300 x 300 DPI)

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

197x231mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 72 of 73

Page 73 of 73

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

Journal of Chemical Theory and Computation

77x44mm (300 x 300 DPI)

ACS Paragon Plus Environment