Conformational Substrate Selection Contributes to the Enzymatic

Nov 23, 2016 - Free-energy surfaces were obtained with the sum-hills tool implemented in the Plumed software and plotted with the matplotlib library i...
1 downloads 10 Views 4MB Size
Article pubs.acs.org/JPCB

Conformational Substrate Selection Contributes to the Enzymatic Catalytic Reaction Mechanism of Pin1 Esteban Vöhringer-Martinez* and Ciro Dörner Departamento de Físico-Química, Facultad de Ciencias Químicas, Universidad de Concepción, 4030000 Concepción, Chile S Supporting Information *

ABSTRACT: Conformational changes in enzymes and their role in their catalytic activity have been thoroughly addressed experimentally and theoretically. There is a vivid discussion in the field of enzyme catalysis on whether conformational changes of the enzyme are coupled to catalysis, or whether transition state stabilization through the preorganized protein and its electrostatic properties govern the catalysis. In this study, an additional contribution to the catalysis of one specific enzyme, Pin1, is proposed, which arises from its conformational selection of the peptide substrate from aqueous solution with the lowest activation barrier. The most stable conformers of the reactant (cis) and product (trans) peptide were identified through molecular dynamics simulations in combination with metadynamics. The cis−trans isomerization reaction was studied with quantum mechanical/molecular mechanical molecular dynamics simulations and density functional theory together with the mean reaction force, which allows us to separate structural and electronic contributions to the activation barrier. Our results show that enzyme Pin1 binds the trans isomer in the conformation of the peptide with the smallest activation barrier and reduces the barrier further through specific substrate−enzyme interactions, as we have shown previously. The activation barrier of the cis peptide is independent of the conformer and remains unchanged in the enzyme.



INTRODUCTION The mechanisms by which enzymes catalyze chemical reactions have been discussed extensively over the last few decades. The most accepted mechanism is the stabilization of the transition state by the preorganization of the enzyme cavity. In recent years, dynamical contributions to catalysis have been addressed but often with diverging definitions of the specific effect. A dynamical contribution is understood as the coupling of a specific enzyme mode with the reaction coordinate, which is different from aqueous solution and alters the Boltzmann distribution along the reaction coordinate. These couplings would violate the transition state assumption and contribute to the deviation of the transmission factor.1 However, until now no direct theoretical or experimental evidence for a nonBoltzmann distribution along the reaction coordinate has been found. Furthermore, the presence of a catalytic network has been proposed as an important factor to consider in enzyme catalysis.2−4 Enzymes and the enzyme−substrate complex undergo conformational changes due to their intrinsic flexibility, which built a network of conformations over which the reaction takes place. In this network picture, two extreme cases are discussed: (1) Conformational induction, where the substrate binds to the enzyme in one conformation that after a specific conformational transition leads to the reaction and (2) conformational selection in which one specific enzyme conformation binds the substrate facilitating the reaction. © XXXX American Chemical Society

Both were discussed as possible contributions to the mechanism by which enzymes catalyze chemical reactions.2,3 Here, an additional contribution to enzyme catalysis is proposed that combines both mechanisms and involves a conformational selection of the substrate by the enzyme. Our hypothesis applied on the specific enzyme Pin1 is that this enzyme selectively binds the substrate conformer in the trans form with the smallest activation barrier in aqueous solution and reduces the barrier further by transition state stabilization to yield the cis isomer, as we have shown in previous studies.5,6 It is anticipated that this contribution to enzyme catalysis is not general and might only be relevant for these types of unimolecular reactions that involve a cis−trans isomerization. The enzyme under study Pin1 is an enzyme that catalyzes the cis−trans isomerization of peptides with a phosphorylated threonine or serine-proline (pThr/pSer-Pro) motif7,8 (see Figure 1). From experimental studies in bulk water, it is known that at room temperature this process is very slow.9 The experimentally determined free-energy barrier for the cis to trans isomerization is 84.1 kJ/mol for peptides with the pThrPro motif.10 The catalytic mechanism of the isomerization reaction by Pin1 has recently been studied experimentally and theoretically. Received: September 11, 2016 Revised: November 22, 2016 Published: November 23, 2016 A

DOI: 10.1021/acs.jpcb.6b09187 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 1. Structure of the peptide undergoing cis−trans isomerization at the amide bond before the proline ring and the backbone dihedral angles ϕ and ψ of each residue together with the χ angle defining the relative position of the phosphate group used to identify conformers in aqueous solution.

barriers, which might have not been observed in the finite molecular dynamics simulation. From the obtained conformers, the reaction was followed restraining a validated reaction coordinate to obtain the free-energy profiles with QM/MM molecular dynamics simulations and density functional theory in aqueous solution for each conformer.

Important interactions among reactant, transition state, and product peptide analogues and residues in the catalytic cavity of Pin1 have been revealed by X-ray crystallography and NMR spectroscopy.11−17 Moreover, further insights were obtained from kinetic isotope studies. 18 From the theory side, accelerated molecular dynamics simulations have also contributed to the qualitative understanding of the reaction mechanism in Pin1 and in solution.19−23 In two recent studies, we were able to identify the role of specific residues in the catalytic reaction mechanism and provide an accurate freeenergy barrier along the reaction coordinate by extensive quantum mechanical/molecular mechanical (QM/MM) molecular dynamics simulations in combination with the mean reaction force (MRF) for the reaction in solution and in the enzyme.5,6 Our main conclusion was that the enzyme catalyzes the reaction from trans to cis, reducing the activation barrier to 20−30 kJ/mol, whereas the inverse reaction from cis to trans presented the same activation barrier as that in aqueous solution. In addition to the effect of direct substrate−enzyme interactions present in Pin1, specific intramolecular hydrogen bonds in the peptide between the nitrogen atom of the rotating amide bond and the N−H hydrogen atom of the next but one residue (phenylalanine in the studied peptide) have been proposed to reduce the activation barrier in the enzyme and in aqueous solution.24,25 The reduction is explained by the stabilization of the sp3 character of the nitrogen atom in the transition state. This possible contribution to catalysis, which could also be present in aqueous solution, is addressed in this study, comparing the activation barrier of conformations of the peptide with and without an intramolecular hydrogen bond. As the stabilization of the sp3 character of the nitrogen atom clearly requires the correct description of the electron distribution in the system, electronic structure methods in a QM/MM framework in combination with molecular dynamics are used. In previous molecular dynamics simulation studies that used common biomolecular force fields, these stabilizations and possible electronic changes in the reactions were not considered.19−23 To test our hypothesis that Pin1 selects one specific trans conformer from solution, which possesses a reduced activation barrier, extensive molecular dynamics simulations of the two configurations (cis and trans) of the peptide shown in Figure 1 in combination with a clustering method to identify the most stable conformers were employed. Additionally, metadynamics simulations with a specific combination of dihedral angles as collective variables identified conformers separated by small



METHODS

Molecular Dynamics Simulations. All molecular dynamics simulations were performed with the Gromacs 4.5.3 software package.26 For a detailed overview of the employed parameters, the reader is referred to our previous studies.5,6 The QM/MM molecular dynamics simulations were performed in combination with the ORCA package27 at the BLYP/SVP level of theory in an electrostatic embedding as incorporated in the Gromacs 4.5.3 package. In these simulations, the electrostatic interactions in the MM region were treated with the reaction field method (cutoff radius = 1.0 nm, switching radius = 0.9 nm, εr = 78) and the time step was reduced to 1 fs. The neighbor list was updated at every step. For the well-tempered metadynamics simulations, the Plumed software28 was used in combination with the Gromacs 4.6.5.26 Two collective variables represented by the ψ dihedral angle of T2p and Pro were considered in each simulation. To explore the optimal set of parameters and simulation time, first standard metadynamics were performed at a fixed value of σ of 0.35 radian and the height was varied between 0.1 and 0.5 kJ mol−1. With a height of 0.5 kJ mol−1 and a frequency of Gaussian deposition every 500 steps (dt = 2 fs), a complete sampling of the two degrees of freedom was achieved in 50 ns of simulation. Additionally, the amide ω angle was restrained at 0 or 180° to maintain the cis or trans configuration with a force constant of 500 kJ mol−1 rad−1. To obtain the final free energy as a function of these two collective variables, well-tempered metadynamics simulations with a bias factor of 10 were used. Free-energy surfaces were obtained with the sum-hills tool implemented in the Plumed software and plotted with the matplotlib library in python. System Setup. The initial structure of the phosphorylated peptide with sequence Ace-Gly-TPO-Pro-Phe-Gln-Nme was taken from the X-Ray structure (PDB-ID: 2Q5A) and solvated in a cubic box with 5 nm side lengths. Sodium and chloride ions were added to reach the physiological concentration of 154 mmol L−1 and neutralize the box. The SPC/E water model was used to describe water and the AMBER99sb force field for the peptide and the ions. For the phosphorylated threonine residue, the bonded parameters of Homeyer et al. were used and the B

DOI: 10.1021/acs.jpcb.6b09187 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 2. (A−C) Ramachandran plot of the backbone dihedral angles of the three relevant residues in the peptide A: T2p, B: Pro, C: Phe, colored as a function of the total simulation time in nanoseconds. D: χ angle of threonine phosphate reflecting the relative position of the phosphate group. At top and bottom, the stable conformers a−e in the cis configuration obtained from the analysis of the dihedral angles. The unique angles of each conformer are assigned in each Ramachandran plot together with the χ angle value. C

DOI: 10.1021/acs.jpcb.6b09187 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B atomic charges were derived using the D-HI method in aqueous solution, as described below. Atomic Charges for the Phosphorylated Threonine Residue. For the calculation of the atomic charges in the phosphorylated threonine residue in aqueous environment, QM/MM molecular dynamics simulations were performed in which the threonine residue up to the α-carbon was described quantum mechanically at the HF/6-31G(d) level of theory and the rest of the peptide with the Amber99sb force field. For the covalent bond connecting the QM α-carbon and the MM atoms, the link atom method was used. The method for the derivation of the atomic charges was the Hirshfeld-E method,29 as implemented in the HORTON 1.2.0 software package.30 To account for the solvent and peptide dynamics, which produce a varying electrostatic field around the threonine residue, the charges were obtained as an average from 3 ns molecular dynamics simulations. This average value was obtained by first equilibrating the peptide and the solvent employing the standard Amber99sb force field and the Sticht parameters for the phosphorylated threonine residue. Then, new charges were calculated “on the fly” in a QM/MM calculation and the old charges of the respective QM atoms were replaced before the continuation of the simulations. This procedure was repeated after a variable number of integration steps in molecular dynamics simulation. For details, see our previous paper.6 The dependence of the charges of the cis and trans configuration of the peptide was also considered, but no significant differences were observed. Finally, the convergence of the charges was investigated by varying the total simulation time from 500 ps to 3 ns. Reaction Mechanism. The reaction mechanism was studied starting from the different conformers of the cis configuration, and the reaction coordinate represented by the angle ζ (see Figure 1) was followed. This dihedral angle ζ is defined by the Cα atom and the oxygen atom in the carbonyl group of the T2p residue and the Cδ and nitrogen atom of the Pro residue. As described in our previous studies, this reaction coordinate is more appropriate than the amide bond defined along the peptide backbone because it is decoupled from the pyramidalization of the amide nitrogen atom.5,6 It is important to note that the rotation was first followed at the force-field level, equilibrating the system 5 ns before 10° of rotation. This results in a proper relaxation of the perpendicular degrees of freedom and proper sampling. From the simulations at the force-field level, three to six structures were taken from the last 4 ns along the reaction coordinate and served as input for the QM/MM simulations. The QM region was the same as in our previous study, involving all atoms in the amino acid sequence starting at the Cα atom of the phosphorylated threonine residue until the N− H group of phenylalanine (in total 20 atoms). Each window was simulated 15 ps from which the last 10 were used to obtain the MRF. The QM calculations were performed at the BLYP/ SVP level of theory, employing the link atom method in an electrostatic embedding. The reaction mechanism was studied by the MRF, which identifies structural and electronic changes along the reaction coordinate and was used in the enzyme reaction before.5,6,31 Briefly, the MRF is defined as the negative derivative of the free energy along the reaction path (at constant pressure and temperature)

⟨F(ξ)⟩ = −

dG dξ

(1)

It divides the reaction path by its minimum and maximum into three regions: the reactant, transition state, and product regions. The reactant and product regions involve mostly structural rearrangements of the atoms to reach the transition state from either the reactant or the product, whereas the transition state region represents mainly changes in the electronic configuration required to make or break chemical bonds. These regions allow the activation energy ΔG⧧ to be divided into two contributions: ΔG⧧ = [G(ξTS) − G(ξR)] = W1 + W2. W1 is mostly structural, and W2 encompasses mostly electronic contributions. By analogy, for the reverse reaction, the activation barrier is separated into a structural contribution W4 and an electronic one W3. As in our previous study, the MRF is obtained with the Umbrella Integration method of Kästner et al.32,33



RESULTS AND DISCUSSION To test our hypothesis of a substrate conformational selection for catalysis by Pin1, first the different conformers of the peptide in aqueous solution in the cis and trans configurations were studied by molecular dynamics simulations and welltempered metadynamics. Identification of Relevant Conformers in Aqueous Solution. The cis configuration of the peptide taken from the crystal structure co-crystallized with the enzyme Pin1 (PDB-ID: 2Q5A) was simulated in aqueous solution under physiological conditions with a deprotonated phosphate group for 1 μs at 298 K and 1 bar (experimentally, only a minor difference of 2− 4 kJ/mol was found for the activation barrier of the partial protonated state). To identify the different conformers of the cis configuration, the backbone dihedral angles (see Figure 1) of the three residues next to the rotating amide bond were analyzed: phosphorylated threonine (T2p), proline (Pro), and phenylalanine (Phe). The dihedral angles of the other residues were omitted because their rotation is limited by friction with the solvent rather by intramolecular interactions. Additionally, also the χ angle (see Figure 1), defining the rotation of the phosphate group with minus two charge, was analyzed to monitor its relative orientation to the rotating amide bond (see Figure 2D). As shown in Figure 2C,D during the total simulation time, the phosphate group and the dihedral angle of Phe display the largest fluctuations, indicating fast transitions between different conformers. The ψ angle of Pro changes from 0 to 150° at the beginning of the simulation and returns to the original value twice, indicating the presence of different conformers separated by larger barriers. As the starting structure was obtained from the crystal structure in the enzyme cavity of Pin1, the change to another ψ value of Pro might indicate that the conformer in the enzyme is not necessarily the most stable one in aqueous solution. The other backbone angle ϕ in Pro is not free to rotate due to the proline ring, which fixes its value at ϕ ≃ −75°. For T2p, only small changes are observed especially in the ϕ angle, which changes to a value of ϕ ≃ −150° from the original value of −60°. These changes are correlated with changes in the ψ angle of Pro and correspond to one conformer. After a careful analysis, five representative conformations a−e were chosen, which had a unique combination of values for the D

DOI: 10.1021/acs.jpcb.6b09187 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 3. Free-energy surfaces obtained from the metadynamics simulations for the cis and trans configuration as a function of the two collective variables: the ψ angle of phosphorylated threonine (T2p) and proline (Pro).

of all conformations of the cis configuration in the equilibrium simulation. These conformations are further divided into two groups differing by the ψ angle of Pro, and in each group, minor differences are attributed to the position of the phosphate group or the backbone dihedral angles of Phe. From the above results, one can conclude that the ψ angle of the proline residue distinguishes between the two main conformers. In our previous study of the isomerization reaction in the enzyme, the relevance of the rotation of the ψ angle in the preceding phosphorylated threonine residue (T2p) in the reaction mechanism was shown.6 The populations of the two main conformers found above, however, may be biased by the total simulation time and the total transitions observed, and conformers divided by larger free-energy barriers might not have been observed at all due to initial conditions. Therefore, well-tempered metadynamics simulations have been carried out with the ψ angle of Pro and T2P as collective variables to obtain the free-energy map for the cis isomer (for details see the Method section). The free-energy landscape is shown on the left of Figure 3 as a function of the ψ angles of the phosphorylated threonine and proline residue. In the free-energy landscape, two minima (blue regions) associated with the ψ backbone angle of proline are identified, which are separated by a small barrier of ≃10 kJ/mol. These two conformers have almost the same free energy although the one with the proline ψ angle at a value larger than 90° is slightly more stable than the one with ψ = 0°, which presents the intramolecular hydrogen bond. This is in agreement with the relative population of the two clusters described above and reflects the probability distribution of each conformer observed in the analysis of the dihedral angles in Figure 2. There is no minimum for values of the ψ backbone angle of T2p smaller than 90° or π , as shown in Figure 3 in this 2 configuration. In summary, the metadynamics simulations confirmed the conformations observed in the equilibrium simulation and their abundance in aqueous solution. Metadynamics have also been applied to study the freeenergy landscape of the trans configuration, which represents

shown angles (for the angle combination of each conformer, see Figure S1 in the Supporting Information). The atoms constituting the peptide backbone of the three residues and the phosphate group of these five conformers a−e (in the order as they appear in the simulation) are shown in Figure 2 at the top and bottom (other atoms of the peptide, solvent, and ions were omitted for clarity). Comparing the five structures and the dihedral angles, two groups are identified: the first formed by a and b, which differ in the position of the phosphate group and the backbone angles of the phenylalanine residue, and the second group by c−e, which present a ψ angle of Pro associated with no intramolecular hydrogen bond to the nitrogen atom of the amide bond to be rotated. Therefore, the first group with the intramolecular hydrogen bond shares values of ψ = 0 of Pro and the second group missing this type of interaction values of 180°. The other three structures of the second group c−e differ by the position of the phosphate group and the backbone angles of Phe, which present small rotational barriers reflected in several transitions in the dihedral analysis. The obtained conformers were compared with a cluster analysis based on the root-mean-square deviation values of the backbone atoms, including the carbonyl oxygen atoms as proposed by Daura et al.34 This method has been used before to identify conformations in hexapeptides. The analysis of the 1 μs trajectory resulted in two clusters, which only differed in the ψ angle of the proline residue and resembled the two groups described above. The cluster with the largest population (69%) was the one with ψ = 130° not forming the intramolecular hydrogen bond (structures c−e) followed by the second cluster with a population of 25%, which corresponds to the conformation with ψ = 0° of the proline residue forming the hydrogen bond (structures a and b). This second lesspopulated conformer is the one which binds in the cis configuration to the Pin1 enzyme as evidenced by the crystal structure. Together, these two clusters made 95.5% of the total population, and therefore, represent most of the conformations (see Figure S2 for more details on the cluster distribution). In summary, the cluster and the dihedral angle analysis identify the five conformations in Figure 2 to be representative E

DOI: 10.1021/acs.jpcb.6b09187 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B the product of the isomerization reaction. The right-hand side of Figure 3 shows that the free-energy landscape changes considerably because minima are also observed for negative values of the ψ backbone angle of T2p. There is even a pronounced minimum at ψT2p ≃ −45° and ψpro ≃ 0°. This conformation belongs to a compact trans conformer, which presents an intramolecular hydrogen bond and is the one that binds to Pin1 as evidenced by crystal structures of peptide analogues11(see Figure 5). The conformation with ψT2p ≥ 90° corresponds to extended structures, in which one group with values of ψpro ≃ 0° presents an intramolecular hydrogen bond contrary to the ones with larger values (see e.g., Figure 5). As evidenced in the next section, all of these conformations of the trans configuration obtained in the metadynamics simulations were also observed as products after the isomerization reaction. Cis−Trans Isomerization Reaction in Aqueous Solution. As shown schematically in Figure 4, the molecular

Figure 5. Products of the isomerization reactions in the trans configuration obtained from conformer c in Figure 2 (compact) and from conformer a in Figure 2 (extended).

(see Figure 2) was the starting structure (for more details, Figure S3 shows the complete set of dihedral angles of all three residues). Starting from this conformation (ζ ≃ −20°) toward positive values of the reaction coordinate, a change in the ψ angle of T2p and Pro is observed. The structure of the obtained product in the trans configuration is shown in Figure 5 on the left-hand side. This trans conformer corresponds to the compact form characterized in the metadynamics simulations above with a negative ψ angle of T2p and Pro. Interestingly, to obtain this conformer in the isomerization reaction, the ψ angle of T2p has to perform a concerted rotation with the amide bond. The ψ angle of T2p starts with positive values in the cis configuration and changes in a concerted manner with the reaction coordinate to reach the required negative values in the product. One has to note that the same concerted motion is observed if the simulation is started from the equilibrated trans conformer (total simulation time 1 μs) to reach the reactant cis configuration. The rotation of all other conformers in the cis configuration shown in Figure 2 resulted in extended conformations in the trans configuration, as shown representatively in Figure 5 for starting structure a in the cis configuration (the evolution of the backbone angles of the three residues as a function of the reaction coordinate for this conformer is shown in Figure S4). The obtained trans conformer was also observed in the metadynamics simulations and corresponded to the minima at positive ψ values of T2p and positive or negative ψ angles in Pro. When the cis isomer is rotated in the opposite direction toward negative values of the ζ angle, all conformations lead to an extended conformation. As this rotation direction is impeded in the enzyme and we have shown previously that the activation barrier is larger in aqueous solution,5 it was not considered further. To include electronic effects as the polarization of the electron density due to the surroundings or intramolecular and intermolecular hydrogen bonds, the isomerization reaction of the representative cis conformers, a and c, was also studied with QM/MM molecular dynamics simulations at the BLYP/SVP level of theory using the Umbrella Integration method.32,33 This method yields the derivative of the free energy at defined values of the reaction coordinate, which resembles the MRF. At the top of Figure 7, the MRF is shown on the left for the rotation of cis conformer a to yield the extended trans conformer and the rotation of cis conformer c to yield the compact trans conformer on the right. The MRF is zero at a value of −20 for the ζ angle corresponding to the stable cis reactant state and decreases until its minimum to become zero at the transition state. After the transition state, the MRF increases to positive values at its maximum and becomes zero

Figure 4. Scheme of the simulation protocol used to study the isomerization reaction.

dynamics simulations of the cis configuration and the subsequent dihedral and cluster analysis resulted in the representative conformers a−e, which were further confirmed by metadynamics simulations. These conformations of the cis configuration in aqueous solution have been used as starting structures to rotate the ζ angle, which is independent of the nitrogen pyramidalization, as shown by us and others.5,6 Details on the procedure are found in the Method section, but it should be emphasized that to observe possible conformational changes along the reaction coordinate, the rotation was first performed at the force-field level with a proper equilibration at each restrained value of the reaction coordinate. From these simulations, three to six representative structures were used to carry out the QM/MM molecular dynamics simulations at the BLYP/SVP level of theory, avoiding possible deficiencies of the force field and accounting for the polarization of the electron density. The isomerization at the force-field level yields structures, which could be assigned to a compact and extended form, as shown in Figure 5 on the left and right-hand sides, respectively. Representatively, Figure 6 displays the change in the relevant dihedral angles of T2p and Pro as a function of the reaction coordinate, where the conformation, c, in the cis configuration F

DOI: 10.1021/acs.jpcb.6b09187 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 6. Average values of the peptide backbone dihedral angles of phosphorylated threonine and proline residues as a function of the reaction coordinate. In blue, the ψ angle is shown and in red, the ϕ angle.

The reaction free-energy is positive for both conformers, which corresponds to a more stable cis rather than a trans isomer. This is in disagreement with the experimental observation of only some percent of cis isomer in aqueous solution, which corresponds to ≃7 kJ/mol more stable trans isomer. The shaded area in the figure, however, clearly shows that it is not possible to obtain this accuracy through the integration of the MRF, which also includes an inherent integration error. To improve the free energy between the two configurations would require more intermediate values for the reaction coordinate with 3−5 QM/MM simulations for each value, which would result in a pronounced computational cost. Finally, one has to stress that to reproduce the experimental free-energy difference between the isomers, the relative free energy of all conformers would also be necessary, which would involve sampling all conformations at the QM/MM level, which is clearly out of reach. For the reverse process from trans to cis, the barrier differs considerably, yielding values that are ≃30 kJ/mol smaller for the reaction that involves the trans conformer in a compact form, as shown in Figure 5. Considering that the compact form is the only one which binds to the enzyme due to steric hinderance of the extended form, one can conclude that the enzyme selects the conformation in the trans configuration with the smallest barrier and reduces this barrier further, as shown in our previous studies.5,6 The MRF also permits elucidating why the barrier of the extended trans configuration is larger than that of the compact form. The absolute values of the MRF are larger between angles of 90 and 180° for the extended trans configuration (left-hand side) than those for the compact conformer, which implies that less free energy is needed to reach the transition state from the latter conformer. Additionally, the MRF naturally divides the reaction in processes associated with structural changes between the reactant and its minimum, and the product and its maximum, and processes where the electronic configuration

again at the product trans configuration. The error bars for each point are obtained from 3 to 5 independent simulations considering the standard deviation of the mean angle displacements. A comparison of the MRF for the two conformers shows that the MRF from the cis configuration until the transition state is very similar with small variations in the position of their minima. After the transition state however, the MRF of the reaction that yields the compact trans isomer (right-hand side in Figure 7) reveals smaller values with a maximum of the MRF displaced to the products. From these results, one concludes that the combined rotation of the ζ angle and the ψ angle of T2p observed for conformer c above and associated with the formation of the compact trans isomer reduces the value of the MRF and therefore should also result in a smaller free energy of activation when the reverse reaction from trans to cis is considered. To obtain the free energy as a function of the reaction coordinate, the MRF was integrated over the entire rotation. This numerical integration has to consider the error of the derivative at each point, which results in the shaded area of the free-energy values and is obtained by integration of the minimal and maximal values of the MRF. The free energy shows that the cis to trans reaction for both conformers a and c presents similar barriers in the range of 110 kJ/mol, which are close to the experimental value of 84.1 kJ/mol. It should be noted that the conformers a and c in the cis configuration differ by the presence of the intramolecular hydrogen bond present in a but not in c. Because both conformers’ activation barriers differ only by ≃10 kJ/mol and the MRF values are similar, one can conclude that the intramolecular hydrogen bond has no or little effect on the activation barrier for these types of peptides. One possible explanation for the negligible effect of the intramolecular hydrogen bond is that the sp3 character of the nitrogen atom at the transition state in aqueous solution is also stabilized by hydrogen bonds to water molecules diminishing the role of the intramolecular hydrogen bond in the reaction. G

DOI: 10.1021/acs.jpcb.6b09187 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 7. Top: MRF along the reaction coordinate for the isomerization reaction of the cis conformers a (left) and c (right) in aqueous solution at the BLYP/SVP level of theory. Bottom: Integrated free-energy profile along the reaction coordinate with shaded area representing the error obtained by integration of the MRF considering the error bars.

of the system is altered, which lie in the transition state region between the two extrema of the MRF. These regions along the reaction coordinate separate the free-energy barrier in structural contributions (W1 and W4) and electronic ones (W2 and W3) for each direction of the reaction. Table 1 displays the contribution to the free-energy barrier for the reaction starting with conformer a in the cis configuration and forming the

extended trans product and the c conformer forming the compact trans configuration. The structural contributions to the activation barrier from the cis configuration outweigh the electronic ones, W1 and W2 respectively, for both conformers a and c. The structural contribution of conformer a is slightly larger, which can be ascribed to an intramolecular hydrogen bond between the backbone atoms of T2p and Phe (see Figure 2). This hydrogen bond is not present in the other conformer and might explain the increase in the structural contribution by 20 kJ/mol. The electronic contributions of the two conformers are the same evidencing no influence of the intramolecular hydrogen bond between the amino group of Phe and the nitrogen atom of the proline ring. If the intramolecular hydrogen bond stabilizes the sp3 character of the nitrogen atom at the transition state, then

Table 1. Contributions to the Activation Barrier for the Cis− Trans and Trans−Cis Reaction in kJ/mol cis conformer

W1

W2

W3

W4

trans conformer

a c

80 64

41 42

−50 −48

−43 −18

extended compact H

DOI: 10.1021/acs.jpcb.6b09187 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

In light of our results and the ones on Cyclophilin A, one can conclude that cis−trans isomerization reactions may be catalyzed by a conformational substrate selection together with an electrostatic stabilization of the transition state. Finally, one has to mention that conformational substrate selection in enzyme catalysis has also been proposed for other types of unimolecular reactions such as the prephenate reaction catalyzed by chorismate mutase.3

one would expect a smaller electronic contribution in c, which presents this type of interaction, than in a. As this is not observed, autocatalysis is ruled out for the reaction in aqueous solution, and the relevance of the surrounding water molecules, which would stabilize both conformers to the same amount and lead to the same electronic contributions for both conformers, is reinforced. For the reverse reaction starting from the trans configuration, a smaller structural contribution W4 is clearly identified for the compact trans configuration formed from the cis conformer c compared to the extended trans conformer obtained after the rotation of cis conformer a. The electronic contribution to the activation barrier is the same for both trans conformers. The origin of the reduced structural contribution lies in the closed form of the peptide, which is maintained during the rotation. This is accomplished by rotating the ψ angle of T2p in a concerted manner with the amide bond in a jump-rope type mechanism (see Figure 5). The reduced structural contribution can also be understood as a smaller entropy difference between the compact trans conformer and the minimum of the MRF compared to the extended one, which loses more flexibility to reach the same value of the reaction coordinate. Finally, the obtained results have to be considered in the context of the reaction taking place in the enzyme. The crystal structure of a trans-locked alkene peptide from Zhang et al.11 shows that Pin1 is only able to bind the compact form in the trans configuration. Therefore, our hypothesis of a conformational selection of the trans configuration is confirmed because the enzyme selects the conformer with the smallest activation free-energy barrier from aqueous solution. Furthermore, our previous results5,6 show that the enzyme reduces the barrier for the trans to cis isomerization further by direct interactions with enzyme residues such as Ser-154. However, the net overall catalytic effect of the enzyme on the reaction can only be addressed if additionally the binding free energies of the cis and trans conformations are considered. If Pin1 catalyzes the reaction by a combination of conformational selection and reduction of the activation barrier by specific substrate−enzyme interactions, then one may ask whether this mechanism also holds for other enzymes, which catalyze cis to trans isomerizations. One extensively studied representative is Cyclophilin A, which has been studied experimentally and theoretically. A recent study of Leone et al. applied several computational methods to study the reaction in the enzyme and aqueous solution concluding that the enzyme specifically catalyzes one conformer of all of the different ones present in solution.35 As they used force-fieldbased molecular dynamics simulations, no quantitative difference between the activation barriers of the different conformers was identified. In a more recent study by Camilloni et al., a combination of molecular dynamics simulation and NMR studies on Cyclophilin A proposed a mechanism based on the electric field formed by the enzyme, which acts on the dipole moment of the rotating carbonyl bond.36 This handle mechanism is perfectly compatible with the specific Ser-154 interactions found in Pin1, favoring the reaction through the jump-rope mechanism. The previous computational study of Leone et al. and the recent results of Camilloni, therefore, suggest that in Cyclophilin A also both mechanisms of conformational selection and specific interactions providing the correct electric field are responsible for the observed catalysis in this enzyme.



CONCLUSIONS To test our hypothesis of a conformational selection by Pin1, molecular dynamics simulations in combination with metadynamics were used to identify the stable conformers of the peptide substrate in aqueous solution. For each conformation in the cis configuration of the peptide, the isomerization reaction was followed along a predefined reaction coordinate with the umbrella integration technique in the QM/MM framework to account for electronic changes. It was found that for the cis configuration, the barriers are independent of the conformation used and no catalytic influence of an intramolecular hydrogen bond proposed previously was observed. For the reverse direction starting with the trans configuration, the compact conformation that binds to the enzyme presented a lower activation free-energy barrier than the extended conformation. The reduced barrier was explained through the MRF as a smaller structural contribution that is needed to reach the transition state or the minimum of the MRF. Altogether, the results suggest that Pin1 selectively binds the trans conformer with the lowest barrier, which contributes to the catalytic mechanism of this enzyme together with the stabilization of the transition state through specific peptide− enzyme interactions.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b09187. Figure S1 displays the peptide backbone dihedral angles of the relevant residues in the equilibrium simulations of the cis isomer; Figure S2 the population distribution of the clusters; Figure S3 the same as S1 but for the trans isomer; and Figures S4 and S5 the average value of all dihedral angles during the isomerization reaction yielding the compact and extended conformation (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: 0056-41-2204986. ORCID

Esteban Vöhringer-Martinez: 0000-0003-1785-4558 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors acknowledge funding support provided by FONDECYT 1160197 and Millenium Nucleus NC No.120082. Computational resources at the Southern GPU Cluster - Fondequip EQM150134 are acknowledged. I

DOI: 10.1021/acs.jpcb.6b09187 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B



(20) Velazquez, H. A.; Hamelberg, D. Conformational selection in the recognition of phosphorylated substrates by the catalytic domain of human Pin1. Biochemistry 2011, 50, 9605−9615. (21) Barman, A.; Hamelberg, D. Cysteine-mediated dynamic hydrogen-bonding network in the active site of Pin1. Biochemistry 2014, 53, 3839−3850. (22) Velazquez, H. A.; Hamelberg, D. Dynamical role of phosphorylation on serine/threonine-proline Pin1 substrates from constant force molecular dynamics simulations. J. Chem. Phys. 2015, 142, No. 075102. (23) Di Martino, G.; Masetti, M.; Cavalli, A.; Recanatini, M. Mechanistic insights into Pin1 peptidyl-prolyl cis-trans isomerization from umbrella sampling simulations. Proteins: Struct., Funct., Bioinf. 2014, 82, 2943−2956. (24) Cox, C.; Young, V.; Lectka, T. Intramolecular catalysis of amide isomerization. J. Am. Chem. Soc. 1997, 119, 2307−2308. (25) Yonezawa, Y.; Nakata, K.; Sakakura, K.; Takada, T.; Nakamura, H. Intra- and Intermolecular Interaction Inducing Pyramidalization on Both Sides of a Proline Dipeptide during Isomerization: An Ab Initio QM/MM Molecular Dynamics Simulation Study in Explicit Water. J. Am. Chem. Soc. 2009, 131, 4535−4540. (26) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (27) Neese, F. The ORCA program system. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 2, 73−78. (28) Tribello, G.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185, 604−613. (29) Verstraelen, T.; Ayers, P. W.; Van Speybroeck, V.; Waroquier, M. Hirshfeld-E Partitioning: AIM Charges with an Improved Trade-off between Robustness and Accurate Electrostatics. J. Chem. Theory Comput. 2013, 9, 2221−2225. (30) Fischer, S.; Dunbrack, R. L.; Karplus, M. Cis-Trans Imide Isomerization of the Proline Dipeptide. J. Am. Chem. Soc. 1994, 116, 11931−11937. (31) Vöhringer-Martinez, E.; Toro-Labbe, A. The mean reaction force: a method to study the influence of the environment on reaction mechanisms. J. Chem. Phys. 2011, 135, No. 064505. (32) Kästner, J.; Thiel, W. Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method: ”Umbrella integration”. J. Chem. Phys. 2005, 123, No. 144104. (33) Kästner, J.; Thiel, W. Analysis of the statistical error in umbrella sampling simulations by umbrella integration. J. Chem. Phys. 2006, 124, No. 234106. (34) Daura, X.; van Gunsteren, W.; Mark, A. Folding−unfolding thermodynamics of a β-heptapeptide from equilibrium simulations. Proteins: Struct., Funct., Bioinf. 1999, 34, 269−280. (35) Leone, V.; Lattanzi, G.; Molteni, C.; Carloni, P. Mechanism of action of cyclophilin a explored by metadynamics simulations. PLoS Comput. Biol. 2009, 5, No. e1000309. (36) Camilloni, C.; Sahakyan, A.; Holliday, M.; Isern, N.; Zhang, F.; Eisenmesser, E.; Vendruscolo, M. Cyclophilin A catalyzes proline isomerization by an electrostatic handle mechanism. Proc. Natl. Acad. Sci. U.S.A. 2014, 111, 10203−10208.

REFERENCES

(1) Kamerlin, S. C. L.; Warshel, A. At the dawn of the 21st century: Is dynamics the missing link for understanding enzyme catalysis? Proteins: Struct., Funct., Bioinf. 2010, 78, 1339−1375. (2) Ma, B.; Nussinov, R. Enzyme dynamics point to stepwise conformational selection in catalysis. Curr. Opin. Chem. Biol. 2010, 14, 652−659. (3) Giraldo, J.; Roche, D.; Rovira, X.; Serra, J. The catalytic power of enzymes: Conformational selection or transition state stabilization? FEBS Lett. 2006, 580, 2170−2177. (4) Wei, G.; Xi, W.; Nussinov, R.; Ma, B. Protein Ensembles: How Does Nature Harness Thermodynamic Fluctuations for Life? The Diverse Functional Roles of Conformational Ensembles in the Cell. Chem. Rev. 2016, 116, 6516−6665. (5) Vöhringer-Martinez, E.; Duarte, F.; Toro-Labbe, A. How Does Pin1 Catalyze the Cis-Trans Prolyl Peptide Bond Isomerization? A QM/MM and Mean Reaction Force Study. J. Phys. Chem. B 2012, 116, 12972−12979. (6) Vöhringer-Martinez, E.; Verstraelen, T.; Ayers, P. W. The influence of Ser-154, Cys-113, and the phosphorylated threonine residue on the catalytic reaction mechanism of Pin1. J. Phys. Chem. B 2014, 118, 9871−9880. (7) Rahfeld, J.-U.; Schierhorn, A.; Mann, K.; Fischer, G. A novel peptidyl-prolyl cis/trans isomerase from Escherichia coli. FEBS Lett. 1994, 343, 65−69. (8) Yaffe, M. B.; Schutkowski, M.; Shen, M.; Zhou, X. Z.; Stukenberg, P. T.; Rahfeld, J.-U.; Xu, J.; Kuang, J.; Kirschner, M. W.; Fischer, G.; et al. Sequence-specific and phosphorylation-dependent proline isomerization: a potential mitotic regulatory mechanism. Science 1997, 278, 1957−1960. (9) Dugave, C.; Demange, L. Cis-trans isomerization of organic molecules and biomolecules: Implications and applications. Chem. Rev. 2003, 103, 2475−2532. (10) Schutkowski, M.; Bernhardt, A.; Zhou, X.; Shen, M.; Reimer, U.; Rahfeld, J.; Lu, K.; Fischer, G. Role of Phosphorylation in Determining the Backbone Dynamics of the Serine/Threonine-Proline Motif and Pin1 Substrate Recognition. Biochemistry 1998, 37, 5566−5575. (11) Zhang, M.; Wang, X.; Chen, X.; Bowman, M.; Luo, Y.; Noel, J.; Ellington, A.; Etzkorn, F.; Zhang, Y. Structural and kinetic analysis of prolyl-isomerization/phosphorylation cross-talk in the CTD code. ACS Chem. Biol. 2012, 7, 1462−1470. (12) Namanja, A. T.; Wang, X. J.; Xu, B.; Mercedes-Camacho, A. Y.; Wilson, K. A.; Etzkorn, F. A.; Peng, J. W. Stereospecific gating of functional motions in Pin1. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 12289−12294. (13) Namanja, A. T.; Peng, T.; Zintsmaster, J.; Elson, A.; Shakour, M.; Peng, J. Substrate recognition reduces side-chain flexibility for conserved hydrophobic residues in human Pin1. Structure 2007, 15, 313−327. (14) Xu, G.; Slebodnick, C.; Etzkorn, F. Cyclohexyl ketone inhibitors of Pin1 dock in a trans-diaxial cyclohexane conformation. PLoS One 2012, 7, No. e44226. (15) Xu, G.; Zhang, Y.; Mercedes-Camacho, A.; Etzkorn, F. A reduced-amide inhibitor of Pin1 binds in a conformation resembling a twisted-amide transition state. Biochemistry 2011, 50, 9545−9550. (16) Labeikovsky, W.; Eisenmesser, E. Z.; Bosco, D. A.; Kern, D. Structure and Dynamics of Pin1 During Catalysis by NMR. J. Mol. Biol. 2007, 367, 1370−1381. (17) Eichner, T.; Kutter, S.; Labeikovsky, W.; Buosi, V.; Kern, D. Molecular Mechanism of Pin1-Tau Recognition and Catalysis. J. Mol. Biol. 2016, 428, 1760−1775. (18) Mercedes-Camacho, A. Y.; Mullins, A.; Mason, M.; Xu, G.; Mahoney, B.; Wang, X.; Peng, J.; Etzkorn, F. Kinetic isotope effects support the twisted amide mechanism of Pin1 peptidyl-prolyl isomerase. Biochemistry 2013, 52, 7707−7713. (19) Velazquez, H. A.; Hamelberg, D. Conformation-Directed Catalysis and Coupled Enzyme-Substrate Dynamics in Pin1 Phosphorylation Dependent Cis-Trans Isomerase. J. Phys. Chem. B 2013, 117, 11509−11517. J

DOI: 10.1021/acs.jpcb.6b09187 J. Phys. Chem. B XXXX, XXX, XXX−XXX