Mechanism of Action of Cyclophilin A Explored by Metadynamics Simulations

Trans/cis prolyl isomerisation is involved in several biological processes, including the development of numerous diseases. In the HIV-1 capsid protein (CA), such a process takes place in the uncoating and recruitment of the virion and is catalyzed by cyclophilin A (CypA). Here, we use metadynamics simulations to investigate the isomerization of CA's model substrate HAGPIA in water and in its target protein CypA. Our results allow us to propose a novel mechanistic hypothesis, which is finally consistent with all of the available molecular biology data.

The best characterized isomerization in vivo and in vitro occurs in the uncoating and recruitment of the human HIV-1 capsid (CA) protein in the virions [7], and it is catalyzed by cyclophilin A isomerase (CypA, [12]). At the structural level, CypA features ahelices flanking a beta-barrel, while CA is made of several ahelices connected by loops ( Figure 1). In the X-ray structure of the complex [13], the targeted proline-containing backbone unit (G89-P90), is accommodated in a hydrophobic pocket of CypA (residues F60, F113, L122, and H126 in Figure 1). Although the backbone unit is in trans conformation in the X-ray structure [13], NMR studies have shown that 45% of conformers are cis for CA-CypA R55A in aqueous solution [12]. The rather significant population of cis conformers could arise not only by the replacement of R55 with Ala, but also by crystal packing forces, different temperature conditions (100 K and 298 K for the X-ray and NMR experiments, respectively), along with different hydration and ionic strength in the two experiments. Free energy calculations further support the hypothesis that CypA promotes a significant population of cis conformation [14,15]. Thus, the cis population is likely to increase substantially from water -where it is ,10% [12] -to the complex in aqueous solution.
In vitro kinetic measurements (Table 1) show that CypA decreases the isomerization free energy barrier (DG { cisRtrans , DF { cisRtrans ) of modified substrate fragments in solution by ,7 kcal/mol [16,17]. Free energy calculations, based either on classical [15] or quantum-mechanical/molecular mechanics (QM/MM) [14], of the 6 aminoacids long substrate fragment point to a similar trend (Table 1). In addition, they suggest that Hbonding between the backbone unit of the targeted glycine-proline peptidyl bond and that of N102@CypA [14,15,18], as well as van der Waals interactions between substrate and the CypA hydrophobic pocket (F60, F113, L122, and H126) stabilize the transition state (TS) [14,15,18]. These hypotheses are consistent with the decrease of k cat /K m associated with the polymorphism of cyclophilins in 102 position (N to T, S, H and R, Table S1) and in the F60A, F113A and H126Q CypA mutants [19,20] (Table S1). NMR studies, along with computations, have further suggested that the motions of several CypA residues (including R55 and N102) are linked to the enzymatic activity [15,[21][22][23]. Free energy calculations characterized a network of protein vibrations in CypA that are associated with its isomerase activity: flexible loops on the surface are connected to the active site by a network of hydrogen bonds [15,22,23]. In spite of the great progress in describing the catalytic process, key mechanistic issues need to be addressed. Kinetic studies [17] suggested that the overall process involves trans and cis forms in solution and in complex with the protein. However, the studies so far consider mostly TS stabilization. Most importantly, the current proposed mechanism cannot explain a plethora of molecular biology data. These studies show that residues not involved in TS stabilization in the proposed mechanisms are important for the function, as their mutations cause a decrease of k cat /K m . Indeed, k cat /K m passes from 1.6 10 7 M 21 s 21 in the wild type to 1.6 10 4 M 21 s 21 by mutating the fully conserved R55 residue with Ala [20]. Although this mutation was proposed originally to destroy an H-bond stabilizing uniquely the TS [13,24], such Hbond was subsequently ruled out in computational works [14,15,18] and so far no functional role has been ascribed to R55. In addition, the H54Q mutation, along with the I57V polymorphism, causes a decrease of k cat /K m , although these residues are not involved in TS stabilization (Table S1).
Here we use molecular simulations to address these issues. We calculate the free energy associated with the isomerisation of the 6 aminoacids long (/HAGPIA/) substrate fragment in water and in complex with CypA. The free energy is calculated as a function of the two reaction coordinates f and y (defined in Figure 2), which have been suggested to describe best the energetics of the process [25], as well as other pairs of different coordinates to cross-check our results. We use here metadynamics [26] (in its bias exchange extension [27]), which has already been employed to predict the energetics of protein/peptide interactions [28]. The potential used is an effective one (specifically the AMBER99 force field [29]). This choice allows a very efficient sampling because of its relatively small computational cost. In spite of its limitations, [30] this force field is expected to be relatively accurate to describe equilibrium conformations, which is a key aspect of our problem. The accuracy of the force field for minima and transition states is assessed by comparing our results with first principle quantum chemistry calculations. Based on this comparison, we find that the force-field biases on the energetics of the minima are negligible, whereas their influence on the barriers is more significant. Therefore, here the free energy differences between minima are used to predict the relative populations of the equilibrium conformations, whilst the calculated barriers are only compared at the qualitative level to discriminate the most likely cis/trans isomerisation path.
The enzyme turns out to dramatically stabilize the populations of one specific cis conformer relative to the trans ones, which instead are by far the most stable in aqueous solution. In addition, it lowers the free energy barrier of a specific, one-way isomerization from trans to cis. These findings allow us to propose a mechanistic hypothesis for the isomerization process that is consistent with all the available experimental data.

Structural Models and MD Calculations
The initial structural models of /HAGPIA/ in the free state (PEPT-WAT) and in complex with CypA (PEPT-CypA) were obtained from the CA N-terminal domain/CypA X-ray structure at 2.0 Å resolution (PDBID:1M9C) [13]. The X-ray data were collected at 100 K and slightly basic pH (8.0) [13]. The peptide and the protein were considered in their zwitterionic form. The protonation states of all histidine residues were determined by pKa calculations based on the H++ server [31,32] and by visual inspection of the hydrogen bond network. All His residues were protonated at Nd1. The result for PEPT-WAT is consistent with the NMR study on CA N-terminal domain at pH 7 [33].
PEPT-WAT and PEPT-CypA were solvated in a 32 Å637 Å638 Å and 57 Å680 Å 60 Å boxes, containing 1,423 water molecules (for a total of 4,349 atoms) and 8,233 water molecules and a Cl 2 ion (for a total of 27,283 atoms), respectively. The ion was added to neutralize the system.
The Amber99 all atom force field [29] was used for the protein and the chlorine atom; the TIP3P model [34] was used for water. After minimization up to a convergence of 10 24 kcal/mol (conjugate gradient algorithm), the system was equilibrated for 1 ns (time step of 1 fs) in an isothermal-isobaric ensemble (1 atm, 310 K) with the Langevin barostat [35] (the oscillation period and the decay coefficient were set to 200 fs and 100 fs, respectively) and thermostat [36] (coupling coefficient 1 ps 21 ). We use the Particle Mesh Ewald scheme [37] with 12 Å cutoff and 0.75 Å -spaced Fourier grid; we assume a dielectric constant of 1. Van der Walls interaction cutoff was set to 12 Å . Minimization and equilibration were performed with the NAMD 2.6 program suite [38].

Free Energy Calculations of Proline cis/trans Isomerisation
Free energies were calculated using the metadynamics method in its bias exchange variant [26,27]. This approach consists in a continuous addition of a history dependent potential energy that enforces the dynamics to explore conformations that were not previously visited. Briefly, the forces acting on each atom, or centroid, are corrected by a history dependent contribution, obtained as the derivative of the history dependent potential energy V t with respect to the atom coordinate. V t is given by the sum of a set of Gaussian functions centered on the values s t of each chosen collective variable (CV) s as time t: The time interval between the addition of two Gaussian functions t, as well as the Gaussian height w and Gaussian width d, are tuned to optimize the ratio between accuracy and computational cost (Table S2). Eventually, after exploring all conformations defined within the CV space, the probability distribution of Gaussians becomes flat, and the free energy profile does not change any more. At this stage, the free energy surface can be easily reconstructed as the opposite of the sum of all Gaussians.
Here we use the bias exchange variant of this method [27], in which several trajectories x a , x b with different history dependent potentials (for instance two CVs a and b), V a G , V b G , are run in parallel. At specific time intervals (of the order of few ps), x a and x b are swapped with a probability P ab evaluated by the standard Metropolis scheme [39].

P ab~1
if Dƒ0 Author Summary Peptidyl prolyl isomerases are ubiquitous enzymes whose actions are crucial in several biological processes, such as, for instance, in cellular signalling and in the onset of several diseases, e.g., HIV infection. Therefore, these isomerases are promising targets for the design of new drugs. For this purpose, we need to understand their molecular mechanism of action. One of the most characterized peptidyl prolyl isomerases is cyclophilin A. Previous studies characterized the roles of several protein regions in isomerase function. However, there are still experimentally identified important portions of the protein whose specific actions in the mechanism are still not known. Here, we address this problem by an extensive computational study of cyclophilin A and a substrate peptide that is part of the HIV-1 capside protein. We present a novel four-step mechanism of the whole enzymatic process, which is consistent with all of the available experimental data. Moreover, these steps can be used as targets for the development of drugs, e.g., for HIV-1 infection. The X-ray structure shows that it is a b-barrel with eight antiparallel b-strands and two a-helices flanking the b-barrel [13]. These secondary structures are connected by several flexible loops. The b-barrel contains the active site for cis/trans prolyl isomerization. The CA protein is composed by two helices connected by a loop. This contains the G89-P90 target peptide bond (Bottom, right). Bottom, left: Transition state of the cis-trans isomerisation, emerging from this and previous [14,18,22]  where k B is the Boltzmann constant and T is the temperature. This method was shown to provide an efficient sampling of the conformational space [27]. The free energies associated with cis/trans isomerization of the G3P4 bond in PEPT-WAT and PEPT-CypA were calculated as a function of pairs of CVs in the canonical ensemble. The free energies were first calculated as a function of the f and y dihedral angles ( Figure 2 and Table S2), because these angles have been shown to be essential to describe the isomerisation process [25]. The resulting free energy plot was used to define the minima and  [17], T = 273 K for [16]).  the transition state regions (R hereafter). Briefly, we performed a cluster analysis of all the structures [40] and evaluated their population using an umbrella sampling-like approach [41]; minima and TSs are characterized by several clusters and only one cluster, respectively (See Text S1 for details). Free energy barriers were then calculated as differences between the TS cluster and the lowest cluster in the basin of the minimum.
The sampling of TSs was higher in PEPT-WAT than in PEPT-CypA, as shown by the few number of structures clusterized in TSs' R (i.e. less than 10 structures). To improve sampling in PEPT-CypA, we performed additional 4 ns-long MD simulations in which the f and y angles were harmonically restrained at the values of the TSs. Restraint center and their associated force constants were fitted so as to keep the entire MD within the region of interest.
To test the robustness of our insights on the enzyme mechanism obtained by the above free energy profiles, we compared our results with free energy plots as a function of other pairs of CVs. These are: (i) the f angle and the proline nitrogen (P4N) pyramidalization p that might have a role in peptidyl prolyl cis/ trans rotation [25]. p is defined as the distance between P4N and the centre of a plane determined by three atoms belonging to G3 and P4 ( Figure 3 and Table 2). (ii) The f and P4N H-bond coordination number (s1) that has been shown to be important for in vacuo prolyl isomerization [25]: where r P4N-j is the distance between P4N and the atom j (excluding solvent molecules) while r 0 , n, m are parameters chosen to obtain s1,1, when the P4N forms one H-bond (Table S2). s1 increases with the number of H-bonds formed by P4N with H-bond acceptors of the peptide (in water) or of the peptide and protein in the complex.
(iii) p and s1, to assess the correlation between the pyramidalization and H-bonding of the P4N atom.
Finally, we calculated the free energy profiles as a function of selected pairs specific for PEPT-WAT or for PEPT-CypA. For PEPT-WAT, these are the following: (i) the f and P4N coordination number with water solvent molecules (s2). s2 is defined as s1 except that j runs over the water oxygen atoms. This pair tests whether water is able to catalyze prolyl isomerization by forming an H-bond with P4N. (ii) p and s2, to evaluate if there is a correlation between P4N pyramidalization and its H-bonding to water molecules.
For PEPT-CypA, we introduced CVs that relate the substrate to the enzyme: i) hydrophobic coordination numbers s3, s4 and s5, i.e. quantities which depend only on non-polar carbon atoms. s3, s4 and s5 describe the hydrophobic interaction between the peptide and the protein: where r ij is the distance between atoms i and j. i runs over the residues involved in the rotating prolyl bond of the peptide (G3P4) for s3, and the C and N terminal of the peptide (H1, A2, I5 and A6) for s4; j runs over any non-polar carbon of CypA belonging to residues located at 4 Å or less from the peptide (Table S2) for both s3 and s4. The parameters n and m were chosen so as to distinguish conformations with and without hydrophobic interactions (Table   S2). This was done by running 1 ns MD test simulations with different n and m parameters. s5 describes a subset of such hydrophobic interactions, which play a particularly important role during the isomerisation, as suggested by NMR studies [21]. These are the interactions of L98 and S99 in CypA and C terminal of the peptide (I5, A6). Thus, for s5, i runs over I5 and A6, while j runs over the atoms of L98 and S99.
Finally, we introduced the reaction coordinate s6: where j runs over all the H-bond donors/acceptors of the peptide or of CypA that were found within 3Å from R55 N atoms in 1 ns MD: those turned out to belong to residues Q63 and N149; r R55Ns-j is the distance from any N atom of R55 and any atom j; r 0 , n and m are chosen as in s1 (Table S2). This reaction coordinate increases with the number of H-bonds formed with R55 and hence measures the ability of this residue to form H-bonds with the peptide or CypA residues (organizing its active site): this has been indicated as crucial for the catalysis [14]. The free energy is calculated as a function of (i) f, s3; (ii) f, s4; (iii) f, s5 and (iv) f, s6 to analyze hydrophobic (s3, s4, s5) and hydrophilic (s6) interactions along with the cis (f,0u)«trans(f,6180u) interconversion.
All metadynamics calculations (17 ns for PEPT-WAT and 40 ns for PEPT-CypA) were performed at the physiological temperature of 310 K. The temperature was controlled by a Nosè-Hoover thermostat [42], with coupling time constant of 0.05 ps. Electrostatic interactions were assessed using the particle mesh Ewald schemes [37] with 34 wave vectors in each dimension and fourth-order cubic interpolation. Van der Waals interactions were evaluated as specified for the equilibration phase. The time-step was 1 fs. The Gaussians were added with a frequency of 2 GHz and they had a height of 2.5 kJ/mol, the width (d) for each CV is reported in Table S2. The exchange trial frequency was 400 MHz. All bonds were constrained with the LINCS algorithm [43]. We removed the protein centre of mass translation every 10 steps of molecular dynamics. The calculations were performed with a locally modified version of Gromacs 3.3.1 [27,44,45].

Calculated Properties at Minima and TSs
The number of hydrogen bonds (N HB ) was calculated in TS and minima R, supposing a Boltzmann distribution within the cluster, as assumed in [46] for enthalpy calculations: where R is the region identifying the minima or the TS (See Figure  S1), N i is the number of structures belonging to cluster i (see Figure  S1), N HB (i) is the number of H-bonds and F(i) is the free energy for cluster i. The H-bond between the acceptor P4N and donors D is assumed to exist if the distance (D-P4N)#3.2 Å and the angle (P4N…D-H)#70u.
The population of puckered conformations P x2(C) is defined as: where C = (up, down or planar puckering) and i is a cluster belonging to this state. The free energy associated with each conformation C is estimated as F C~{ where i is the index of the carbon atoms and j runs over the carbon atoms of CypA hydrophobic residues within 4 Å of G3 and P4 residues.
The internal energy/entropy contributions are evaluated for the restrained MD simulations of each minima and TS. The internal energy contribution to the free energy is calculated as the potential energy averaged over all conformations of the cluster [46], while the entropic contribution is obtained by the standard thermodynamic relation: TDS~DU{DF .
Principal component analysis (PCA) [48,49] was used to identify large-scale collective fluctuations in minima and TSs restrained dynamics. This analysis was performed on Ca of PEPT-CypA complex, using Gromacs 3.3.1 [44,45]

Errors
The statistical error of the free energy for minima and TSs are estimated as the largest difference between the F(f) values calculated from different 2D free energy profiles (for more details see text S1). The statistical error on the enthalpy of minima and TSs was estimated similarly. The accuracy of the force-field used (Amber [29]) was established by a comparison with quantum chemical results on a model system. We compared the cis«trans isomerization potential energy of a N-acetyl proline methylamide using the Amber force field [29] and DFT calculations at the B3LYP/6-31G(d) [50][51][52][53] level of theory (see text S1 for more details on these calculations). Classical calculations were performed with Gromacs 3.3.1 [44,45]. DFT calculations were performed with Gaussian98 [54].

Results
We calculated the free energy F associated with the isomerization of the G3P4 bond in the /HAGPIA/ peptide as a function of the angles f and y (defined in Figure 2), which are fundamental collective variables (CVs) to describe the isomerization process [25].
The peptide is either solvated in water (PEPT-WAT) or in complex with CypA, in aqueous solution (PEPT-CypA). Forcefield based metadynamics in its bias-exchange variant, in the canonical ensemble, provided the F = F(f, y) profile, from which we defined minima and TS regions and clusters within these regions (See Methods Section).

PEPT-WAT
The free energy profile plot shows that the trans 0 conformer, in which f,6180u and y,0u, is the absolute minimum ( Figure 4 and Table 2). The second trans minimum is trans 180 (f,6180u and y,6180u), whose free energy is higher by 3 kcal/mol (Table 2,    (Table S3, see Figure 2 and Figure 1 for atom labelling). However, the number of waters around the peptide is larger for trans 180 than for trans 0 , since the former has a more extended structure ( Figure 4). This water-peptide interaction may not stabilize enough trans 180 against trans 0 to overcome the stabilization of the latter by the P4N…I5N intramolecular H-bond.
The free energies of the two cis conformers, (cis 0 : f,0u and y,0u and cis 180 : f,0u and y,180u) are 3 and 5 kcal/mol higher than trans 0 , respectively. The stabilization of trans relative to cis has been ascribed to steric clashes in the trans conformations (see, e.g. [10]), and is affected by intramolecular interactions. In our case, in trans A6 H-binds to H1, A2 and G3 with higher persistence than cis (Table   S3). This is expected to stabilize further the trans conformation. The reason why cis 0 is more stable than cis 180 is again the presence of the P4N…I5N H-bond only in the first conformation.
As in the case of trans conformers, there are more waters around the peptide in the cis 180 than in the cis 0 conformation, but waterpeptide interactions do not stabilize significantly the cis 180 conformation compared to the intramolecular H-bond interaction that stabilize cis 0 .
The difference in hydration between trans 180 and cis 180 is reproduced between trans 0 and cis 0 . A direct comparison of our calculation with classical and QM/MM umbrella sampling studies based on the reaction coordinates v [15] and t [14], respectively, is not possible, because these calculations do not distinguish between cis and trans conformations at y = 0u and y = 180u.
The trans X Rcis X and cis X Rtrans X (X = 0u,180u) isomerization free energy barriers range between 4-18 and 11-15 kcal/mol, respectively (Table 3). These values are similar to those obtained by force-field based umbrella sampling calculations performed on the v variable (Table 1) [15]. However, early quantum gas phase calculations showed that, at variance with the f angle, the v angle alone does not describe properly the proximity to a saddle point conformation [25].
The water shell does not change along the trans X «TS«cis X pathways ( Figure 4). This is consistent with the suggested poor solvent reorganization during isomerisation in peptides of similar size [55,56]. However, the solvent could play an important role for the isomerization of peptides smaller than that considered here.  The isomerisation is enthalpy-driven ( Figure S2) similarly to that found experimentally on different proline containing peptides [52,53].
Furthermore, proline puckering populations are fully consistent with statistical distributions across peptidyl prolyl bonds in the Brookhaven Protein Databank [47] (see Text S1 and Figure S3).
The lowest pathway is the cis 0 RTS1Rtrans 0 pathway with an associated barrier of 11 kcal/mol (Figure 4), also shown by in vacuum DFT potential energy calculations of N-acetylproline methylamide isomerization (see Text S1 and Table S4). This value is of the same order of the experimental values for smaller peptides (,4 residues, Table 1). The stabilization of TS1 (as opposed to TS2-TS4) may be caused, at least in part, by the larger persistency of the P4N-I5N H-bond (Table S3 and Table S5). The formation of this H-bond, not only for TS1 but also for TS3, has been previously predicted [25]. P4N instead does not interact with the solvent (See Text S1).
Obviously, the cis 0 Rtrans 0 pathway is favored over the reverse one, since trans 0 is the most stable state.
The isomerizations involving states with y,6180u are higher in free energy because the P4N-I5N intramolecular H-bond is not present in these conformation ( Figure 4).
Again, a comparison with previous free energy calculations in aqueous solution is not possible since previous studies did not discriminate between trans/cis 0 and trans/cis 180 .
A similar picture emerges from calculations of the free energy as a function of other CV pairs ( Table 2). These include the P4N pyramidalization, the key H-bond between P4N and the H-bond donors of the peptide (Figure 3), and P4N hydration (See Methods, Text S1 and Figure S4).
We conclude that trans 0 is the most populated state (Figure 4 and Table 2) and that the fastest kinetic process produces this isomer, starting from the most populated cis isomer (cis 0 ).

PEPT-CypA
The presence of the protein alters dramatically the population of the four minima in the F = F(f, y) plot ( Figure 5). The global minimum, and therefore the most populated state, is not a trans configuration: it is cis 180 (representative structure in dataset S4). Cis stabilization has been already found by AMBER-based and QM/ MM free energy calculations as a function of v [15] and t [14] (t dihedral angle -defined as C (i-1) -O (i-1) -C d(i) -C a(i-1) -is similar to f -C a (i-1) -O (i-1) -C (i) d -C a (i-1) -used in this work). Trans 0 (representative structure in dataset S1) and trans 180 (representative structure in Dataset S2) are 1 kcal/mol higher in free energy than cis 180 ; cis 0 (representative structure in Dataset S3) is scarcely populated ( Figure 5 and Table 2).
N102 stabilizes cis 0 , as cis 180 , forming an H-bond to P4O. Several residues of CypA form hydrophobic interactions with cis 0 . Most of these residues stabilize also other minima (Q63, A101, F113). However, H54@CypA and A103@CypA stabilize only cis 0 , the minimum with the highest free energy. The lowest free energy barrier is associated with the trans 180 RTS2Rcis 180 (counterclockwise) path (Tables 2 and 3). This is also the only isomerization pathway catalyzed by the enzyme relative to PEPT-WAT (Tables 2 and 3). Preferential lowering of the counterclockwise N-terminal rotational free energy barrier is consistent with previous classic MD and free energy calculations [15,57,58].
The enzyme does not decrease the barrier of the reverse pathway (cis 180 RTS2Rtrans 180 ), as CypA stabilizes similarly cis 180 and TS2 (the latter is a bit more stabilized given the higher persistence of hydrophobic and hydrophilic interactions).
In addition, the H-bonds and hydrophobic contacts stabilize TS3 (representative structure in Dataset S7) and TS4 (representative structure in Dataset S8) along with their connected minima (Table 4, Table 2, and Figure 5). We further notice that the Hbonds in TS3 are the same as in TS2, although less persistent (Table 4). TS1 (representative structure in Dataset S5) is not stabilized by H-bonds and interacts weakly with CypA hydrophobic residues (in particular, F60, L122 and W121, Table 4). Therefore the peptide in this conformation does not form a tight complex with the protein and it is exposed to the solvent ( Figure 5). Thus, this isomerization pathway (trans 0 «TS1«cis 0 ) is not too dissimilar from that in water and, indeed, the barriers for the two processes are practically identical ( Table 2).
In each minimum and TS we identify large collective motions of the PEPT-CypA complex using PCA. We observe that significant modes involve almost the same residues as the motions found in Ref. [15](see Figure S6).
As for the peptide in water, we used other CV pairs (in addition to f, y) to describe the cis/trans isomerization. These include the pyramidalyzation as well as the number of H-bonds formed by P4N and residues of the peptide or of CypA. In addition, we included coordinates to take into account (i) the hydrophobic interaction between the peptide and the enzyme, and (ii) the interactions of R55 with both the peptide and the CypA active site (for a discussion of the other CV pairs see Figure S5, Text S1, and Table S7). All these calculations provide a consistent picture that leads us to conclude that the enzyme catalyzes only one pathway, the one from the most populated trans conformation, trans 180 , to the most stable minimum, cis 180 . Table 4. Number of H-bonds of CypA residues with G3P4@PEPT along prolyl cis/trans isomerization.  Table 5. Interface coordination number a of CypA residues with G3P4@PEPT along prolyl cis/trans isomerization.

Discussion
Our calculations produced the energetics of minima and TS of PEPT-WAT and PEPT-CypA. A critical discussion of the accuracy of the calculations is necessary to proceed with possible insights on the enzyme working mechanism. The source of errors in a free energy calculations are basically of two types [59].
The first is caused by a non-complete sampling (statistical error): our estimated statistical error turns out to be relatively small, less than 1.5 kcal/mol (see Text S1).
The second is caused by the force field. This is expected to be particularly relevant for the TSs, because P4N passes from almost sp 2 hybridization (resonance of prolyl-peptidyl bond) in minima to sp 3 at the TSs (the so-called pyramidalization). Notice that since the Amber99 force-field is not designed for having all bonds fixed [29], the isomerization barriers might have been overestimated. However, we expect this effect to bias all barriers in the same direction, so that our conclusions on free-energy differences should not be affected by this choice.
To estimate the error due to the force field, quantum chemical calculations are compared with force field based ones on model systems. Previous work has shown that DFT-B3LYP with 6-31G(d) basis set provide similar results as more expensive (albeit usually more accurate) MP2 calculations on the N-acetyl-N-N9methyl proline amide [60]. Therefore, here we limit ourselves to a DFT-B3LYP/6-31G(d) calculation on the N-acetyl proline methylamide system, in vacuum, at 0 K, and compare our results with AMBER calculations.
We find small differences for the minima (of the same magnitude of the statistical error, 1.5 kcal/mol, see Text S1). Instead, the differences are 4 kcal/mol or lower for the barriers. We further notice that in three of the reversible pathways (cis 0 «TS1«trans 0 ; cis 180 «TS2«trans 180 ; cis 180 «TS4«trans 180 ), DFT barriers are higher than force field ones, while QM barriers of cis 0 «TS3«trans 0 are slightly lower than the force-field ones (Table S4).
We conclude that our calculations may provide a relatively accurate estimate of the Boltzmann populations of cis and trans free energy minima. In addition, as the errors associated to the barriers are likely to be similar for all pathways, our calculations can be used for qualitative comparisons of the free energy barriers, i.e. to identify the lowest free energy pathways.

Population of Free Energy Minima
In water, trans 0 is by far the most populated specie ( Table 2). The same result for the potential energy has been obtained with DFT (see Text S1) gas phase calculations on N-acetylproline methylamide.
The protein environment changes dramatically the populations of the proline conformers, stabilizing cis with respect to trans, as previously reported by CypA-SUC-Ala-Phe-Pro-Phe-pNA complex [17] and CypA R55A -CA [12]. However, our study provides a quantitative estimate of the populations of conformers (Table 2): in presence of the enzyme, three conformers (trans 0 , trans 180 and cis 180 ) are significantly populated. The most populated one corresponds to cis 180 , that is exactly the most unfavorable state in aqueous environment. This finding has never been reported in literature.

Qualitative Comparison of Free Energy Barriers
CypA accelerates only one interconversion, namely the one from trans 180 to TS2 to cis 180 (Table 3), as already suggested by previous calculations [15]): our calculated free energy barrier decreases by 4 kcal/mol with respect to the peptide in water.
Indeed, as well as in the other pathways, almost no effect is found on the reverse pathway (cis 180 RTS2Rtrans 180 , 1 kcal/mol, i.e. within the error of the calculations), again in agreement with Refs. [14,15] (Table 3). The trans 180 RTS2Rcis 180 is also associated with the lowest free energy barrier. We conclude that this pathway is the most likely in the enzyme, although we cannot establish with high accuracy the free energy barrier of the enzymatic reaction.
Based on our calculation, we propose the following mechanism ( Figure 6): 1. The most abundant conformation of the peptide in water is, as it is well known [12], trans conformation. Specifically, in our system it is trans 0 . 2. CypA sequesters trans 0 and rapidly interconverts it into trans 180 , because the population of the latter is significantly greater in the protein complex. Notice that the process is much faster than the conversion to cis 180 or to cis 0 , as the activation free energy associated with them is one order of magnitude smaller. 3. The enzyme catalyzes the isomerization to the most populated minimum, cis 180 , (trans 180 RTS2Rcis 180 pathway). This minimum is strongly stabilized by the same PEPT-CypA interactions as those observed in the most stabilized transition state (TS2), but less persistently ( Table 4). Notice that the reverse interconversion from cis 180 to trans 180 is not likely because: (i) the free energy barrier associated with it is larger (indeed, the trans 180 RTS2Rcis 180 pathway is the only pathway catalyzed by the enzyme) and (ii) the population of trans 180 is a few times smaller than cis 180 (Table 2). This contrasts with isomerization in water solution, where the population of the most stable trans conformation (trans 0 ) is several orders of magnitude higher than the one of the preponderant cis (cis 0 ). 4. The enzyme can now either detach cis 180 or interconvert it to cis 0 . However, only in the latter case H54 does interact with the peptide, consistently with the experimental fact that the H54Q mutation decreases the enzymatic activity; the k cat /K m is reduced to the 15% of its wild-type value [20]. Notice that all the enzymatic activity measures, for mutants of polymorphism, do not provide information on k cat or K m alone but only on their ratio. Therefore one cannot conclude whether a mutation affects either k cat or K m or both. Thus, it is plausible to assume that the enzyme interconverts cis 180 to cis 0 . This assumption may be validated against further molecular biology experiments: in particular, it would be highly useful to measure the activity of A103X mutants (X = polar or charged residues), because A103 stabilizes only cis 0 . 5. cis 0 is released.
Our mechanism provides a first rationale for mutational data which have not been explained so far. The decreased k cat /K m values (Table S1) of the W121A mutant [20] and the I57V variant [19] cannot be explained in terms of loss of TS stabilization. In our mechanism, these two mutations are likely to affect substrate stabilization in the trans 0 conformation, reducing CypA ability to capture the substrate. In addition, the dramatic decrease of enzyme efficiency in the R55A mutant (0.10% residual activity [20]), whose effect on TS stabilization is controversial, may be, at least in part, a consequence of the reduced population of trans 180 , the reactant of the CypA catalyzed pathway (trans 180 R TS2Rcis 180 ).
Next, we analyze the effect of several mutations and polymorphism, so far ascribed only to TS stabilization [14,15,18], as also found here, that might also have an impact on ground state populations.
Indeed, reduced k cat /K m values on enzymes where N102 exchanges to T, H and R residues [19] may be not only due to destabilization of the TS, but also of the cis 180 conformation. Similarly, some CypA mutations also decrease k cat /K m : (i) F60A [20] may destabilize TS2, trans 0 and cis 180 ; (ii) F113A and H126Q [20] may affect all the species along the trans 180 RTS2Rcis 180 pathway (Table S1).
cis 0 , proposed to be the most probable final step in our mechanism based on H54Q mutant, is stabilized exclusively by interactions of H54 and A103 with G3P4. Figure 6. Proposed mechanism of action for CypA. CypA sequesters the most abundant conformation in water, trans 0 , that is rapidly interconverted into the most abundant conformer trans 180 . Then, CypA catalyzes the isomerization of the peptide along TS2, producing the mostly populated minimum, cis 180 . The peptide most probably detaches in the cis 180 conformation. This readily interconverts to cis 0 once the peptide is in aqueous solution. CypA residues that form important H-bonds (hhb) and hydrophobic interactions (nnb) with the G3P4 moiety are shown. Residues with almost exclusive relevance to each conformation are highlighted in red. doi:10.1371/journal.pcbi.1000309.g006 We conclude that our mechanism is consistent with all the mutational and polymorphism data and provides a structural basis for most of them.

Supporting Information
Dataset S1 Representative structure of trans0 conformation (center of the lower free energy cluster within trans0 region) in PDB format.  Figure S2 Non catalyzed /HAGPIA/ peptide cis/trans prolyl isomerization: free energy (kcal/mol) and enthalpy (kcal/mol) as a function of the dihedral angle f (in degrees) obtained by integrating out the variable p. The entropic contribution is given by the difference between the two curves. The isomerization process is mainly enthalpy-driven.
Found at: doi:10.1371/journal.pcbi.1000309.s010 (0.15 MB TIF) Figure S3 (Top) Up, planar and down proline ring puckering conformations. Puckering is defined based on the value of the x2 dihedral angle (Ca-C b-C c-C d): up-puckering for x2.10u, planar puckering for 210u,x2,10u, down-puckering for x2,210u. (Bottom) x2 distribution within clusters classified as up, down and planar puckering. x2 has a bimodal distribution in all clusters with a largest maximum at x2 = +40u for up clusters, at x2 = 240u for down clusters and with two even peaks at 240u and +40u for planar clusters. Found at: doi:10.1371/journal.pcbi.1000309.s011 (0.27 MB TIF) Figure S4 (A) Non catalyzed /HAGPIA/ peptide cis/trans prolyl isomerization: free energy (kcal/mol) as a function of the dihedral angles f (in degrees, showed in the inset) and the pyramidalization p. (B) Non catalyzed /HAGPIA/ peptide cis/ trans prolyl isomerization: free energy (kcal/mol) as a function of the dihedral angles f (in degrees, showed in the inset) and the coordination of P4N with peptide H-bond donors (s1). (C) Non catalyzed /HAGPIA/ peptide cis/trans prolyl isomerization: free energy (kcal/mol) as a function of the dihedral angles f (in degrees, showed in the inset) and the coordination of P4N with all water molecules (s2). (D) Non catalyzed /HAGPIA/ peptide cis/trans prolyl isomerization: free energy (kcal/mol) as a function of the pyramidalization (p) and the coordination of P4N with peptide Hbond donors (s1). (E) Non catalyzed /HAGPIA/ peptide cis/trans prolyl isomerization: free energy (kcal/mol) as a function of the pyramidalization (p) and the coordination of P4N with all water molecules (s2). Found at: doi:10.1371/journal.pcbi.1000309.s012 (2.08 MB TIF) Figure S5 (A) Enzyme catalyzed /HAGPIA/ peptide cis/trans prolyl isomerization: free energy (kcal/mol) as a function of the dihedral angles f and the pyramidalization p. (B) Enzyme catalyzed /HAGPIA/ peptide cis/trans prolyl isomerization: free energy (kcal/mol) as a function of the dihedral angles f and the coordination of P4N with peptide and enzyme H-bond donors (s1). (C) Enzyme catalyzed /HAGPIA/ peptide cis/trans prolyl isomerization: free energy (kcal/mol) as a function of the dihedral angles f and the hydrophobic coordination of G3P4 with CypA (s3). (D) Enzyme catalyzed /HAGPIA/ peptide cis/trans prolyl isomerization: free energy (kcal/mol) as a function of the dihedral angles f and the hydrophobic coordination of pepide C-and Ntermini (H1, A2, I5, A6) with CypA (s4). (E) Enzyme catalyzed / HAGPIA/ peptide cis/trans prolyl isomerization: free energy (kcal/mol) as a function of the dihedral angles f and the hydrophobic coordination of peptide C-term (I5A6) with L89 and S90 (s5). (F) Enzyme catalyzed /HAGPIA/ peptide cis/trans prolyl isomerization: free energy (kcal/mol) as a function of the dihedral angles f and the coordination of R55 with P4N and protein H-bond donors (s6). (G) Enzyme catalyzed /HAGPIA/ peptide cis/trans prolyl isomerization: free energy (kcal/mol) as a function of the pyramidalization p and the coordination of P4N with peptide H-bond donors (s1).       Text S1 Clustering, error evaluation, puckering of peptide in water and analysis of secondary free energy profiles obtained in the bias exchange metadynamics simulations. Found at: doi:10.1371/journal.pcbi.1000309.s022 (0.10 MB DOC)