Probing the Nanosecond Dynamics of a Designed Three-Stranded Beta-Sheet with a Massively Parallel Molecular Dynamics Simulation

Recently a temperature-jump FTIR study of a designed three-stranded sheet showing a fast relaxation time of ~140 ± 20 ns was published. We performed massively parallel molecular dynamics simulations in explicit solvent to probe the structural events involved in this relaxation. While our simulations produce similar relaxation rates, the structural ensemble is broad. We observe the formation of turn structure, but only very weak interaction in the strand regions, which is consistent with the lack of strong backbone-backbone NOEs in previous structural NMR studies. These results suggest that either DPDP-II folds at time scales longer than 240 ns, or that DPDP-II is not a well-defined three-stranded β-sheet. This work also provides an opportunity to compare the performance of several popular forcefield models against one another.


Introduction
Xu et al. have recently studied the nanosecond time scale folding dynamics of a designed threestranded sheet mini-protein [1]. This peptide, called D P D P-II, is one of many peptide sequences originally designed by the Gellman group for the purposes of elucidating the sources of thermodynamic stability and folding cooperativity of beta-hairpin and beta-sheet structures [2]. The most stable of these beta-sheet designs have a scaffold that incorporates successive D-proline and glycine residues ( D PG) in the turn regions, a motif shown to form a stable Type-II´ turn [3].
Of great interest as model systems have been several three-and four-stranded beta-sheet designs from the Gellman group. The first of these, D P D P, was studied by NMR and shown to have cross-strand NOEs and chemical shifts indicative of beta-sheet populations [4]. Syud et al. built upon D P D P, producing (among others): D P D P-II, a three-stranded sheet whose C-terminal hairpin is identical to the N-terminal hairpin of D P D P ( Figure 1); and D P D P D P, a four-stranded composite of D P D P and D P D P-II [5]. The stability of these designs was assessed in a similar fashion by NMR. Figure 1. Designed beta-sheet peptides designed by the Gellman group: D P D P [4], D P D P-II [5], and D P D P D P [5]. D P denotes D-Proline, and O denotes Ornithine.
Recently, designed D PG-turn beta-sheet peptides have become interesting candidates for ultrafast folding beta-sheet systems. Many proteins have been engineered to fold quickly [6][7][8], close to the "speed limit" of folding [9,10]. Upper limits on protein folding rates are thought ultimately to be controlled by the conformational search rate for forming intermolecular contacts [11,12]. For betasheet proteins, the entropic barriers of turn formation are rate limiting [13]. Indeed, a designed variant of human Pin1 WW domain, with a D PG substitution in the turn region, shows a 10-fold increase in folding rate compared to native sequence, up to (~10 s) -1 , becoming one of the fastest folding betasheet proteins to date [13].
It had been hypothesized that because of the reduced conformational entropy of the D PG turn regions, the folding landscape of D P D P-II might not have activation barriers to folding, and shown to be that of a "downhill" folder [8], with kinetics shaped mainly by landscape roughness. Using temperature-jump FTIR, Xu et al. showed that D P D P-II has "the fastest T-jump relaxation rate observed for a beta-sheet system so far" of (~140  20 ns) -1 , with single-exponential relaxation kinetics [1]. More recently, T-jump FTIR studies of the related four-stranded peptide, D P D P D P, show similar singleexponential kinetics, but with a folding time of ~440 ns [14].
While the single-exponential kinetics of D P D P-II can be fit to a two-state Arrhenius-type model, Xu et al. showed that one-dimensional Langevin models of dynamics over a rough free-energy surface [15,16] explain the data equally well, which is their preferred interpretation. In the case of the four-stranded D P D P D P, Xu et al. suspect that many parallel but degenerate refolding pathways may be present [14]. One reason to prefer the "downhill" interpretation is the lack of features typical of activated folding kinetics [15,16]. The fast rate of the relaxation (on the time scale of helix-coil transitions) imply that unfolded and folded ensembles must have a similar degree of compactness, and end-to-end distances as measured by FRET show little sensitivity to temperature. Xu et al. suggest a reason for this is the reduced accessible conformational space imposed by the rigid D PG turns. Smith and Tokmakoff used time-resolved infrared spectroscopy along with site-specific isotopic labeling techniques to show that the D PG turn region of de novo hairpin peptide PG12 does not undergo significant rearrangement upon a temperature-jump, where as the mid-strand regions rearrange on a ~130-ns time scale [17]. Their results support a model where the unfolded state is an expanded but native-like ensemble.
Simulation studies have shed some light onto the thermodynamics of D PG-turn beta-sheet proteins. While there have been no previous simulations of D P D P-II, several groups have simulated the related D P D P peptide, which shares an 11-residue stretch of hairpin residues ( Figure 1). Wang and Sung simulated a 100 ns molecular dynamics trajectory of D P D P using an implicit solvent model, starting from an extended conformation [18]. Their results show D P D P folding to beta-sheet structures, and agree with experimental findings that the D PG turn is more stable than designed three-stranded peptides with NG or GS turns. Roe et al. used replica exchange molecular dynamics in a modified AMBER99 forcefield with an implicit solvation model to sample the thermodynamics of D P D P [19]. REMD (12 replica trajectories each of ~130 ns) dramatically enhanced the convergence of the free energy landscape compared to single-replica MD. The two hairpins of D P D P show simulated populations of ~50% and ~75%, respectively, consistent with NMR and CD studies [4,20], and estimates of thermodynamic cooperativity of -1 to -3 kcal/mol. The less stable of the two D P D P hairpins comprises the C-terminal sequence of D P D P-II.
We have been interested in D P D P-II as a target for molecular simulation for several reasons. It appears to be the fastest-folding beta-sheet protein so far, with relaxation kinetics within the time scale range that can be effectively addressed with all-atom molecular simulation. Moreover, the experimental kinetics remain ambiguous as to whether activation barriers exist for this peptide. To investigate the underlying conformational dynamics, we perform massively parallel molecular dynamics simulations of D P D P-II in explicit solvent.
As we report below, the reaction coordinates of average radius of gyration and solvent-accessible surface area of backbone C=O over time show good agreement with the experimentally measured relaxation rates, but we observe very few three-stranded sheet structures folded within 240 ns, regardless of the forcefield model used. These results suggest that either D P D P-II folds at time scales longer than 240 ns, or that D P D P-II is not a stable well-defined β-sheet, which is consistent with previous NMR spectroscopic data [5].

Results and Discussion
The Folding@Home distributed computing platform [21] was used to simulate molecular dynamics (MD) trajectories, each up to 240 ns in length, for five different forcefields, for a total of ~8.2 ms of simulation. Simulations were performed using the GROMACS simulation package [22], with AMBER94 [23], AMBER96 [24], AMBER99 [25], AMBER99 [26], and AMBER03 [27] forcefields (see Methods section). 1000 total trajectories were generated for the AMBER99 simulations, and 10,000 trajectories each were generated for the other forcefield simulations. Figure 2 shows the distribution of trajectory lengths for each forcefield tested. Figure 2. The distribution of trajectories achieving a given trajectory length, shown for the forcefields tested in this study.

Simulated relaxation kinetics
Time-resolved FTIR measurements cannot directly determine whether a protein is folded, but instead report the status of backbone amide groups, which may be closely related. To best connect with the relaxation rates experimentally measured using FTIR [1], we therefore analyzed the ensemble time course of the total solvent-accessible surface area of backbone C=O groups, as well as the average radius of gyration of the entire molecule. In general, reaction coordinates must be carefully chosen because a poor choice of can yield projection-dependent results [28,29]. The C=O solvent-accessible surface area is a measure that closely connects with the measured amide I band, which is known to be sensitive to hydration status [30]. The radius of gyration is a global quantity that does a good job of characterizing the structural distribution and compactness of a conformational ensemble.
We simulated 1,000 trajectories each (100 each for AMBER99) from 10 different starting configurations ( Figure 3) taken randomly from a high-temperature equilibration trajectory of D P D P-II started from a semi-extended state (see Methods). Figure 4a shows a typical trace of the ensembleaverage radius of gyration over time (for a particular combination of forcefield and starting conformation), which fits well to a bi-exponential curve. Figure 6a shows a typical trace of the ensemble-average C=O solvent-accessible surface area over time. The kinetics also fit well to a biexponential curve. In both cases, the kinetics show a fast equilibration phase (usually  1 ~1-10 ns) and a slower relaxation phase ( 2 ~100 ns). Similar kinetics were computed across all forcefields and starting conformations (Figure 4b). In most cases, the fast phase corresponds to fast equilibration of the starting conformation. Alternatively, in some cases, the fitted values of  1 were extremely short (~0.1 ns), less than the snapshot frequency, indicating that the kinetics may be better described as a singleexponential process with rate constant  2 . Numerical values for all fitted kinetic parameters are shown in the Supplementary Material (Tables S1 and S2).  Regardless of forcefield choice, the slow relaxation times estimated from our simulations are consistent, ranging from about ~60 to ~100 ns for the radius of gyration reaction coordinate (Figures 4b,5), and ~80 to ~150 ns (Figures 6b, 7) for the solvent-accessible surface area. Both compare very favorably to the experimentally measured relaxation time of ~140  20 ns obtained by T-jump infrared spectroscopy [1].   The agreement between simulated and experimental rates is comparable to other contemporary examples of physical kinetics simulations [31]. The slightly faster relaxation rates observed in the simulations may in part reflect the anomalously high diffusion constant of the TIP3P water model [32].
The average radius of gyration at 240 ns across all forcefields is 8.32Å ± 0.47Å, and the average value of the exponential baseline, C, is 8.28Å ± 0.47Å. This reflects a more conformationally expanded ensemble than seen in simulations of D P D P, which showed radius of gyration of ~7Å for unfolded states, ~6.5Å for partially unfolded states, and ~5.5Å for a fully strand-paired native-state conformation [19].

Secondary structures over time
The per-residue secondary structure over time for each forcefield was calculated using the DSSP algorithm [33]. The general features observed across the different forcefields include fast formation of the D PG turn regions, and negligible amounts of sheet formation as quantified by the amount of backbone hydrogen-bonded strand content (see Supplementary Material). It should be noted that strand content may be somewhat underestimated due to the stringent definition required by DSSP.
The amounts of secondary structure across different forcefields reproduce previously noted secondary structural biases [26]. For example, AMBER94 is slightly biased toward more helical conformations and has more populated turn regions (as defined by DSSP), while AMBER96 biased toward beta-sheet conformations, which detectable populations of strand (see Supplementary Material). The more modern forcefields of AMBER99, AMBER99phi, and AMBER03 all show comparable amounts of secondary structural propensities intermediate between AMBER94 and AMBER96. In all cases the DSSP populations are relatively static after ~100 ns.

QH1 and QH2 at 100-150 ns and 200-230 ns over time
To examine hairpin formation, we computed two quantities, Q H1 and Q H2 , reporting the fraction of "native" contacts in (N-terminal) hairpin 1 and (C-terminal) hairpin 2, respectively (see Methods). The quantities Q H1 and Q H2 were used as reaction coordinates to compute the landscape of sampled conformations at two time slices: 100-140 ns and 200-240 ns (Figure 8). Regardless of the choice of forcefield, the conformational landscape mostly disfavors the formation hairpin. Recall that the less stable of the two D P D P hairpins comprises the C-terminal sequence of D P D P-II, corresponding to hairpin 2. With the exception of the AMBER96 simulations, only the formation of hairpin 2 is mostly observed, and only then with a population of 3% or less at 200 ns. For the AMBER96 trajectories, formation of both hairpin 1 and hairpin 2 is observed. For all the simulations, comparisons of the conformational landscape at 100 ns and 200 ns shows very little change in hairpin populations on the ~100 ns time scale (Figure 8).
Folding to a three-stranded sheet is observed for only two out of a total of 10,000 AMBER96 trajectories (Figure 9). One of these two trajectories shows a fully hydrogen-bonded three-stranded sheet structure, while the other shows only hairpin 2 with defined hydrogen bonds, but is otherwise "native" according to inter-residue contacts defined by the Q H1 and Q H2 reaction coordinates. Figure 9. Only two of 10,000 AMBER96 trajectories show folding events for D P D P-II within 240 ns. Shown is the time course of reaction coordinates Q H1 and Q H2 , which monitor the fraction of hairpin 1 and hairpin 2, with conformational snapshots. The second of the two trajectories is "native" by our reaction-coordinate definition, although hairpin 1 does not have a fully hydrogen-bonded structure.
Across all of the forcefields we studied, most all of our simulations do not produce stable threestanded hairpin conformations. We think that this result is very unlikely to be due to poor sampling. With as many as 10,000 simulation replicas per forcefield, there should be a strong likelihood of observing at least some trajectories reaching the folded state [34]. It is possible that forcefield deficiencies may be at work here, but we tested a wide range forcefields, and consistently found negligible amounts of three-stranded. Parallel simulation techniques to accelerate kinetic sampling also has its limits on short timescales where first-passage times are short compared to the folding time [35], but that is not the case here. If the experimentally observed relaxation does indeed correspond to folding, then the overlap in simulated and experimental relaxation time scales should be very favorable for observing transitions to native conformational ensembles. Is D P D P-II a stable folded three-stranded sheet? While Syud et al. reported qualitative NOE data for D P D P-II, this peptide was the least well-folded compared to the other designed sequences in this paper [5]. The measured NMR resonances were weak, and aggravated by poor dispersion, so only key interresidue contacts hinting at the designed structure were reported (Syud and Gellman, personal communication). Combined with our simulation results, this suggests that perhaps D P D P-II is unstable as a three-stranded sheet, and may not be a very relevant model system for studying beta sheet peptides.
Similar plasticity has observed in another designed three-stranded sheet, the betanova peptide [36]. Both betanova and the D PG-turn peptides of Gellman et al. were designed with stable turns and hydrogen-bonded strand regions, to be used as model systems to study beta-sheet cooperativity. WW domains, by contrast, are three-stranded beta-sheet proteins found in nature, whose structures are welldefined [37]. Unlike designed beta-sheet peptides, WW domains additionally possess a conserved network of hydrophobic interactions between their termini. Thus, in general, beta-sheet model systems such as D P D P-II may not have the necessary amount of long-range cooperative interactions needed to fully stabilize their structure.

Conformational clustering and Markov State Model (MSM) analysis
Kinetics-based conformational clustering was performed for all snapshots from the AMBER96 trajectories. The AMBER96 trajectories were chosen as they contained the greatest extent of beta-sheet structure, and the only observed folding events. Our clustering procedure was used to identify five macrostate clusters calculated to be the most metastable, which were used to construct Markov State Models [38][39][40] of the dynamics (see Methods).
We constructed a series of MSMs from matrices of macrostate transition counts, using different lag times ranging from 8 to 240 ns. The performance of these models reveals much about the underlying folding landscape.
The most striking result of the MSM-building procedure was our failure to identify well-separated metastable states that would indicate large activation barriers on the folding landscape. The first indication of this comes from our clustering algorithm, designed to identify the most kinetically metastable states. Only five metastable states were identified, and each contained a broad ensemble of microstate conformations, with average RMSD between any two microstates ranging from 7.3-8.0 Å ( Figure 10).
Regardless of lag time, the spectrum of relaxation rates predicted by the MSM is broad, without a large gap that would indicate a pronounced separation of time scales ( Figure S2, Supplementary Material). These results, at least within the time scale of our simulations, are not inconsistent with either multistate folding or the "downhill" folding interpretation of Xu et al. Moreover, as we increase the lag time used to build the models, the longest implied timescale also increases. If clear activation barriers were present, such that metastable dynamics on the > 10 ns time scale resulted, the implied timescales should level off as the lag time increases. This result is not simply a consequence of poor state definitions, because the kinetic clustering procedure we use should insure that the macrostates are the most metastable basins. Figure 10. Kinetics-based clustering was used to find five maximally metastable macrostates (see Methods). The representative conformations shown for each state are the most probable conformations in that state. Shown next to each representative conformation is the average RMSD between microstates in that cluster, a measure of the compactness of the conformational ensemble, and the number of microstates (of 4000 total) comprising each macrostate.
The variability of simulated relaxations across the ten starting conformations offer an additional indication of the absence of large barriers. This is not only evident from the bi-exponential fits of average radius of gyration and average solvent accessible surface area over time, but also from individual MSMs we built using trajectory data generated from each conformation ( Figures S3 and S4, Supplementary Material). Similar kinds of heterogeneity in relaxation dynamics for different starting conformations have been observed in previous parallel simulations of ultrafast folders [9].
The other striking result of our MSM-building procedure is the unexpected sensitivity of the average C=O solvent-accessible surface area (SAS) to expanded states ( Figure 11). When we use the average SAS of each macrostate to compute a projection of the time evolution of the SAS observable, the effects of averaging over each macrostate is severe enough to produce a signal that increases over time instead of decreasing. When the average SAS is projected onto each microstate, this effect is less severe, yet is still present. The simulation data suggest that the average SAS is more sensitively dependent on expanded conformations that quickly collapse, as compared to more compact conformations. Given our good overall results in recapitulating experimentally observed relaxation rates, we remain confident that our representative set of starting conformations is a useful ensemble to compare with FTIR T-jump experiments. However, the SAS projections underscore the importance of choosing experimental observables that overlap well with the reaction coordinate of interest (in this case, the folding reaction) as to best report the underlying dynamics. Figure 11. Average C=O solvent-accessible surface area (SAS) over time computed from simulation snapshot data (blue), compared to the average SAS of each microstate projected onto the 4000 microstate populations over time (red), and average SAS of each macrostate projected onto 5 macrostate populations over time (green). The differential effects produced by averaging over microstates and macrostates indicate a sensitive dependence of the SAS observable on short-lived expanded conformations.
One question partially addressed by our work is how experiment and simulation might be used to distinguish between two-state vs. "downhill" folding. As we have shown, one potential indication of so-called "downhill" folding from simulations might be a failure to build a Markovian kinetic model able to describe dynamics as transitions between well-defined metastable states. However, our work suggests that perhaps D P D P-II is not a well-defined beta-sheet structure, which brings into question what is meant by "folding" in this case.
Can simulations help suggest experiments that could discriminate downhill vs. activated folding? This is a challenging task, as the observed experimental kinetics for "downhill" folders may depend on many factors. Liu and Gruebele, using one-dimensional Langevin models, present an excellent elucidation of the possible experimental outcomes that can arise from slight differences in folding landscapes (such as native-biases, roughness, and barrier heights) and the reaction coordinatedependence of reporter probes [16]. Using simulations to identify observables that connect well with folding reaction-coordinates may be particularly useful. For example, our simulations of D P D P-II suggest that the C=O solvent-accessible surface area (SAS) is more sensitive to expanded versus compact states. The insensitivity may be in part because the SAS is an aggregate measure across all peptide residues. To the extent that the SAS correlates with the amide I band spectroscopic observable in FTIR T-jump experiments, we suggest that multiple time-resolved FTIR experiments using isotopic labeling of specific residues, combined with microscopic information about peptide conformations from simulation, would help to better resolve folding landscapes for ultrafast folding proteins.

System preparation and simulation protocol
Ten initial starting conformations were selected iteratively from a 1 ns stochastic dynamics (SD) simulation at 3000K, with 9 Å cutoffs for Coulomb and vdW interactions, integration time step of 1 fs, neighbor searching on a grid every 10 steps, at solvent (shear) viscosity of 10 ps -1 . Ten conformations were picked iteratively from a collection of snapshots saved every 1 ps. After picking the first conformation, the most diverse structure (as measured by RMSD) was picked as the next. This procedure was repeated to create a structurally diverse starting set. Each chosen (nearly random) structure was then minimized and equilibrated for the production runs.
Production runs were performed using the TIP3P water model [32] for explicit solvation. A rhombic dodecahedral box of largest dimension 58.7Å was used with periodic boundary conditions. The box contained a D P D P-II molecule with uncapped termini, approximately 4,650 water molecules (this number varied slightly with starting conformation) and two chloride counterions to achieve a net neutral charge. Molecular dynamics (MD) simulations were ran at 308 K in the NVT ensemble with a 2 fs integration time step. The same cut-off and neighbor-list settings above were employed, along with a reaction-field electrostatics model, Berendsen temperature coupling, and constrained bonds with the LINCS algorithm. Trajectory snapshots were recorded every 100 ps. Total C=O solvent-accessible surface area was calculated for each snapshot from the set of all carbon and oxygen atoms in the backbone carbonyl groups, using a solvent probe radius of 1.4Å.

Exponential curve fitting
Best-fit parameters *=(A, B, C,  1 ,  2 ) for bi-exponential curves of the form f(t) = Aexp(-t/ 1 ) + Bexp(-t/ 2 ) + C were calculated for time series of the average radius of gyration and C=O solventaccessible surface area, by using a simulated annealing protocol to minimize the sum of squared errors.
The first 5 ns of the time series were omitted from the fitting procedure. Variances  i  in average radius of gyration at each time point i were calculated by non-parametric bootstrap of 100 samples.
Errors in parameter estimates for each  j were calculated as diagonal elements of the covariance matrix and W is an N x N diagonal matrix of inverse variances: W ij = 1/ i  for i=j, W ij = 0 for i≠j [41].

Secondary structure and "native" hairpin contacts
The DSSP algorithm was used to assess the extent of helix, strand (sheet), turn, and loop secondary structures [33]. DSSP recognizes eight types of secondary structures based on hydrogen bonding patterns: G (3 10 helix), H (alpha helix), I (pi helix), B (beta bridge), E (extended sheet), T (turn), S (loop). We monitor helix content as the total of G, H, I, the strand content as the total of B and E. Q H1 and Q H2 report the fraction of native contacts present for (N-terminal) hairpin 1 and (Cterminal) hairpin 2, respectively. We use the same criteria derived by Roe et al., who used a model of the native conformation to define "native" contacts in each of the two possible hairpins [19]. For hairpin 1, the set of native sidechain contacts (C  for glycine) is defined as residue pairs (R1,I3) A contact between sidechains is defined when centroid distances < 6.5Å and a backbone contact is defined when hydrogen donor-acceptor distance < 2.5Å.

Kinetics-based clustering for building Markov State Models (MSM)
Representative conformations were extracted from the simulation data using a procedure previously described [42], though constant temperature simulations were used. This method uses Markov State Models (MSMs) to identity kinetically related regions of phase space. Thus, two conformations will be found in the same state if a simulation can move between them quickly but will be grouped into different states if transitioning between them is slow. The definitions of fast and slow are based on the timescales observed in the simulations [39,42].
The first step in building such an MSM is to group conformations with a high degree of structural similarity into small sets called microstates. In this study 4,000 microstates were generated based on their all-atom RMSD using a k-centers clustering algorithm [43]. A desirable feature of this algorithm is that the resulting microstates have approximately equal volumes so their populations are directly related to their densities, or free energies. If each microstate is sufficiently small then it is assumed that structural similarity is equivalent to kinetic similarity since it should take a very short time to transition between very similar conformations. Kinetically related microstates, as judged by the number of transitions between them observed in the data, are then grouped together using the PCCA algorithm and this lumping is refined using a simulated annealing scheme [38,44,45]. The center of the most populated microstate from each macrostate is then selected as the representative conformation for that macro state as it is the most probable.

Markov State Model (MSM) construction
The matrix of transition probabilities T between the five macrostates was computed from the trajectory data. The entries of this matrix T ij contain the probability of transitioning from state i to state j in time , which ranged from 8 ns to 240 ns. Diagonalization of (T T -1) produces a set of eigenvalues  k and corresponding eigenvectors e k which describe the dynamics of state populations p(t) as a linear combination of relaxation processes: where  k =[ln  k ]/, and the  i are determined by the initial state populations p(0) [38,39]. Thus  k -1 are the set of implied timescales involved in the relaxation dynamics.

Conclusions
We performed massively parallel folding simulations of D P D P-II to investigate the conformational dynamics underlying its nanosecond refolding dynamics. The simulated relaxation rates, as monitored by average radius of gyration and average C=O solvent-accessible surface area, agree well with the single-exponential relaxation rates experimentally measured by T-jump FTIR. Furthermore, Markov state models built from the trajectory data do not show a separation of metastable timescales consistent with large activation barriers. These results, at least within the time scale of our simulations, are not inconsistent with either multistate folding or the "downhill" folding interpretation of Xu et al. However, despite the agreement with experimental kinetics, we observe very few trajectories that fold to stable three-stranded beta-sheet structures. These results suggest that either D P D P-II folds at time scales longer than 240 ns, or that D P D P-II is not a well-defined three-stranded β-sheet. The latter interpretation is consistent with previous NMR spectroscopic data [5].