Free-Energy Landscape and Rate Estimation of the Aromatic Ring Flips in Basic Pancreatic Trypsin Inhibitors Using Metadynamics

Aromatic side chains (phenylalanine and tyrosine) of a protein flip by 180° around the Cβ–Cγ axis (χ2 dihedral of the side chain), producing two symmetry-equivalent states. The study of ring flip dynamics with nuclear magnetic resonance (NMR) experiments helps to understand local conformational fluctuations. Ring flips are categorized as slow (milliseconds and onward) or fast (nanoseconds to near milliseconds) based on timescales accessible to NMR experiments. In this study, we investigated the ability of the infrequent metadynamics approach to estimate the flip rate and discriminate between slow and fast ring flips for eight individual aromatic side chains (F4, Y10, Y21, F22, Y23, F33, Y35, and F45) of the basic pancreatic trypsin inhibitor. Well-tempered metadynamics simulations were performed to estimate the ring-flipping free-energy surfaces for all eight aromatic residues. The results indicate that χ2 as a standalone collective variable (CV) is not sufficient to obtain computationally consistent results. Inclusion of a complementary CV, such as χ1(Cα–Cβ), solved the problem for most residues and enabled us to classify fast and slow ring flips. This indicates the importance of librational motions in ring flips. Multiple pathways and mechanisms were observed for residues F4, Y10, and F22. Recrossing events were observed for residues F22 and F33, indicating a possible role of friction effects in ring flipping. The results demonstrate the successful application of infrequent metadynamics to estimate ring flip rates and identify certain limitations of the approach.


Circular Correlation of Side-chain torsions
The circular correlation coefficient 1 for dihedrals is calculated using astropy python package 2 .Here, we report the results of infrequent metadynamics simulations starting from stable states A, I1 and directed symmetry-equivalent A' and I1' states of residue F22.The trajectories are directed to specific intermediates using harmonic wall bias (see Equation 3 in main manuscript).The details of wall bias parameters are provided below: The transition rate matrix based on these simulations is reported in Table S5.The diagonal values [3][4] of this matrix are set to negative of the sum of row elements.All eigenvalues are negative real values as suggested in reference 3. The dominant eigenvalue 83.0  −1 was observed.

Figure S3 .
Figure S3.The time series of 2 and 1 during  2 −  1 metadynamics simulation showed in figures (a) and (b), respectively.(c) the projection of free energy on  2 in time intervals of 100 ns.(d) The free energy difference between two symmetric states along  2 .

Figure S4 .
Figure S4.The diffusion of  2 and  1 variables during  2 −  1 metadynamics simulation of residue F33 (figures (a) and (b), respectively).(c) the projection of free energy on  2 in time intervals of 100 ns.(d) The free energy difference between two symmetric states along  2 .

Figure S5 .S6Figure S6 .
Figure S5.The diffusion of  2 and  1 collective variables during 2 metadynamics simulation of residue Y21 (figures (a) and (b), respectively).(c) the free energy along  2 in the time intervals of 100 ns.(d) The free energy difference between two symmetric states along  2 .

Figure S7 .
Figure S7.The time series of  2 and contact map (CMAP) collective variables during  2 -CMAP metadynamics simulation of residue Y23 in figures (a) and (b), respectively.(c) the projection of free energy along  2 in time intervals of 100 ns.(d) the free energy difference between two symmetric states along  2 .

Figure S8 .
Figure S8.The time series of  2 and  1 during  2 −  1 metadynamics simulation of residue Y35 in figures (a) and (b), respectively.(c) the projection of free energy on  2 in time intervals of 200 ns.(d) The free energy difference between two symmetric states along  2 .

Figure S9 .
Figure S9.The time series of  2 and  1 during  2 metadynamics simulation of residue Y35 in figures (a) and (b), respectively.

Figure S10 .
Figure S10.The diffusion of  2 and  1 during  2 −  1 metadynamics simulation of residue F45 in figures (a) and (b), respectively.(c) the projection of free energy on  2 in time intervals of 100 ns.(d) The free energy difference between two symmetric states along  2 .

Figure S16 .
Figure S16.Errors in p-value and transition timescales from  2 −  1 infrequent metadynamics simulations of residue F4 obtained with bootstrapping.

Figure S18 .
Figure S18.Errors in p-value and transition timescales from  2 −  1 infrequent metadynamics simulations of residue Y10 obtained with bootstrapping.

Figure S20 .
Figure S20.The illustration of bias deposition in the barrier region along  1 during infrequent metadynamics simulation along  2 −  1 of F33 in some of the trajectories in panels (a)-(c).Panel (d) shows an example where bias is not deposited in the TS region.

Figure S21 .
Figure S21.Theoretical Poisson fit to the cumulative distribution of transition timescales obtained from  2 −  1 infrequent metadynamics simulations of residue F33 (a) without repulsive walls, and (b) with repulsive wall potential (UPPER_WALL) at 1 = 2.4 rad to avoid bias deposition on transition state (12 simulations).

Figure S22 .
Figure S22.Errors in p-value and transition timescales from  2 −  1 infrequent metadynamics simulations of residue F33 without walls (upper two panels) and with wall bias (lower two panels) obtained with bootstrapping.

Figure S23 .
Figure S23.Theoretical Poisson fit to the cumulative distribution of transition timescales obtained from  2 −  1 infrequent metadynamics simulations of (a) residue Y21 and (b) F45.

Figure S25 .
Figure S25. 1 and  2 torsion values for residue F22 and F33 in the metadynamics simulation of residue F22.

Figure S26 .
Figure S26.Theoretical Poisson fit to the cumulative distribution of transition timescales obtained from  2 −  1 infrequent metadynamics simulations of residue F22.

Figure S28 .
Figure S28.Representative structures observed in each minimum of  2 -CMAP free energy surface of residue Y23.

Figure S30 .
Figure S30.Errors in p-value and transition timescales from  2 -CMAP biased infrequent metadynamics simulations of residue Y23.

Figure S32 .
Figure S32.The number of hydrogen bonds between residues (a) Y23-R1, (b) Y23-D3, (c) Y23-G56, and (d) representative structure showing the hydrogen bonding of Y23 with G56 during the ring flip.The native orientation of the Y23 aromatic side-chain before flipping is shown in yellow color.

Figure S33 .Figure S34 .
Figure S33.Free energy landscape of  2 −  1 biased metadynamics simulation of residue Y23.White circles indicate the 26 trajectory endpoints immediate after the ring flip.

Figure S35 .
Figure S35.The distribution of side-chain torsions of residues C14 and C38 (disulfide bridge C14-C38) in case of non-converged  2 −  1 metadynamics simulation of residue Y35.The free energy in kcal/mol.

Figure S36 .
Figure S36.The conformations of residue Y35 and C14-C38 disulfide bridge in the native state (red) and unusual conformations (light and dark blue) during well-tempered metadynamics simulations with  2 ,  1 of Y35, and  1 dihedral of a C14-C38 disulfide bridge.

Figure S40 .
Figure S40.The circular correlation coefficient values of torsions  2 and  1 of aromatic residues of BPTI obtained from 1-ms-long trajectory.

Figure S41 .
Figure S41.Theoretical Poisson fit to the cumulative distribution of transition timescales obtained from   infrequent metadynamics simulations of residue F4, Y10 and F33.

Figure S42 .
Figure S42.Theoretical Poisson fit to the cumulative distribution of transition timescales obtained from   infrequent metadynamics simulations of residue F22, Y23, and F45.

Figure S43 .
Figure S43.Results of  2 metadynamics of residue Y10: (a)  2 vs. simulation time, (b)  1 vs. simulation time, and (c) test of convergence along  2 and (d) free energy difference between two identical states.The free energy difference between these two states should be ideally zero as they are identical states.

Figure S45 .Figure S46 .
Figure S45.Diffusion of  2 and  1 during  2 −  1 metadynamics simulation of residue F22 (with CHARMM36m/charmm TIP3P water model) is shown in figures (a) and (b), respectively.(c) Projection of free energy on  2 in time intervals of 100 ns.(d) The free energy difference between two symmetric states along  2 .

Figure S47 .
Figure S47.The decay of eigenvalues with respect to lag time values for F22 and F33 (figures (a) and (b), respectively).The coefficients of descriptors defining eigenvectors  1 22 and  1 33

Table S1 .
The values of  0 and  0 for   atom pairs in CMAP (Equation2)   atom pair ()

Table S3 .
Results of infrequent metadynamics simulations of residue F33 with wall bias and without wall bias.

Table S5 .
and  0 for state-to-state directed InMetaD simulations of residue F22.The transitions to nearest neighbors are considered and five simulations were performed in each case.The symmetry of stable state positions on free energy surface allowed to build a complete transition rate matrix with fewer simulations.The infrequent metadynamics simulations are performed for (i) A to I1 transitions, (ii) I1 to A transitions, and (iii) I1 to I1' transitions.The less stable state I2 (and I2') are ignored in the present calculation.The results are reported in TableS4below.