Unraveling the Addition–Elimination Mechanism of EPSP Synthase

Aug 22, 2017 - Institute of Biological Sciences, Federal University of Pará, Belém, PA 66075-110, Brazil. ‡ Institute of Exact and Natural Science...
0 downloads 0 Views 2MB Size
Subscriber access provided by UNIV OF ARIZONA

Article

Unraveling the Addition-Elimination Mechanism of EPSP Synthase through Computer Modeling Alberto Monteiro dos Santos, Anderson Henrique Lima, Claudio Nahum Alves, and Jeronimo Lameira J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b05063 • Publication Date (Web): 22 Aug 2017 Downloaded from http://pubs.acs.org on August 25, 2017

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

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

Page 1 of 34

The Journal of Physical Chemistry

1 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

Unraveling the Addition-Elimination Mechanism of EPSP Synthase through Computer Modeling

Alberto M. dos Santos1, Anderson H. Lima1, Cláudio Nahum Alves,2 Jerônimo Lameira1* 1

Institute of Biological Sciences, Federal University of Pará, 66075-110, Belém, PA, Brazil.

2

Institute of Exact and Natural Sciences, Federal University of Pará, 66075-110, Belém, PA, Brazil.

*Corresponding author: Jerônimo Lameira; Phone: +55 91-32018905; email: [email protected]

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 2 of 34

2 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

ABSTRACT

Enolpyruvyl transfer from phosphoenolpyruvate (PEP) to the hydroxyl group of shikimate-5-OH-3phosphate (S3P) is catalyzed by 5-enolpyruvylshikimate 3-phosphate (EPSP) synthase in a reaction that involves breaking the C-O bond of PEP. Catalysis involves an addition–elimination mechanism with the formation of a tetrahedral intermediate (THI). Experiments have elucidated the mechanism of THI formation and breakdown. However, the catalytic action of EPSP synthase and the individual roles of catalytic residues Asp313 and Glu341 remains unclear. We have used a hybrid quantum mechanical/molecular mechanical (QM/MM) approach to explore the free energy surface in a reaction catalyzed by EPSP synthase. The Glu341 was the most favorable acid/base catalyst. Our results indicate that the protonation of PEP C3 precedes the nucleophilic attack on PEP C2 in the addition mechanism. Also, the breaking of the C-O bond of THI to form an EPSP cation intermediate must occur before proton transfer from PEP C3 to Glu341 in the elimination mechanism. Analysis of the FES supports cationic intermediate formation during the reaction catalyzed by EPSP synthase. Finally, the computational model indicates a proton transfer shift (Hammond shift) from Glu341 to C3 for an enzyme-based reaction with the shifted transition state, earlier than in the reference reaction in water.

ACS Paragon Plus Environment

Page 3 of 34

The Journal of Physical Chemistry

3 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

INTRODUCTION

The shikimate pathway is essential to the synthesis of aromatic amino acids and almost all other aromatic compounds in algae, vascular plants, bacteria, fungi, and parasites of the phylum Apicomplexa.1 The enzymes of the shikimate pathway are absent in mammals making them an attractive target for the development of new antimicrobial agents.2,3 The 5-enolpyruvylshikimate 3phosphate (EPSP) synthase, also known as AroA, participates in the biosynthesis of the aromatic amino acids via the shikimate pathway. AroA forms EPSP and inorganic phosphate (Pi) from shikimate 5-OH-3-phosphate (S3P) and phosphoenolpyruvate (PEP).4 The reaction catalyzed by EPSP synthase occurs in an addition–elimination mechanism (Scheme 1). Addition consists of the PEP activation by C3 protonation associated with a nucleophilic attack on carbon 2 (C2). The 5-hydroxyl group of S3P and its covalent link to C2 of PEP forms a tetrahedral intermediate (THI). According to Scheme 1, the elimination step consists of THI C3 deprotonation and the C-O bond cleavage that leads to Pi departure.5 EPSP synthase accelerates THI breakdown by a factor of >105, relative to nonenzymatic breakdown.6 The biochemical reaction catalyzed by this enzyme is unusual, since the enzymatic synthesis of EPSP from PEP and S3P involves C-O cleavage of the pyruvate ester.7 Normally, enzymes that utilize phosphate as a substrate break the high-energy P-O bond.

Scheme 1. Reaction catalyzed by EPSP synthase and structures of shikimate 5-OH-3-phosphate (S3P), phosphoenolpyruvate (PEP), tetrahedral intermediate (THI), and 5-enolpyruvyl shikimate-3-phosphate (EPSP).

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 4 of 34

4 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

EPSP synthase is a member of the enolpyruvyl transferase family of enzymes that has only one homolog, that is, UDP-N-acetyl-glucosamine enolpyruvyl transferase (MurA), an enzyme that catalyzes the first committed step in peptidoglycan synthesis. The enzymes involved in peptidoglycan biosynthesis are among the best-known targets in the search for new antibiotics.8,9 AroA and MurA have been the targets of antibiotic development for four decades, but inhibitor design remains challenging.10 The shikimate pathway also presents an attractive target for malaria chemotherapy because shikimate analogs have been shown to inhibit the growth of Plasmodium falciparum.11 Some potent inhibitors of AroA have been created based on THI structure.12 On the other hand, analogs of the transition state can bind millions of times tighter than substrates and show promise for drug development for several targets.13 In light of the biological importance of EPSP synthase, a considerable body of experimental data aimed at understanding the details of its catalyzed reaction has been accumulated over the years.1–13 However, despite these data, there is no general agreement on the molecular mechanisms for AroA-catalyzed EPSP hydrolysis. THI partitioning experiments identified Lys22 and Glu341 residues as acid/base catalytic residues14 in the proposed mechanism (Scheme 1). Other experiments showed that catalytic residues that promote THI breakdown are the same in both forward and reverse reactions, where both Asp313 and Glu341 can deprotonate the C3 of THI.15 Three possible AroA mechanisms can be considered.16 The reaction could begin with (i) C3 protonation in a stepwise AE+AN mechanism or (ii) concerted C3 protonation and nucleophilic attack at C2 (ANAE) or (iii) non-activated nucleophilic attack at C2 stepwise AN+AE mechanism.10 Scheme 2 shows three possible mechanisms for a reaction catalyzed by EPSP synthase.

ACS Paragon Plus Environment

Page 5 of 34

The Journal of Physical Chemistry

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

Scheme 2. EPSP synthase catalysis proposals. (a) Stepwise AN+AE mechanism (b) concerted C3 protonation and nucleophilic attack on C2 (ANAE) mechanism. (c) Non-activated nucleophilic attack at C2 stepwise AN +AE mechanism. The point at which S3P corresponds to shikimate-3-phosphate and X could be an Asp or Glu residue.

The modeling of complex molecular systems is crucial to advancing our understanding of biological systems and the rational design of new drugs.17 Recently, we have used a free energy surface and quantum mechanics/molecular mechanics (QM/MM) approach18 to investigate biomolecular systems with an emphasis on enzymatic reactions.19–21 In the QM/MM approach, a small part of the system (ligand/substrate species) is described by quantum mechanics, while the protein and solvent environment is represented by MM force fields.18 Knowledge of the molecular mechanism leads to the definition of transition states and other species (intermediates) involved along the reaction path that provides a step toward the design and synthesis of new inhibitors that could bind tightly to AroA. Herein, we report on the first QM/MM free energy surface study of the addition and elimination mechanism of EPSP synthase. In the present study, we examined the accuracy of the AM1/d-PhoT/MM and M06-2X/6-31G(d)/MM methods for the calculation of the activation free energy of the reaction

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 6 of 34

6 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

catalyzed by this enzyme. Finally, we also decipher the specific roles of the Asp313 and Glu341 catalytic residues and show that the AroA attains a Hammond shift relative to the reference reaction in solution.

METHODS

The system In the present study, the computational model for the QM/MM MD calculations was based on the crystal structure of the THI-EPSP synthase complex, with PDB code 1Q36,22 for which we performed an Ala313Asp mutation that reversed the system to the wild-type enzyme (EPSPS wildtype). To this end, AM1/d-PhoT Hamiltonian23 was employed to describe the QM part while the other atoms of the system were described using the CHARMM32/TIP3P24 force fields. Before starting the molecular dynamic simulations, an assignment of the protonation states of the residues at pH = 7 was carried using the empirical propKa program25, excluding those residues directly involved in the mechanism. After adding the hydrogen atoms to the full structure, a series of optimization algorithms26 was applied. All the heavy atoms of the protein and substrate were restrained by means of a Cartesian harmonic umbrella with a force constant of 239 kcal · mol-1 ·Å-2 (1000 kJ mol-1 ·Å-2). Subsequently, the system was fully relaxed, but the peptide backbone was restrained with a lower constant of 100 kJ · mol-1 ·Å-2. Then, the optimized protein was placed in an octahedral box of pre-equilibrated water (80 Å side), using the center of mass of the complex as the geometrical center. Then, 100 ps of hybrid QM/MM Langevin-Verlet MD at 300 K and in a canonical thermodynamic ensemble (NVT) were used to equilibrate the system. For the hybrid QM/MM calculations, the atoms of the substrate and side chain of the Lys22, Lys411, Asp313 and Glu341 residues were selected for treatment by QM as described in the Supporting Information (Figure S1), using a semi-empirical AM1/d-PhoT Hamiltonian.23 The other atoms of the system, protein, and water molecules, were described using the ACS Paragon Plus Environment

Page 7 of 34

The Journal of Physical Chemistry

7 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

CHARMM32/TIP3P24 force fields, respectively. The number of QM atoms then proved to be 39, while the final system consisted of 58,012 atoms. A recent study by of Martínez et al.27 proposed the use of a charge shift analysis to explore the minimum number of protein residues needed to attain a quantitative agreement with large-QM simulations. Based on geometry optimizations and transition-state searches, the authors have proposed that approximately 200–300 atoms are needed to fully describe the mechanism catalyzed by catechol O-methyltransferase. They concluded that the properties of enzyme catalysis as obtained by simulations depend on the size of the QM region in the QM/MM calculations.27 However, using geometry optimizations, Liao and Thiel28 proposed that the size of the QM region has little effect on the single-point QM/MM relative energies (typically 1–2 kcal.mol-1), and even when using a large QM region (408 atoms) uncertainties arise from the choice of the model. In addition, Liao and Thiel28 recommended that, in QM/MM studies of enzyme reaction mechanisms, a relatively small QM region should first be used to establish the reaction mechanism. Jindal and Warshel29 indicated that the size dependence is not as crucial as one might suspect. This recent study points out that the residues that are in direct contact with the reacting atoms can be involved in a major charge transfer or even in another reaction path. However, the QM description of residues that are not in direct contact do not result in a major change to the activation barriers to chemical transformation.29 Therefore, we are assuming that the QM region selected in the present work can fully describe the effect of the protein environment on the mechanism being studied. Due to the number of degrees of freedom, any residue that was 20 Å from any of the atoms of the initial inhibitor was kept frozen in the subsequent calculations. Cut-offs for the non-bonding interactions were applied using a switching scheme, within a range radius of 14.5 to 16 Å. Next, the system was equilibrated by QM/MM MD at a temperature of 300 K for 1.3 ns. The computed RMSD for the protein during the last 500 ps is always below 0.9 Å. Furthermore, the RMS of the temperature along the different equilibration steps was always less than 2.5 K and the variation coefficient of the potential energy during the dynamics simulations was never more than 0.3%. ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 8 of 34

8 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

Free Energy Surface Increased computer power, new parallel computing approaches and advances in electronic structure calculations have resulted in significant progress in ab initio QM/MM simulation methodology.30 These developments have possibly given rise to new means of overcoming the challenges facing ab initio calculations of free energy profiles in solutions and proteins.31–34 However, it remains extremely difficult to obtain computationally converging sampling for higher-level QM/MM (QM(HL)/MM) energy corrections for free energy surfaces in condensed phases (due the large number of gradient values evaluation). In the present case, we described the quantum region using AM1/dPhoT Hamiltonian23 to compute the activation free energies at a significantly reduced computational cost while incorporating a d-extension for the phosphorus atom and modified AM1 parameters for the oxygen and hydrogen atoms. The AM1/d-PhoT potential has been designed to accurately reproduce high-level DFT results such as geometries, dipole moments, proton affinities, and relative energies for a large set of molecules, complexes, and chemical reactions relevant to biological phosphoryl transfer.23 The results obtained by Nam et al. suggest that the AM1/d-PhoT Hamiltonian provides significantly higher accuracy with results that are in good agreement with the density-functional calculations at the B3LYP level with a reduction in the computational cost of 3–4 orders of magnitude.23 To study the addition mechanism, we constructed 2D-PMF using a combination of the R1–R2 antisymmetric distance and R3 distance. R1–R2 corresponds to the proton transfer from Glu341 to C3, while R3 corresponds to the nucleophilic attack of S3P- 5'O on the activated PEP C2. For elimination, we constructed 2D-PMF using a combination of the R4 distance and R5–R6 antisymmetric distance. R4 corresponds to the phosphate release from THI, while R5–R6 corresponds to the proton transfer from THI C3 to Glu341 (see Scheme 3).

ACS Paragon Plus Environment

Page 9 of 34

The Journal of Physical Chemistry

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

Scheme 3. Proposed mechanism for the addition and elimination steps of the reaction catalyzed by EPSP synthase. *Key distances employed to explore the molecular mechanism are depicted in blue.

We also carried out a 1D-PMF simulation for proton transfer from Lys22 to the phosphate bridging oxygen (Figure S2). All the PMFs were calculated using the weighted histogram analysis method (WHAM) combined with an umbrella sampling approach.35,36 All the simulations were carried out using the pDynamo program.37 The PMF calculation requires a series of molecular dynamics simulations in which the distinguished reaction coordinate variable, ξ, is constrained around particular values.36 The values of the variables sampled during the simulations are then pieced together to construct a distribution function from which the PMF is obtained as a function of the distinguished reaction coordinate (W(ξ)). The PMF is related to the normalized probability of finding the system at a particular value of the chosen coordinate, given by equation (1):38

( )(

)

W (ξ ) = C − kT ln ∫ ρ r N δ ξ (r N ) − ξ dr N (1) The activation free energy can be then expressed as:34

[

∆G ‡ (ξ ) = W (ξ ‡ ) − W (ξ R ) + Gξ (ξ R )

]

(2)

where thesuperscripts indicate the value of the reaction coordinate at the reactants (R) and TS ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 10 of 34

10 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

and G ξ (ξ R ) indicate the free energy associated with setting the reaction coordinate to a specific value at the reactant state. Normally this last term makes only a small contribution and the activation free energy is directly estimated from the PMF change between the maximum of the profile and the reactant’s minimum:

∆G‡ (ξ ) ≈ W (ξ ‡ ) − W (ξ R )= ∆W ‡ (ξ )

(3)

A total of 41 simulations were performed with different values of the R1–R2 antisymmetric combination of distances (-2.0 Å–2.0 Å), with an umbrella force constant of 119.5 (500 kJ mol-1 ·Å-2) kcal.mol-1 applied to this reaction coordinate. In addition, 41 simulations were performed at different values of R3 (1.0 Å–3.0 Å), again with an umbrella force constant of 119.5 kcal.mol-1 for this antisymmetric combination of distances. Consequently, 1681 simulation windows were needed to obtain the 2D-PMF value for the addition mechanism. The values of the variables sampled during the simulations were then pieced together to construct a full distribution function from which the 2D-PMF was obtained. In each window, a 10-ps relaxation time was followed by 10 ps of production with a time step of 0.5 fs, given the nature of the chemical step involving hydrogen transfer. The Verlet algorithm was used to update the velocities. The same protocol as that used in the simulation was used to describe the 2D-PMF of the elimination mechanism, using the R4 and R5–R6 distances, in which 1681 simulation windows were again needed to obtain the value of 2D-PMF. Three one-dimensional PMFs were also computed using separate 1-D reaction coordinates: the coordinates of the antisymmetric combination of distances R1–R2, R3, R4, and R5–R6. A total of 41 simulations were performed using different values of R1–R2 (-2.0 Å–2.0 Å), R3 (1.0 Å–3.0 Å), R4 (1.0 Å–3.0 Å), and R5–R6 (-2.0 Å–2.0 Å) with an umbrella force constant of 119.5 kcal.mol-1. In each window the protocol of the simulation was as described above for 2D-PMF.

ACS Paragon Plus Environment

Page 11 of 34

The Journal of Physical Chemistry

11 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

Activation free energy at M06-2X/6-31G(d)/MM level The paradynamics (PD) model is a recent advance that enables the incorporation of QM/MM into the free energy calculation.39–41 The PD model is based on a number of earlier works, including those that have established means of efficiently sampling the reference potential (RP) using the Empirical Valence Bond42 (EVB) method, as well as ideas for refining EVB RP43 to minimize the energy difference between the EVB RP and the QM/MM target potential (TP). In the present study, we have used AM1/d-PhoT/MM as RP and M06-2X/6-31G(d)/MM as TP to further explore the challenge of evaluating higher-level energy corrections for the QM/MM activation free energies of a reaction catalyzed by EPSP synthase. M06-2X is shown to describe noncovalent interactions (important to biomolecules such as proteins) better than the density functionals which are currently in common use.44,45 Moreover, M06-2X has been successfully used to study the energetics of the hydrolysis of phosphate esters46 and the catalytic effect of enzymes.47 Basically, we have obtained a single-point at a M06-2X/6-31G(d)/MM potential using representative snapshots from RS, INT, and TS ensembles from 2D-PMF QM/MM simulations computed at the AM1/d-PhoT potential, for which a total of 2000 representative snapshots for each state were used in the simulations. In this way, the extensive sampling required to calculate the QM(HL)/MM free energy barrier can be done on a computationally inexpensive reference potential rather than directly on the expensive target potential.48 The free energy required to move from AM1/d-PhoT/MM to M06-2X/6-31G(d)/MM at the RS position was determined using a straightforward one-step free energy perturbation: ∆∆/ →  = −  〈 −!" − "/ #$〉/ &

(4)

The same idea was used to compute the free energy required to move from AM1/d-PhoT/MM to M06-2XG(d)/MM at the TS position that was obtained using: ∆∆/ → ' = −  〈 −!" − "/ #$〉/ &

(5)

It is important to emphasize that these equations are based on an FEP approach and reference

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 12 of 34

12 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

potential concept, which is used in the PD approach. On the other hand, we are not interested in refining RP as is normally done in PD, but our goal is to compute the energy gap between RP and TP instead. Then, we use it to compute the activation energy at M06-2X/6-31G(d)/MM level through: ‡ ‡ = ∆∆(/ + ∆∆/ → + − ∆∆/ →-+ ∆∆(

(6)

‡ where the activation free energy computed at the AM1/d-PhoT/MM level (∆∆(/ 

was obtained from the PMF simulations, as explained above.

RESULTS AND DISCUSSION

EPSP Synthase Addition Step Despite a consensus based on many experimental data,4,5,10,15 a clear statement about the molecular details of a reaction mechanism catalyzed by EPSP synthase has yet to be determined. Computational descriptions of chemical reactions in enzymatic environments are usually based on the selection of distinguished reaction coordinates.47,49 Therefore, in the present study we used free energy surface simulations to explore and test the catalytic mechanism of EPSP synthase. This section presents results for FES simulations for the addition mechanism of reaction catalyzed by EPSP synthase. We have traced a 2D-PMF using the combination R1–R2 antisymmetric distance and R3 distance (see Scheme 3), where R1–R2 corresponds to the proton transfer (PT) from O of the carbonyl of Glu341 to PEP C3, while R3 corresponds to the nucleophilic attack of S3P- 5'O on C2 of the activated PEP. The addition step in the EPSP synthase requires the protonation of C3 of PEP and the deprotonation of S3P O5’H.14 Thus, the addition mechanism should start with deprotonated S3P (S3P-, Scheme 3), where S3P- is likely to form a stable salt bridge with Lys22. The reaction free energy for this step can be estimated by: ∆ = 2.3'12 3445 − 12 677845

ACS Paragon Plus Environment

(7)

Page 13 of 34

The Journal of Physical Chemistry

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

The above equation and pKa values of the side chain of Lys (pKa = 10.5) and the O5’H group of S3P (pKa = 15.5) can be used to evaluate the free energy associated with proton transfer from the O5’H group of S3P to lysine amine as being close to 6.8 kcal.mol-1 in aqueous solution. The reaction free energy obtained from PMF corresponds to 4.9 kcal.mol-1 (see Figure S2) in solution, which is in reasonable agreement with the experimental estimates using pKa values. In this case, we observed an exothermic process in the proteins and an endothermic process in a solution (see Figure S2). This option is supported by partitioning experiments14 and by our PMF results for the THI-EPSP synthase complex that suggests that Lys22 is possibly protonated (see Figure S3). Figure 1a shows the 2D-PMF obtained for the addition mechanism of in the EPSP enzyme at the AM1/d-PhoT/MM level as a function of a combination of the distances (R1–R2, and R3, see Scheme 3) describing the breaking and forming of bonds. In this PMF, the RS structure corresponds to S3P- and PEP. In the TS1 (the transition state for PT), the bond distance for R1 corresponds to 2.2 Å, whereas R2 corresponds to 1.8 Å, indicating a very advanced stage for the proton transfer, while the bond distance for R3 remains constant at around 2.8 Å during PT (Table 1). On this FES (Figure 1a), INT1 corresponds to the first cationic intermediate, a structure where a nucleophilic attack has not yet occurred. In the TS2 structure (in the transition state of a nucleophilic attack) the bond distance for R3 correspond to 2.3 Å, indicating C-O formation (see Table 1 and Figure 2).

Figure 1. 2D-PMF (in kcal.mol-1) calculated at AM1/d-PhoT/MM level (a) R1–R2 and R3 antisymmetric combinations

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 14 of 34

14 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

of distances describing addition mechanism (b) R4 and R5–R6 antisymmetric combinations of distances describing elimination mechanism. Values of contour potential energy lines are reported in kcal.mol-1 and coordinates in Å. Values on the isoenergetic lines are reported for each 1 kcal.mol-1

Table 1. Average Key Distances (Å) for structures obtained from 2D-FESsa Addition Distance

RS

TS1

INT1

TS2

THI

R1: (Glu341)O-H(Glu341)

1.04(0.03)

1.58(0.05)

2.82(0.06)

2.78(0.04)

2.32(0.11)

R2: C3(PEP)-H(Glu341)

2.65(0.04)

1.47(0.04)

1.20(0.05)

1.17(0.03)

1.10(0.03)

R3: O5’(S3P)-C2(PEP)

2.82(0.04)

2.86(0.04)

2.76(0.03)

2.29(0.04)

1.44(0.02)

R4: C2(PEP)-O14(PEP)

1.43(0.03)

1.44(0.03)

1.34(0.02)

1.34(0.02)

1.46(0.02)

Elimination Distance

a

THI*

TS3

INT2

TS4

PS

R5: (Glu341)O-H(THI)

1.12(0.05)

1.10(0.02)

1.19(0.03)

1.39(0.05)

2.14(0.05)

R6: C3(THI)-H(THI)

4.44(0.15)

2.39(0.16)

2.21(0.04)

1.29(0.05 )

1.09(0.04)

R3: O5’(THI)-C2(THI)

1.44(0.02)

1.42(0.02)

1.34(0.02)

1.37(0.02)

1.40(0.02)

R4: C2(THI)-O14(THI)

1.41(0.02)

1.54(0.02)

2.76(0.05)

2.77(0.05)

2.76(0.05)

Reaction coordinates R1−R6 defined in Scheme 3.

* Protonated tetrahedral intermediate

The resulting FES path has an exothermic reaction with low activation barriers thus confirming the viability of a stepwise mechanism that leads to tetrahedral intermediate (THI) formation. Therefore, simulations reveal that the addition step proceeds via a stepwise mechanism, starting with the transfer of protons from Glu341 to C3 of PEP, for which the activation free energy at AM1/d-PhoT level is 6.1 kcal.mol-1. Then, the nucleophilic attack at PEP C2 takes place with an activation barrier of 4.7 kcal.mol-1 (see Table 2). In accordance with the obtained free energy surface (Figure 1a) the PT from

ACS Paragon Plus Environment

Page 15 of 34

The Journal of Physical Chemistry

15 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

Glu341 to PEP C3 is a rate-limiting step for the addition mechanism.

Figure 2. Representative structures along the 2D-PMF obtained from the QM/MM free energy simulations for addition step. Structures are named according to the locations along the 2D PMF shown in Figure 1A for the addition step simulation. The steps for elimination are as follows: (R) to (TS1) corresponds to PT from Glu341 to C3 of PEP; (TS1) to (INT1) corresponds to the formation of a PEP cationic transition state; (INT1) to (TS2) corresponds to the nucleophilic attack at C2. The initial deprotonation of S3P is not shown.

Table 2. Activation free energy calculations for EPSP synthase at AM1/d-PhoT/MM and M06-2X/6-31G(d)/MM levels (in kcal.mol-1)

Addition C3

Elimination

Nucleophil

C-O bond

Method

a

C3 Deprotonation Protonation

ic Attack

cleavage

AM1/d 2D PMF

6.1

4.7

9.3

7.9

M06-2X/6-31G(d)

16.3a

11.3a

18.6a

19.2a

The values on the table were obtained considering the perturbation from AM1/d-PhoT/MM reference potential to M06-

2X/6-31G(d)/MM target potential.

As mentioned in the method section, we used AM1/d-PhoT/MM for RP and M06-2X/6ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 16 of 34

16 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

31G(d)/MM for TP to determine the higher-level energy corrections for QM/MM activation free energies of a reaction catalyzed by EPSP synthase. Representative snapshots from the RS, TS1, INT1, and TS2 ensembles from AM1/d-PhoT 2D-PMF were used as the start point to evaluate the activation free energy of EPSP synthase at the M06-2X/6-31G(d)/MM level by following the approach described in the method section. Note that estimating the average difference between TP and RP allows us to compute the activation free energies at a significantly reduced computational cost, using only two trajectories: one at the reactants and one at the transition state. Therefore, the reference potential concept appears to offer a feasible approach in terms of computational cost for calculating the M062X/6-31G(d)/MM free-energy barrier for an enzymatic reaction.40,41 At this point, it is important to clarify that we used paradynamics to correct the AM1d-PhoT energies, in which we computed the energy gap for 2000 configurations of RS and TS. Here, we assume that AM1/d reproduces the DFT geometry accurately and, consequently, we should be matching the correct stationary points in both potential AM1/d and DFT. Corrections using SCS-MP2 would incur a considerable computational time and we would not be able to guarantee that we could match the stationary points. After M06-2X/6-31G(d)/MM corrections, the resulting activation barrier for the transfer of protons from Glu341 to PEP C3 was 16.3 kcal.mol-1, which is in good agreement with the experimental estimates15 of >15 kcal.mol-1 and with the activation barrier (14.8 kcal.mol-1) calculated at the M062X/6-311++G(d,p)/MM level for this step. Therefore, we used M06-2X/6-31G(d)/MM, which is computationally less demanding when studying the other steps of the reaction. Applying M06-2X corrections, the barrier for nucleophilic attack is 11.3 kcal.mol-1. Both the M06-2X/6-31G(d) and semiempirical barriers indicate that the transfer of protons from Glu341 to C3 of PEP is the rate-limiting step of the EPSP synthase addition mechanism. Thus, within the information on the full AM1/d-PhoT free energy surface and M06-2X/6-31G(d)/MM energies of this mechanism, we can suggest that the rate-limiting step of the addition mechanism is that controlled by the TS1, which corresponds to the PT from Glu341 to C3 of PEP. Therefore, the protonation of C3 is the greater catalytic imperative in THI ACS Paragon Plus Environment

Page 17 of 34

The Journal of Physical Chemistry

17 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

formation which is in very good agreement with the proposal based on experimental data.15

EPSP Synthase Elimination Step In the present study, we have also explored the elimination of the reaction catalyzed by EPSP synthase. Before presenting the results for this step of the reaction, we should point out that the X-ray structure (PDB code 1Q36) of the THI-EPSP synthase complex suggests that there is a strong hydrogen bond interaction between the Nζ atom of Lys22 and the ester oxygen of THI. It is important to remember that the protein changes the environment of the proton donor and acceptor and consequently the apparent pKa of the residues are modified by the local environment. Here, we did not directly evaluate the free energy of moving the proton in every case, since in most cases we were interested in finding the protonation state that leads to the most probable mechanism (lower activation energy). In this sense, we presented the results of an exhaustive hybrid QM/MM exploration of the catalytic mechanism used by EPSP synthase (PMF for the phosphate release step), in which we started our study by exploring different possibilities of protonation (a total of 24 simulations, see Table S1 and Figure S3 in Supporting Information) for residues in the active site of EPSP including two key residues, namely, Glu341 and Asp313, where the combination Lys22/Lys411 protonated and Asp313/Glu341 deprotonated exhibited a lower activation barrier for phosphate release. For this reason, we obtained a pKa value of 3.5 for the side chain of Lys22 using the empirical propka25 for the THI-EPSP synthase complex. The suppression of the pKa value for the side chain of Lys22 indicates that a neutral form of Lys22 is possible since the ionizable groups can be shifted by local protein environments. It is important to clarify that the Lys22 forms salt bridge with S3P- before the nucleophilic attack of S3P- 5'O on the activated PEP C2 (addition step) and that the negative charge of the phosphate group is stabilized by another lysine residue (Lys44) in both the addition and elimination steps. Figure S5 describes the proton transfer from Lys22 to the bridging oxygen of THI

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 18 of 34

18 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

and the phosphate release calculated at the AM1/d-PhoT/MM level. The result suggests that the proton transfer may occur before the C-O bond breakdown of THI. Therefore, the protonated tetrahedral intermediate was used as a starting point for the simulations. It is important to note that the protonation of the THI phosphate allows a barrier that is more energetically favorable for the C-O bond cleavage of THI. To explore the elimination step of a reaction catalyzed by the studied enzyme, we traced a 2DFES using the R4 distance and R5–R6 antisymmetric distance combination, where R4 corresponds to the C-O bond breakdown of THI, while R5–R6 corresponds to the proton transfer from C3 of the PEP moiety to Glu341 (see Scheme 3). Furthermore, the FES results suggest that the activation barrier for proton transfer from C3 to Asp313 is greater than that involving proton transfer from PEP C3 to Glu341 (see Figure S4) which is in agreement with experimental study that identified Glu341 as a general acid/base catalytic residue.14 Figure 1b shows a 2D free energy surface (2D-FES) for the elimination step using a combination of distances R4 and R5–R6, as described above. The reaction starts with the protonated tetrahedral intermediate (THI*) and first leads to the formation of a cationic intermediate named INT2. The product of the reaction is EPSP on the 2D-FES (Figure 1b), thus completing the reaction cycle. The TS4 observed on the FES corresponds to the proton transfer from PEP C3 to Glu341, where the bond distance for R5 and R6 is 1.4 and 1.3 Å, respectively (see Table 1 and Figure 3), and R4 is around 2.8 Å, indicating the advanced progress of the breaking of the C-O bond such that, consequently, the phosphate is released from THI. This poses a new view of this mechanism in that elimination appears to occur not continuously but stepwise with no antiperiplanar requirement.

ACS Paragon Plus Environment

Page 19 of 34

The Journal of Physical Chemistry

19 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

Figure 3. Representative structures along the 2D-PMF obtained from the QM/MM free energy simulations for elimination. Structures named according to the locations along the 2D PMF shown in Figure 1B for the elimination step simulation. The steps for elimination are as follows: (THI*) to (TS3) corresponds to C-O bond cleavage followed by phosphate departure; (TS3) to (INT2) corresponds to the formation of an EPSP cationic transition state; (TS4) to (P) corresponds to C3 deprotonation from the C3 methyl group to Glu341 carboxyl group. The mechanism associated with (THI*) to (PS) is stepwise.

The computed barrier is about 9.3 kcal.mol-1 for phosphate release and 7.9 kcal.mol-1 for proton transfer (see Table 2). Clark and Berti4 reported experimental values for EPSP hydrolysis in solution with rate constants of 7.7 × 10-8 s-1. As it has been proposed that AroA-catalyzed EPSP hydrolysis is a factor of (1–5) × 103 faster than the non-enzymatic rate at pH 7, corresponding to 4–5 kcal.mol-1 of transition state stabilization, we estimated that ∆ ‡ = 18 kcal.mol-1 for EPSP hydrolysis. The M062X/6-31G(d)/MM calculated activation free energy for the C-O bond cleavage is 19.2 kcal.mol-1, which is in excellent agreement with the activation free energy value derived from the experimental estimation based on the rate constant, suggesting that the barrier estimated by the M06-2X/631G(d)/MM TP approach is useful for studying this mechanism since the predicted activation free

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 20 of 34

20 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

energy is in very good agreement with the catalytic activation observed experimentally. Obviously, the proposed mechanism requires further theoretical and experimental study. At present, we cannot exclude the possibility of the catalytic reaction used by EPSP synthase proceeding via another pathway.

Catalytic effect of EPSP To explore the origin of the catalytic power of EPSP synthase, it is important to have a clear idea of the activation barrier for the reaction in both enzymatic and aqueous solution environments. In this sense, we have determined the AM1/d-PhoT/MM 1D-FES for the proton transfer from Glu341 to C3 (PT in addition mechanism) and from C3 to Glu341 (PT in elimination mechanism) in both water and protein environments (see Figure 4). The calculated difference in the barrier height for the proton transfers in the two environments of 25 kcal.mol-1 at the AM1/d-PhoT level is significantly overestimated relative to the experimental value of 14 kcal.mol-1.50 The origin of the deviation between the catalytic effect described by the experimental data and the PMF could be attributed to the size of the QM region. For the reaction in the protein we have included the atoms of the substrate and the side chains of Lys22, Lys411, Asp313, and Glu341 in the QM region (see SI) while for reaction in water we have only the atoms of the substrate and side chain of Glu341. Other source of the deviation could be attributed to limitations of the AM1/d-PhoT method. On the other hand, these results also suggest that AM1/d-PhoT is useful for qualitatively describing the free energy barrier in both the enzymatic and solution reactions.

ACS Paragon Plus Environment

Page 21 of 34

The Journal of Physical Chemistry

21 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

Figure 4. Reaction profiles showing the addition and elimination of the EPSP synthase mechanism in enzyme (red line) and water (blue line) environments. A) Free energies for PT in addition mechanism B) Free energies for PT in elimination mechanism, where INT2 was used as the start point of the simulation.

At this point it is useful to address how the enzyme stabilizes in essentially the same transition state as a non-enzymatic reaction, thus lowering the activation energy. In accordance with the 1D-FES (Figure 4) results and experimental data,10,50 the proton transfer (PT) is significantly higher in water than in protein, given that an earlier TS in the protein is observed. As expected, our results also indicate a shift in the proton transfer from C3 of PEP to Glu341 (and vice versa) with the earlier TS in the protein. This observation agrees with the previous computational work of Lou et al., who proposed that a significant part of AroA’s catalytic strategy is to stabilize the positive charge in the EPSP cation by shifting the transition state prior to the acid-catalyzed reaction, where AroA accomplishes a Hammond shift with a transition state earlier in EPSP synthase than in an acid-catalyzed reaction in solution.10 Enzymes involved in the shikimate pathway have been investigated using computational methods, such as shikimate kinase,51,52 chorismate mutase,53,54 and chorismate synthase.55,56 Lawan et al.57 estimated the potential energy barrier for the conversion of 5-enolpyruvylshikimate-3-phosphate (EPSP) to chorismate catalysed by chorismate synthase. Their results showed that the enzyme stabilizes the negative charge of the flavin cofactor and EPSP substrate via hydrogen bonds with positively

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 22 of 34

22 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

charged arginine residues. Recently, computational studies have suggested that the stability, and consequently catalysis, can mostly be assigned to the energy cost associated with the reorganization of the environment in solution.47,58 Indeed, it is well known that electrostatic stabilization is the primary factor affecting the rate acceleration of the chemical reactions catalyzed by many enzymes.59 Krzeminska, Moliner, and Swiderek found that the electric field created by the protein in the active site of the protease from HIV-1 emerges as being critical to the electronic reorganization required during catalysis.60 Another work showed that the promiscuous activities of the phosphonate monoester hydrolases from Burkholderia caryophylli and Rhizobium leguminosarum are due to the compensatory electrostatic effects of different residues, allowing sufficient flexibility in the electrostatic environment of the active site to accommodate multiple substrates with distinct transition states and charge distributions.61 In this sense, it is important to emphasize that the catalytic effect is in fact determined by electrostatic preoganization,59 in which the electrostatic stabilization of the RS and TS in the protein is larger than that in water, thus leading to a smaller activation barrier in the enzymatic environment and the shift required for proton transfer in EPSP synthase is also a result of the electrostatic preorganization. Herein, to evaluate the key residues that contribute to the electrostatic stabilization of the TS in the protein environment, we perform an analysis of the electrostatic contributions of individual residues to the protein–substrate interactions (Figures S6 and S7). The origin of the substrate–enzyme interaction can be analyzed by decomposing the total interaction energy into a sum of the contributions due to each of the enzyme residues.62–64 At this point, note that only the substrate is included in the QM region in this analysis. The electrostatic contributions of individual residues exhibit an interaction with Glu341, this being the most important interaction in the addition mechanism, given that the residue makes substantial electrostatic contributions to the binding of the TS into an enzyme. Another important residue interacting with TS in the elimination mechanism is Asp313. These two residues are well positioned to stabilize any partial positive charge that develops during the reaction, which is in ACS Paragon Plus Environment

Page 23 of 34

The Journal of Physical Chemistry

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

agreement with Berti and Chindemi’s proposal that identified Asp313 and Glu341 as being key residues for stabilizing cationic intermediate/transition states in the AroA catalysis.15 Obviously, there are other important residues involved in the AroA catalysis to provide the key active site interactions that optimally stabilize the cationic intermediates/transition states for the enzyme reaction. However, Asp313 and Glu341 play large catalytic roles in the catalytic strategy of AroA, which is in agreement with the previous experimental data.15 We have also evaluated the influence of EPSP synthase on the electrostatic potential surfaces of the substrate along the reaction coordinates, determining the molecular electrostatic potential (MEP) surfaces for the substrate in both a protein and solution environment (see Figure S5 for all MEPs). MEP surfaces were obtained from the M06-2X/6-31G(d) level using single-point structures obtained from FESs. These surfaces correspond to an isodensity value of 0.002 au. The most nucleophilic regions (negative electronic potential) are shown in red, while the most electrophilic regions (positive electrostatic potential) are shown in blue. Our results show that regions of positive charge predominantly on C2 of the substrate are more positive in INT2 (see Figure S10 and S11), which corresponds to cationic intermediates (EPSP cations). This structure also displays more positive regions around the H atom of C3 of THI matching the position of Glu341 which can abstract the proton from this carbon (see SI). The hydrogen bound to C3 is ready to be transferred to Glu341, and can thus provide electrostatic stabilization acting as a general base accepting the proton from INT2. In addition, the positive charge on C2 of INT2 and TS4 matches the positions of Asp313 and Glu341, which is consistent with the experimental data that suggest that carboxylate side chains of Asp313 and Glu341 form an “electrostatic sandwich”15 to stabilize the cationic center at C2. It is important to note that the EPSP cations resemble the oxacarbenium ion intermediates (or oxacarbenium ion-like transition states) formed by the glycosylases65 that have long been targets of inhibitor design, and many cationic inhibitors that mimic the oxacarbenium ion have been characterized.66,67 This may provide a basis for the design of AroA inhibitors, since analogs of TS and INT can be used as inhibitors of EPSP. In ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 24 of 34

24 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

addition, note that INT2 and TS4 have a more positive charge around C2 in the protein than in a solution (see Figure 5). This positive potential on C2 is a consequence of both the electrostatic interaction between the enzyme and substrate and the electrostatic interaction with the substrate itself. Thus, as expected, EPSP synthase devotes its catalytic power to stabilizing the EPSP cationic intermediate and TS4. These results are in good agreement with the electrostatic preorganization concept19,59 that has been shown to be the primary factor accelerating the rate of the chemical reactions catalyzed by enzymes.

Figure 5. Structures along the 1D-PMF obtained from the QM/MM free energy simulations for addition and elimination. A) Structures named according to the locations along the 1D PMF shown in Figure 4A for the addition simulation and their respective MEP surfaces. B) Structures named according to the locations along the 1D PMF shown in Figure 4B for the elimination simulation and their respective MEP surfaces. The increase in the negative charges goes from blue (positive) to red (negative).

CONCLUSION In the present study, we took a QM/MM approach to the calculation of the free energy surfaces for the addition and elimination mechanisms used by EPSP synthase (AroA). An important observation for the FESs obtained for the reaction catalyzed by EPSP synthase is that the energetic profile for the addition mechanism is the exact converse of the elimination mechanism. Furthermore, the Glu341 was ACS Paragon Plus Environment

Page 25 of 34

The Journal of Physical Chemistry

25 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

identified as an acid/base catalytic residue for the addition and elimination mechanism. This conclusion is reinforced by the fact that the barrier for proton transfer from Glu341 to C3 of PEP is more energetically favorable than that from Asp313 to C3 of PEP. These results account for the experimental data that was not investigated in previous theoretical studies involving free energy surfaces. Considering the FES for the addition mechanism, any concerted transition state would be highly oxacarbenium ion-like, with C3 protonation almost complete and nucleophilic attack having just begun. In other words, the FES reveals that addition proceeds via a stepwise mechanism, starting with proton transfer (PT) from Glu341 to C3 of PEP, with its rate limiting step being one controlled by the PT. Mechanistically, the FES also suggests that the energetically preferable mechanism for elimination starts with the C-O breakdown of PEP during EPSP hydrolysis, after which the PT from PEP C3 to Glu341 was observed, with the preferred mechanism being a stepwise reaction with EPSP cation intermediate formation. Note also that the results obtained for the M06-2X/6-31G(d)/MM calculated activation free energy are in excellent agreement with the activation free energy value derived from the experimental estimation. In addition, the results also indicate a shift in the proton transfer from Glu341 to C3 for the reaction in enzyme and highlight the importance of the electrostatic interactions of Asp313 and Glu341, which form an “electrostatic sandwich”15 to stabilize the cationic center at C2. In summary, we have achieved the first theoretical energetic analysis based on a free energy surface for a reaction catalyzed by EPSP synthase. This study may be important for providing insight into enzymatic reactions and the design of EPSP synthase inhibitors.

Supporting Information Definition of the protonation models considered in this work and the PMFs calculated at AM1/d-phoT/MM level. Per-residue interaction energy contributions of individual amino acids. Map electrostatic potential (MEP) surfaces derived from M062X/6-31G* calculations for the RS, TS1 and INT1 polarized by the charges in the protein and water environments in the corresponding states. ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 26 of 34

26 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

Acknowledgments The authors thank the Coordenacão de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Conselho Nacional de Desenvolvimento Cientifico e Tecnológico (CNPq), Fundação de Amparo e Desenvolvimento da Pesquisa (FADESP) and Pró-reitoria de Pesquisa e Pós-graduação of Universidade Federal do Pará (PROPESP-UFPA) for their financial support.

REFERENCES (1)

Roberts, F.; Roberts, C. W.; Johnson, J. J.; Kyle, D. E.; Krell, T.; Coggins, J. R.; Coombs, G. H.; Milhous, W. K.; Tzipori, S.; Ferguson, D. J. P.; et al. Evidence for the Shikimate Pathway in Apicomplexan Parasites. Nature 1998, 393, 801–805.

(2)

Jana, S.; Paliwal, J. Novel Molecular Targets for Antimalarial Chemotherapy. Int. J. Antimicrob. Agents 2007, 30, 4–10.

(3)

Roberts, C. W.; Roberts, F.; Lyons, R. E.; Kirisits, M. J.; Mui, E. J.; Finnerty, J.; Johnson, J. J.; Ferguson, D. J. P.; Coggins, J. R.; Krell, T.; et al. The Shikimate Pathway and Its Branches in Apicomplexan Parasites. J. Infect. Dis. 2002, 185, S25–S36.

(4)

Clark, M. E.; Berti, P. J. Enolpyruvyl Activation by Enolpyruvylshikimate-3-Phosphate Synthase. Biochemistry 2007, 46, 1933–1940.

(5)

Anderson, K. S.; Johnson, K. A. Kinect and Structural-Analysis of Enzyme Intermediates Lessons from Epsp Synthase. Chem. Rev. 1990, 90, 1131–1149.

(6)

Byczynski, B.; Mizyed, S.; Berti, P. J. Nonenzymatic Breakdown of the Tetrahedral (AlphaCarboxyketal Phosphate) Intermediates of MurA and AroA, Two Carboxyvinyl Transferases. Protonation of Different Functional Groups Controls the Rate and Fate of Breakdown. J. Am. Chem. Soc. 2003, 125, 12541–12550.

(7)

Bondinell, W. E.; Vnek, J.; Knowles, P. F.; Sprecher, M.; Sprinson, D. B. Mechanism of 5ACS Paragon Plus Environment

Page 27 of 34

The Journal of Physical Chemistry

27 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

Enolpyruvylshikimate 3-Phosphate Synthetase. J. Biol. Chem. 1971, 246, 6191–6196. (8)

Bugg, T. D. H.; Walsh, C. T. Intracellular Steps of Bacterial Cell Wall Peptidoglycan Biosynthesis: Enzymology, Antibiotics, and Antibiotic Resistance. Nat. Prod. Rep. 1992, 9, 199–215.

(9)

Yang, Y.; Severin, A.; Chopra, R.; Krishnamurthy, G.; Singh, G.; Hu, W.; Keeney, D.; Svenson, K.; Petersen, P. J.; Labthavikul, P.; et al. 3,5-Dioxopyrazolidines, Novel Inhibitors of UDP-NAcetylenolpyruvylglucosamine Reductase (MurB) with Activity against Gram-Positive Bacteria. Antimicrob. Agents Chemother. 2006, 50, 556–564.

(10)

Lou, M. Y.; Burger, S. K.; Gilpin, M. E.; Gawuga, V.; Capretta, A.; Berti, P. J. Transition State Analysis of Enolpyruvylshikimate 3-Phosphate (EPSP) Synthase (AroA)-Catalyzed EPSP Hydrolysis. J. Am. Chem. Soc. 2012, 134, 12958–12969.

(11)

McConkey, G. A. Targeting the Shikimate Pathway in the Malaria Parasite Plasmodium Falciparum. Antimicrob. Agents Chemother. 1999, 43, 175–177.

(12)

Priestman, M. A.; Healy, M. L.; Becker, A.; Alberg, D. G.; Bartlett, P. A.; Lushington, G. H.; Schonbrunn, E. Interaction of Phosphonate Analogues of the Tetrahedral Reaction Intermediate with 5-Enolpyruvylshikimate-3-Phosphate Synthase in Atomic Detail. Biochemistry 2005, 44, 3241–3248.

(13)

Schramm, V. L. Transition States, Analogues, and Drug Development. Acs Chem. Biol. 2013, 8, 71–81.

(14)

Mizyed, S.; Wright, J. E. I.; Byczynski, B.; Berti, P. J. Identification of the Catalytic Residues of AroA (Enolpyruvylshikimate 3-Phosphate Synthase) Using Partitioning Analysis. Biochemistry 2003, 42, 6986–6995.

(15)

Berti, P. J.; Chindemi, P. Catalytic Residues and an Electrostatic Sandwich That Promote Enolpyruvyl Shikimate 3-Phosphate Synthase (AroA) Catalysis. Biochemistry 2009, 48, 3699– 3707. ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 28 of 34

28 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

(16)

Lou, M. Y.; Gilpin, M. E.; Burger, S. K.; Malik, A. M.; Gawuga, V.; Popovic, V.; Capretta, A.; Berti, P. J. Transition State Analysis of Acid-Catalyzed Hydrolysis of an Enol Ether, Enolpyruvylshikimate 3-Phosphate (EPSP). J. Am. Chem. Soc. 2012, 134, 12947–12957.

(17)

Warshel, A. Multiscale Modeling of Biological Functions: From Enzymes to Molecular Machines (Nobel Lecture). Angew. Chemie-International Ed. 2014, 53, 10020–10031.

(18)

Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme. J. Mol. Biol. 1976, 103, 227–249.

(19)

Lameira, J.; Bora, R. P.; Chu, Z. T.; Warshel, A. Methyltransferases Do Not Work by Compression, Cratic, or Desolvation Effects, but by Electrostatic Preorganization. Proteins Struct. Funct. Bioinforma. 2015, 83, 318–330.

(20)

Silva, J. R. A.; Govender, T.; Maguire, G. E. M.; Kruger, H. G.; Lameira, J.; Roitberg, A. E.; Alves, C. N. Simulating the Inhibition Reaction of Mycobacterium Tuberculosis L,DTranspeptidase 2 by Carbapenems. Chem. Commun. 2015, 51, 12560–12562.

(21)

Lameira, J.; Alves, C. N.; Tunon, I.; Marti, S.; Moliner, V. Enzyme Molecular Mechanism as a Starting Point to Design New Inhibitors: A Theoretical Study of O-GIcNAcase. J. Phys. Chem. B 2011, 115, 6764–6775.

(22)

Eschenburg, S.; Kabsch, W.; Healy, M. L.; Schonbrunn, E. A New View of the Mechanisms of UDP-N-Acetylglucosamine Enolpyruvyl Transferase (MurA) and 5-Enolpyruvylshikimate-3Phosphate Synthase (AroA) Derived from X-Ray Structures of Their Tetrahedral Reaction Intermediate States. J. Biol. Chem. 2003, 278, 49215–49222.

(23)

Nam, K.; Cui, Q.; Gao, J. L.; York, D. M. Specific Reaction Parametrization of the AM1/d Hamiltonian for Phosphoryl Transfer Reactions: H, O, and P Atoms. J. Chem. Theory Comput. 2007, 3, 486–504.

(24)

Huang, J.; MacKerell, A. D. CHARMM36 All-Atom Additive Protein Force Field: Validation ACS Paragon Plus Environment

Page 29 of 34

The Journal of Physical Chemistry

29 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

Based on Comparison to NMR Data. J. Comput. Chem. 2013, 34, 2135–2145. (25)

Sondergaard, C. R.; Olsson, M. H. M.; Rostkowski, M.; Jensen, J. H. Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pK(a) Values. J. Chem. Theory Comput. 2011, 7, 2284–2295.

(26)

Byrd, R. H.; Lu, P. H.; Nocedal, J.; Zhu, C. Y. A Limited Memory Algorithm for Bound Constrained Optimization. Siam J. Sci. Comput. 1995, 16, 1190–1208.

(27)

Kulik, H. J.; Zhang, J. Y.; Klinman, J. P.; Martinez, T. J. How Large Should the QM Region Be in QM/MM Calculations? The Case of Catechol O-Methyltransferase. J. Phys. Chem. B 2016, 120, 11381–11394.

(28)

Liao, R. Z.; Thiel, W. Convergence in the QM-Only and QM/MM Modeling of Enzymatic Reactions: A Case Study for Acetylene Hydratase. J. Comput. Chem. 2013, 34, 2389–2397.

(29)

Jindal, G.; Warshel, A. Exploring the Dependence of QM/MM Calculations of Enzyme Catalysis on the Size of the QM Region. J. Phys. Chem. B 2016, 120, 9913–9921.

(30)

Claeyssens, F.; Harvey, J. N.; Manby, F. R.; Mata, R. A.; Mulholland, A. J.; Ranaghan, K. E.; Schütz, M.; Thiel, S.; Thiel, W.; Werner, H. J. High-Accuracy Computation of Reaction Barriers in Enzymes. Angew. Chemie - Int. Ed. 2006, 45, 6856–6859.

(31)

Kamerlin, S. C. L.; Warshel, A. Multiscale Modeling of Biological Functions. Phys. Chem. Chem. Phys. 2011, 13, 10401–10411.

(32)

Kamerlin, S. C. L.; Haranczyk, M.; Warshel, A. Progress in Ab Initio QM/MM Free-Energy Simulations of Electrostatic Energies in Proteins: Accelerated QM/MM Studies of pKa, Redox Reactions and Solvation Free Energies†. J. Phys. Chem. B 2008, 113, 1253–1272.

(33)

Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem. Int. Ed. Engl. 2009, 48, 1198–1229.

(34)

Hu, H.; Yang, W. Free Energies of Chemical Reactions in Solution and in Enzymes with Ab Initio Quantum Mechanics/molecular Mechanics Methods. Annu. Rev. Phys. Chem. 2008, 59, ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 30 of 34

30 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

573–601. (35)

Torrie, G. M.; Valleau, J. P. Non-Physical Sampling Distributions in Monte-Carlo Free-Energy Estimation - Umbrella Sampling. J. Comput. Phys. 1977, 23, 187–199.

(36)

Roux, B. The Calculation of the Potential of Mean Force Using Computer Simulations. Comput. Phys. Commun. 1995, 91, 275–282.

(37)

Field, M. J. The pDynamo Program for Molecular Simulations Using Hybrid Quantum Chemical and Molecular Mechanical Potentials. J. Chem. Theory Comput. 2008, 4, 1151–1161.

(38)

Rivail, J. L.; Ruiz-Lopez, M. F.; Assfeld, X. Quantum Modeling of Complex Molecular Systems; Springer: Cham, 2015.

(39)

Lameira, J.; Kupchencko, I.; Warshel, A. Enhancing Paradynamics for QM/MM Sampling of Enzymatic Reactions. J. Phys. Chem. B 2016, 120, 2155–2164.

(40)

Plotnikov, N. V; Kamerlin, S. C. L.; Warshel, A. Paradynamics: An Effective and Reliable Model for Ab Initio QM/MM Free-Energy Calculations and Related Tasks. J. Phys. Chem. B 2011, 115, 7950.

(41)

Plotnikov, N. V. Computing the Free Energy Barriers for Less by Sampling with a Coarse Reference Potential While Retaining Accuracy of the Target Fine Model. J. Chem. Theory Comput. 2014, 10, 2987–3001.

(42)

Warshel, A.; Weiss, R. M. An Empirical Valence Bond Approach for Comparing Reactions in Solutions and in Enzymes. J. Am. Chem. Soc. 1980, 102, 6218–6226.

(43)

Bentzien, J.; Muller, R. P.; Florian, J.; Warshel, A. Hybrid Ab Initio Quantum Mechanics/Molecular Mechanics Calculations of Free Energy Surfaces for Enzymatic Reactions: The Nucleophilic Attack in Subtilisin. J. Phys. Chem. B 1998, 102, 2293–2301.

(44)

Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06 Functionals ACS Paragon Plus Environment

Page 31 of 34

The Journal of Physical Chemistry

31 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

and 12 Other Functionals. Theor. Chem. Acc. 2008, 119, 525. (45)

Hohenstein, E. G.; Chill, S. T.; Sherrill, C. D. Assessment of the Performance of the M05−2X and M06−2X Exchange-Correlation Functionals for Noncovalent Interactions in Biomolecules. J. Chem. Theory Comput. 2008, 4, 1996–2000.

(46)

Duarte, F.; Barrozo, A.; Åqvist, J.; Williams, N. H.; Kamerlin, S. C. L. The Competing Mechanisms of Phosphate Monoester Dianion Hydrolysis. J. Am. Chem. Soc. 2016, 138, 10664– 10673.

(47)

Garcia-Meseguer, R.; Marti, S.; Ruiz-Pernia, J. J.; Moliner, V.; Tunon, I. Studying the Role of Protein Dynamics in an SN2 Enzyme Reaction Using Free-Energy Surfaces and Solvent Coordinates. Nat. Chem. 2013, 5, 566–571.

(48)

Plotnikov, N. V; Warshel, A. Exploring, Refining, and Validating the Paradynamics QM/MM Sampling. J. Phys. Chem. B 2012, 116, 10342–10356.

(49)

Pu, J. Z.; Gao, J. L.; Truhlar, D. G. Multidimensional Tunneling, Recrossing, and the Transmission Coefficient for Enzymatic Reactions. Chem. Rev. 2006, 106, 3140–3169.

(50)

Anderson, K. S.; Johnson, K. A. Kinect Competence of the 5-Enolpyruvoylshukimate-3Phosphate Synthase Tetrahedral Intermediate. J. Biol. Chem. 1990, 265, 5567–5572.

(51)

Blanco, B.; Prado, V.; Lence, E.; Otero, J. M.; Garcia-Doval, C.; van Raaij, M. J.; Llamas-Saiz, A. L.; Lamb, H.; Hawkins, A. R.; González-Bello, C. Mycobacterium Tuberculosis Shikimate Kinase Inhibitors: Design and Simulation Studies of the Catalytic Turnover. J. Am. Chem. Soc. 2013, 135, 12366–12376.

(52)

Prado, V.; Lence, E.; Maneiro, M.; Vázquez-Ucha, J. C.; Beceiro, A.; Thompson, P.; Hawkins, A. R.; González-Bello, C. Targeting the Motion of Shikimate Kinase: Development of Competitive Inhibitors That Stabilize an Inactive Open Conformation of the Enzyme. J. Med. Chem. 2016, 59, 5471–5487.

(53)

Crespo, A.; Martí, M. A.; Estrin, D. A.; Roitberg, A. E. Multiple-Steering QM−MM Calculation ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 32 of 34

32 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

of the Free Energy Profile in Chorismate Mutase. J. Am. Chem. Soc. 2005, 127, 6940–6941. (54)

Zhang, X.; Zhang, X.; Bruice, T. C. A Definitive Mechanism for Chorismate Mutase. Biochemistry 2005, 44, 10443–10448.

(55)

Dmitrenko, O.; Wood

Harold B.; Bach, R. D.; Ganem, B. A Theoretical Study of the

Chorismate Synthase Reaction. Org. Lett. 2001, 3, 4137–4140. (56)

Rauch, G.; Ehammer, H.; Bornemann, S.; Macheroux, P. Mutagenic Analysis of an Invariant Aspartate Residue in Chorismate Synthase Supports Its Role as an Active Site Base. Biochemistry 2007, 46, 3768–3774.

(57)

Lawan, N.; Ranaghan, K. E.; Manby, F. R.; Mulholland, A. J. Comparison of DFT and Ab Initio QM/MM Methods for Modelling Reaction in Chorismate Synthase. Chem. Phys. Lett. 2014, 608, 380–385.

(58)

Ruiz-Pernia, J. J.; Marti, S.; Moliner, V.; Tunon, I. A Novel Strategy to Study Electrostatic Effects in Chemical Reactions: Differences between the Role of Solvent and the Active Site of Chalcone Isomerase in a Michael Addition. J. Chem. Theory Comput. 2012, 8, 1532–1535.

(59)

Warshel, A.; Sharma, P. K.; Kato, M.; Xiang, Y.; Liu, H.; Olsson, M. H. M. Electrostatic Basis for Enzyme Catalysis. Chem. Rev. 2006, 106, 3210–3235.

(60)

Krzeminska, A.; Moliner, V.; Swiderek, K. Dynamic and Electrostatic Effects on the Reaction Catalyzed by HIV-1 Protease. J. Am. Chem. Soc. 2016, 138, 16283–16298.

(61)

Barrozo, A.; Duarte, F.; Bauer, P.; Carvalho, A. T. P.; Kamerlin, S. C. L. Cooperative Electrostatic Interactions Drive Functional Evolution in the Alkaline Phosphatase Superfamily. J. Am. Chem. Soc. 2015, 137, 9061–9076.

(62)

Sousa, P. R. M.; De Alencar, N. A. N.; Lima, A. H.; Lameira, J.; Alves, C. N. Protein-Ligand Interaction Study of CpOGA in Complex with GlcNAcstatin. Chem. Biol. Drug Des. 2013, 81, 284–290.

(63)

Phipps, M. J. S.; Fox, T.; Tautermann, C. S.; Skylaris, C.-K. Energy Decomposition Analysis ACS Paragon Plus Environment

Page 33 of 34

The Journal of Physical Chemistry

33 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

Approaches and Their Evaluation on Prototypical Protein–drug Interaction Patterns. Chem. Soc. Rev. 2015, 44, 3177–3211. (64)

Reis, M.; Alves, C. N.; Lameira, J.; Tuñón, I.; Martí, S.; Moliner, V. The Catalytic Mechanism of Glyceraldehyde 3-Phosphate Dehydrogenase from Trypanosoma Cruzi Elucidated via the QM/MM Approach. Phys. Chem. Chem. Phys. 2013, 15, 3772–3785.

(65)

Zechel, D. L.; Withers, S. G. Glycosidase Mechanisms: Anatomy of a Finely Tuned Catalyst. Acc. Chem. Res. 2000, 33, 11–18.

(66)

Bols, M. 1-Aza Sugars, Apparent Transition State Analogues of Equatorial Glycoside Formation/cleavage. Acc. Chem. Res. 1998, 31, 1–8.

(67)

Lillelund, V. H.; Jensen, H. H.; Liang, X. F.; Bols, M. Recent Developments of Transition-State Analogue Glycosidase Inhibitors of Non-Natural Product Origin. Chem. Rev. 2002, 102, 515– 553.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 34 of 34

34 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

TOC Graphic

ACS Paragon Plus Environment