GaMD simulations (Miao et al., 2015) were conducted on three different peptide systems: a KPTP tetrapeptide with the Amber ff14SB force field (default barrier), KPTP with a modified Amber ff14SB force field that has a lowered potential barrier between the cis and trans conformers of the peptide bond in the dihedral angle potential energy term (lowered barrier), and a 12-residue segment of ArkA with the lowered barrier (Scheme 1). The residue numbering for ArkA is based on the standard system (Lim et al., 1994). The ArkA peptide simulations were run with a fully extended peptide starting structure, or with a starting structure from the NMR structure of ArkA bound to the Abp1p SH3 domain (Stollar et al., 2009) (PDB: 2RPN). Both starting structures were capped with an acetyl group at the N-terminus and amino group at the C-terminus to neutralize terminal charges. The KPTP simulations were conducted before the ArkA simulations to determine which potential energy function would be used for the peptide bond dihedral in the ArkA simulations to increase the likelihood of capturing proline isomerization. Overall simulation lengths and relevant properties are listed in Table 1.

Simulations were conducted on a Linux cluster using Amber 18 pmemd CUDA for GPU functionality (Case et al., 2018). GaMD was used for enhanced sampling and to overcome the barrier between the cis and trans states of the ω angle (Miao et al., 2015). The Amber LEaP module was used to create topology and initial coordinate files for all simulated starting structures. The simulations were run with the Amber ff14SB force field (Maier et al., 2015). For the lowered barrier simulations, a modification was made to the V ₂ constant in the peptide bond dihedral potential function. This adjustment was made similarly to the way that Doshi and Hamelberg increased the V ₂ constant to increase the barrier between the cis and trans states of proline (Doshi and Hamelberg, 2009). The V ₂ constant in these force fields is solely responsible for determining the energy barrier height for peptide bond (ω angle) isomerization and does not affect other dihedral angles or any other conformational sampling coordinate. V ₂ was decreased from the default of 20 kcal/mol to 15 kcal/mol. Simulations were solvated with the TIP3P-FB water model (Wang et al., 2014), such that the edge of the box was at least 15 Å from any atom in the peptide. Chloride ions were added to neutralize system (1 for KPTP and 3 for ArkA).

Systems were minimized in two rounds with the first one placing a harmonic potential restraint on the peptide of 10 kcal/mol/Ų. Systems were then heated from 100 to 300 K using a 2 fs integration time step for 20,000 steps (40 ps). Equilibration was also performed in two rounds: the first with a 10 kcal/mol/Ų restraint and pressure and density convergence to 1.013 bar with the Berendsen barostat for 50 ps; the second was run without the restraint for 500 ps and the barostat was switched to a Monte Carlo barostat to maintain isobaric conditions. A final equilibration was performed for 52 ns to calculate the GaMD potentials that dictate how the potential energy landscape will be boosted for enhanced conformational sampling during the production stage. Boosts were calculated and applied with respect to the total potential of the system (igamd = 3) with an energy threshold boost biasing value equal to the upper bound of the potential energy (iE = 2). A 9 kcal/mol upper limit was set for the total potential acceleration standard deviation for the purposes of reweighting the produced ensembles (σ _0P). Productions were conducted under isobaric conditions with a 2 fs time step and a 0.1 ps frame output in order to ensure the ensembles could be appropriately reweighted once completed. Independent simulations were started with new random velocities. Bonds to hydrogen were constrained using the SHAKE algorithm during all simulations. The particle-mesh Ewald procedure was used to handle long-range electrostatic interactions with a non-bonded cutoff of 9 Å for the direct space sum.

The NMR structure GaMD equilibration step experienced a random jump in potential energy (Supplementary Figure S1), slightly increasing the applied GaMD energy boost relative to the extended structure simulations (a known problem with the Amber 18 implementation of GaMD). Because of this, the extended structure production steps were rerun with the boost potentials calculated in the NMR GaMD equilibration to maintain consistency between the two sets of simulations and improve sampling of the peptide bond isomerization. The random boost was not substantial enough to create a physically irrelevant potential energy landscape, so any abnormally high energy structure sampled was adjusted for during the reweighting process.

The cpptraj module in AmberTools 18 was used for all trajectory analysis (Case et al., 2018). In house Python scripts were used for all data processing and plots. Visual Molecular Dynamics was used to create structural figures (Humphrey et al., 1996).

Simulations were reweighted to create a potential of mean force (PMF) for the proline ω dihedral angles using both the cumulant expansion on the second order and Maclaurin series expansion to the 10th order provided by the PyReweighting tool kit (Miao et al., 2014). The cumulant expansion on the second order more accurately reproduced the free energy barriers to isomerization, but the Maclaurin expansion method allowed us to reweight independent of the reaction coordinate. Therefore, all other equilibrium data presented from the GaMD simulations have been reweighted using the Maclaurin expansion series to the 10th order. The PyReweighting tool kit was also used to calculate the anharmonicity along the proline ω dihedral for the GaMD simulations. An average ω angle anharmonicity value < 10⁻³ kcal/mol was used as the cutoff to determine whether adequate sampling had been reached (Miao et al., 2015).

Histograms for the ω angles in the KPTP peptide simulations were created to visualize the distributions of angles sampled during the simulations to define ranges for cis and trans values for the subsequent ArkA analysis. Two ω angle states close to the canonical 0° and 180° cis and trans values were observed in the KPTP system (Supplementary Figure S2). Based on the distribution of the ω angles in KPTP, the ranges ω = -90° to +50° and ω = +100° to +240° were used to define the cis and trans states, respectively. Percentages of time each proline spent either in cis or trans as well as the flipping frequencies between them were calculated for KPTP and ArkA simulations.

Proline correlation analysis was performed to compare the likelihood of multiple prolines being physically dependent on one another to isomerize to either cis or trans. All possible combinations of prolines were analyzed for correlation in the KPTP and ArkA simulations. Expected percentages of time each group of prolines sampled cis together were calculated by assuming such events were independent of each other and thus equal to the product of their individual cis sampling times. These expected values were then compared to the actual percentages of time the proline combinations sampled cis simultaneously within the simulations.

Bend, turn, and 3₁₀ helix secondary structures for ArkA were calculated using the DSSP algorithm in cpptraj. PPI and PPII helix structures were calculated based on canonical values of ϕ, ψ, and ω angle ranges as shown in Table 2 (Mansiaux et al., 2011; Radhakrishnan et al., 2012). Running averages of bend, turn, 3₁₀ helix, and end-to-end distance were used to test structure convergence between extended and NMR starting structure systems of ArkA. Every 10th frame (every 1 ps) of the ArkA trajectories was used for testing convergence. Proline ring pucker states were calculated using the χ2 dihedral angle (Cα—Cβ—Cγ—Cδ) and defined by canonical ranges shown in Table 3 (Radhakrishnan et al., 2012).

Replica Exchange MD (REMD) ArkA simulations previously run as a starting point for binding simulations were used to construct a comparison structural ensemble that contains only trans proline conformations (Gerlach et al., 2020). REMD simulation lengths and structures are listed in Table 4. Structures were saved every 25 ps and the first 50 ns of each simulation was not used for analysis. Data from NMR and extended structure starting systems were combined and analyzed together.

CD can be used to differentiate between different secondary structure compositions of a protein. The SESCA CD program (v0.95) was used to calculate theoretical CD spectra for the ArkA simulation ensembles (Nagy et al., 2019). The DS5-4SC1 basis-set was chosen to calculate the theoretical spectra for ArkA, primarily because of its Turn 1 secondary structure definition containing a PPII component, which allows the calculated spectrum to differentiate between ensembles with different amounts of PPII content.

An ensemble containing every 100th frame (10 ps) of trajectory data from the ArkA GaMD simulations was used for the SESCA CD spectra calculations. The REMD simulations were also analyzed with SESCA. Percentages of basis-set secondary structure compositions were also calculated. A scaling factor for the in vitro ArkA CD spectrum was calculated from the REMD ensemble to correct for measurement cell concentration errors that may have influenced the measured CD spectrum. However, to compare all theoretical spectra to the experimental one, and to avoid modifying the experimentally measured values, the inverse of this scaling factor (0.510) was applied to all theoretically calculated CD spectra.

PPII content percentages were also extracted from the CD data using an approximation that calculates the composition of PPII in a peptide, % PPII = θ m a x + 6100 13700 . [1] based on the CD ellipticity (10³ deg cm2/dmol) at the 228 nm wavelength ( θ m a x ) (Kelly et al., 2001). Average PPII percentages from the in vitro and calculated CD spectra were compared to average PPII percentages calculated directly from the simulation trajectory data.

The 12-residue ArkA peptide (KPTPPPKPSHLK) was synthesized by Peptide 2.0 with acetylated N-terminus and amidated C-terminus and purified to >95% purity. The peptide was resuspended and dialyzed against 5 mM Sodium Phosphate buffer, pH 7.0. CD spectra were collected using a Chirascan™-plus CD spectrometer with a 250 µl 70 µM sample at a pathlength of 1 mm. Measurements were obtained at 2 s/point for wavelengths between 180–260 nm with a step size of 1 nm. All spectra were blank corrected with buffer only (5 mM Sodium Phosphate buffer, pH 7.0). Data were collected in millidegrees and were converted to degrees and mean residue ellipticity, [ θ ] M R , according to [ θ ] M R =   100 ∗ ( e l l i p t i c i t y     i n   d e g )   ( c M R ) × ( p a t h l e n g t h   i n   c m ) , [2] where the mean residue molar concentration, c M R , is calculated based on the molar concentration, c, c M R =   ( n ) × ( c ) . [3]

The number of residues, n, equals 12 in the case of ArkA.

Initially, GaMD simulations were run on KPTP using the Amber ff14SB force field (Maier et al., 2015), which contains a peptide bond torsional potential term with a 20 kcal/mol barrier between the cis and trans minima. This potential barrier has been shown to be too low to reproduce the experimentally observed free energy barrier between the cis and trans states (Doshi and Hamelberg, 2009). However, it is still a large barrier to overcome in MD simulations, and even using the GaMD method we only saw isomerization occur on a time scale of around once every 8 ns(Figure 1A). Therefore, to increase sampling and because we were much more interested in the conformations in the cis and trans energy minima than in the energies during the transition, we implemented a lowered potential barrier for peptide bond isomerization of 15 kcal/mol (see GaMD simulations in methods). GaMD simulations of KPTP with the lowered barrier resulted in an increase of proline isomerization frequency of more than an order of magnitude (Figure 1B). Both the default barrier and lowered barrier simulations showed that the proline peptide bond ω dihedral angle can clearly occupy two distinct states that are similar to the canonical cis and trans values (Supplementary Figure S2).

Cumulant expansion on the second order was used to reweight the GaMD ensembles and construct a PMF along the ω angle reaction coordinate (Figure 2A). Very high free energy values at the barriers between the minima for the default barrier simulations indicate that sampling of the transition is insufficient to accurately assess the barrier height. However, with the increased isomerization in the lowered barrier simulations, the barrier height falls closer to where we expect, between 10 and 15 kcal/mol, and the cis and trans minima for both sets of simulations remained relatively equivalent in energy. Average anharmonicity (Anharm_Def = −0.0013_, Anharm_Low = −0.0014) values for both systems were below the cutoff of 10⁻³, indicating that accurate reweighting of the boosted ensemble is possible using cumulant expansion on the second order (Supplementary Figure S3). We also used Maclaurin series expansion to the 10th order for reweighting, which results in a smoother PMF plot, but tends to underestimate free energy barrier heights and can also result in slight shifts in the location of minima (Figure 2B). While the free energy barriers in this PMF are underestimated, the depth and location of the cis and trans minima are similar to the PMF found using cumulant expansion on the second order. Maclaurin series expansion reweighting also has the advantage that the resulting weights are reaction coordinate independent and can be applied to all further analyses of the ensemble. Therefore, we used the Maclaurin series expansion reweighting for all subsequent analyses of both the KPTP and ArkA conformational ensembles. Although the isomerization rates were different for the default and lowered barrier simulations, the percent of the ensemble in cis for each proline (after reweighting) was similar in both sets of simulations (first proline: ∼14% for both, second proline: ∼10% for both) (Supplementary Table S1).

To simulate proline isomerization for the 12-residue ArkA peptide, we chose to use the lowered potential energy barrier to represent the peptide bond dihedral angle since it allowed more frequent isomerization with GaMD and therefore improved sampling of the isomer states. Simulations were conducted using two different starting structures, an extended structure and an NMR structure taken from the AbpSH3 bound state (PDB: 2RPN) (Stollar et al., 2009), to assess sampling by comparing the resulting ensembles. Both starting structures contained all peptide bonds in the trans state. Recently published work using the same extended and NMR starting structures, but using the REMD enhanced sampling method, did not capture any cis sampling (Gerlach et al., 2020), which is consistent with other studies that show REMD is not an effective method for observing proline isomerization (Yedvabny et al., 2014; Neale et al., 2016; Zosel et al., 2018).

In the GaMD simulations, all five ArkA prolines isomerized with a high frequency, sampling both the cis and trans state (Supplementary Table S2). Turn and 3₁₀ helix secondary structure, as well as end-to-end distance, showed that the two sets of simulations started to converge over time, though additional simulation time would be needed for complete sampling (Supplementary Figure S4). Because the extended and NMR starting structure simulations were not completely converged, the extended structure simulations have slightly higher cis content for each proline (Figure 3). However, the cis percentage differences were small enough that we combined the data from all ArkA simulations into one ensemble for further analysis. The PMF and anharmonicity values for both GaMD extended and NMR simulations show sampling was adequate for accurate reweighting of the ensemble (Supplementary Figure S5). These results show that using GaMD with a lowered peptide bond dihedral angle potential energy barrier allows for adequate sampling of the proline cis state for all prolines in ArkA on the nanosecond to microsecond time scale and that this method would likely work for other biologically relevant proline-rich peptides. In contrast to a previous Monte Carlo simulation of the 81-residue disordered transcription factor Ash1 in implicit solvent (Martin et al., 2016), several of the ArkA prolines are in cis in more than 10% of the ensemble. The higher trans occupancy for the Ash1 prolines might be due to the larger size of the disordered region, differences in adjacent residues to the prolines, or differences in the simulation method.

In order to visualize the ArkA conformational ensemble, we plotted the probability density based on the number of prolines in cis and the peptide end-to-end distance (Figure 4). Using these reaction coordinates, the most populated ArkA conformational state has an end-to-end distance of ∼24 Å with no cis prolines (Figure 4). However, the ensemble also contains many conformations with multiple prolines in the cis conformation simultaneously. In fact, with five proline residues that can sample both cis and trans isomers, ArkA only spends 36% (±7%) of the GaMD ensemble with all peptide bonds in trans. There is also an inverse relationship between the number of prolines in cis and how extended the peptide structure is, since the cis conformation introduces a kink into the peptide chain, as has been previously observed for polyproline peptides (Moradi et al., 2009; Radhakrishnan et al., 2012). These kinks and slight end-to-end shortening can be observed in the example peptide snapshots (Figure 4).

The peptide bond dihedral angle of proline affects the overall structure of the ArkA peptide. In the ensemble generated using GaMD, ArkA samples more bend and PPI helix structure than in the previously published REMD ensemble (Gerlach et al., 2020), while spending less time in PPII and 3₁₀ helix (Figure 5). PPI helix was not observed in the REMD systems because it requires a cis ω angle (Figure 5D), but is a stable structure for polyproline sequences when all the prolines are in cis and has been observed in simulations of polyproline (Moradi et al., 2009; William J Wedemeyer et al., 2002). The increase in the bend character of ArkA in the GaMD simulations compared to the REMD ensemble is also likely due to the kinks in the peptide structure caused by the presence of cis prolines. Looking at the dihedral angles of the proline side chain, we also see that the down-pucker state of proline is more prevalent than the up-pucker state in the GaMD ensemble (Supplementary Table S3), consistent with previous studies that show cis isomers favor this pucker state while trans isomers have no preference (Radhakrishnan et al., 2012). The overall PPII helix content of 29 ± 4% is lower than would be expected based on other studies that have quantified PPII propensity based on sequence (Jha et al., 2005; Elam et al., 2013); however, these studies typically do not account for the fact that proline residues in disordered regions will not always be in PPII (trans) conformation since they will sample some cis conformation which cannot be PPII. While our results indicate that PPII content may be lower than previously thought, further experimental studies are needed to determine precise levels of PPII helix content.

Because ArkA contains multiple prolines that can isomerize, we also wanted to test whether there is any correlation between their isomeric states, such as two adjacent prolines being more likely to sample cis at the same time. To do this, we calculated the expected joint probabilities of each pair of prolines being in cis simultaneously assuming that every proline’s isomeric state were independent of the other residues. These calculated probabilities were then compared to the observed simultaneous cis occurrences (Figure 6), but no significant difference was found. We also looked for correlations among more than two prolines, but still no significant correlations were detected (Supplementary Figure S7). This indicates that there is no strong cooperative effect based on proline ring stacking or other interactions that make it more favorable for consecutive prolines to adopt the same structure, and instead, the isomerization is purely stochastic. A similar result has been observed for polyproline sequences in Monte Carlo simulations (Vila et al., 2004; Radhakrishnan et al., 2012). It is possible that other proline-rich peptides might exhibit more correlations between the proline residues depending on the adjacent sequence. It has been shown in several studies that the identity of surrounding amino acids will affect the relative stability of the proline cis and trans isomers (MacArthur and Thornton, 1991; Hamelberg et al., 2005; Brown and Zondlo, 2012). Additionally, more distal sequence effects could be present in the context of the full-length protein. Regardless, the combinatorial effect of uncorrelated proline isomerization suggests that biologically relevant proline-rich disordered sequences will, in general, contain a significant number of structures with cis isomers based on the combined independent probabilities of cis for each proline residue.

CD was used to compare the ArkA secondary structures sampled in our simulated ensembles to those observed in vitro. Although CD can be used to determine the amount of α-helix and β-sheet in a protein, it is difficult to calculate the exact amount of specific secondary structure content in an IDP ensemble based on CD data because there are multiple possible structural ensembles that would be consistent with the same CD spectrum, especially when more secondary structure types besides α-helix and β-sheet are considered (Johnson, 1999; Sreerama et al., 2000; Greenfield, 2007; Louis-Jeune et al., 2012; Micsonai et al., 2015). Therefore, we chose to calculate theoretical CD spectra based on our simulated ensembles using the SESCA program (Nagy et al., 2019). These theoretical spectra can then be compared to in vitro data. The SESCA DS5-4SC1 basis set was chosen for the calculation because it includes a secondary structure category, Turn 1, that includes PPII helix turns, and therefore can be used to differentiate between structures with and without cis prolines.

We calculated theoretical CD spectra using both our ArkA GaMD ensemble and the all-trans REMD ensemble and compared both to in vitro data on the peptide (Figure 7). The calculated spectrum from the GaMD ensemble shows a closer resemblance to the in vitro data in both the wavelength and ellipticity intensities compared to the REMD ensemble spectrum, suggesting that cis proline isomers are present in the in vitro ArkA conformational ensemble. To understand why the spectrum calculated from our GaMD ensemble is closer than the REMD ensemble to the in vitro CD spectrum, we examined the secondary structure percentages that were used to generate the calculated spectra. For both the GaMD and REMD ensembles, “Turn 1” is the most common secondary structure (Table 5). The average Turn 1 character for the REMD ensemble is much larger than the GaMD ensemble since it contains more PPII helix with all of the prolines in trans. Instead, the GaMD ensemble contains more “Turn 2” and “Other” structure, which may include PPI helix and other conformations with one or more prolines in cis.

Ellipticity intensity minimums for the GaMD and REMD ensemble are shifted to a lower wavelength (196 nm) compared to the in vitro ensemble (202 nm). This can be explained by examining the basis set spectra used to construct the theoretical ensembles. According to the basis spectra for DS5-4SC1 in Supplementary Figure S12 of (Nagy et al., 2019), the Turn 1 spectrum has a minimum at 195 nm and maximum at 218 nm, while the Turn 2 spectrum has a minimum at 202 nm and no maximum at 218 nm. The calculated GaMD ensemble spectrum has a smaller minimum at 196 nm and smaller maximum at 218 nm than the REMD ensemble spectrum because it has less Turn 1 content; however, it may still contain more Turn 1 content than the in vitro spectrum. Since Turn 1 includes PPII conformations, this indicates that our calculated spectra for the GaMD ensemble may still contain more PPII conformations and less cis proline conformations than the true in vitro ArkA ensemble. Additionally, the SESCA basis sets were all developed for proteins with peptide bonds in trans and a basis set that explicitly takes into account cis prolines and PPI structure might be able to more clearly identify differences between the in vitro and MD ensembles.

We also used another method to directly calculate the amount of PPII structure in the ensemble from the CD data developed by Kelly and coworkers (Kelly et al., 2001). This empirical method based on data from model proteins with varying degrees of PPII character uses CD ellipticity values at wavelengths of 228 nm to calculate the PPII percentage of the sample (Eq. 1) (Kelly et al., 2001). This equation was used to obtain the PPII percentage from the in vitro CD spectrum and the calculated CD spectra from our GaMD and REMD ensembles. For the MD ensembles, this number was compared to the actual PPII percentage calculated directly from the simulated structures. Based on the ellipticity at 228 nm, the PPII percentage from the GaMD ensemble is relatively consistent with the PPII composition calculated from both the GaMD and the in vitro CD spectra, while the REMD ensemble PPII composition is considerably higher (Table 6). Although Kelly and coworkers caution that their empirically derived equation may not be accurate for all peptides, ArkA resembles the type of sequences (short and proline-rich) that Kelly et al. used in their derivation, and the location of the maximum was consistent at 228 nm in the ArkA CD spectrum. We therefore apply this analysis to our ArkA data, but with the knowledge that it does not constitute a precise quantitative result. Overall, our CD analysis suggests that an ArkA ensemble that contains a significant population of cis proline conformations (as calculated by GaMD) more closely resembles the true in vitro ArkA ensemble.