Can Peptide Folding Simulations Provide Predictive Information for

Aug 23, 2010 - Using replica exchange molecular dynamics folding simulations, we compute equilibrium conformational ensembles of 142 hexa- and decapep...
0 downloads 12 Views 2MB Size
Download PDF
J. Phys. Chem. B 2010, 114, 11899–11908

11899

Can Peptide Folding Simulations Provide Predictive Information for Aggregation Propensity? Edmund I. Lin and M. Scott Shell* Department of Chemical Engineering, UniVersity of California Santa Barbara, Santa Barbara, California 93106-5080 ReceiVed: May 5, 2010; ReVised Manuscript ReceiVed: August 3, 2010

Nonnative peptide aggregation underlies many diseases and is a major problem in the development of peptidebased therapeutics. Efforts in the past decade have revealed remarkable correlations between aggregation rates or propensities and very simple sequence metrics like hydrophobicity and charge. Here, we investigate the extent to which a molecular picture of peptide folding bears out similar relationships. Using replica exchange molecular dynamics folding simulations, we compute equilibrium conformational ensembles of 142 hexaand decapeptide systems, of which about half readily form amyloid fibrils and half do not. The simulations are used to compute a variety of ensemble-based properties, and we investigate the extent to which these metrics provide molecular clues about fibril formation. To assess whether multiple metrics together are useful in understanding aggregation, we also develop a number of logistic regression models, some of which predict fibril formers with 70-80% accuracy and identify aggregation-prone regions in larger proteins. Importantly, these models quantify the importance of different molecular properties in aggregation driving forces; notably, they suggest that hydrophobic interactions play a dominant role. Introduction Peptides in solution can assemble into a variety of complex nanoscale morphologies.1 On one hand, their aggregation into amyloid structures is implicated in several famous diseases, including type-II diabetes and Alzheimer’s,2–7 and aggregation remains a major problem in the development of biotherapeutics.8,9 On the other, these systems are becoming exceptionally versatile platforms for new nanomaterials that self-assemble in response to cues spanning temperature, pH, salt, solvent, and chemical additives.1,10–15 Though a specific need for inhibition strategies motivates many research efforts in this area, a number of findings have been remarkably general; amyloid aggregate phases have almost universal features for a wide range of different polypeptide sequences at nondilute concentrations, often with remarkably consistent structural patterns of interpeptide β-sheets.16,17 Moreover, experiments bear out strong correlations between aggregation rates and simple single-peptide properties like secondary structure preference, hydrophobicity, and net charge.18–20 Such universality and simple relationships suggests a striking and robust connection between physiochemical properties of single peptides and aggregation driving forces. Indeed, it might be argued that the development of aggregation prediction algorithms has greatly benefited from these relationships.21 The seminal work of Chiti et al. showed that a simple linear correlation involving changes in hydrophobicity, secondary structure propensity, and net charge could be used to predict the relative aggregation rates of single-point mutants of several amyloidogenic systems.18 This work was later extended to a model for predicting the absolute aggregation rate based on six parameters; three peptide properties included hydrophobicity, charge, and β-strand hydrophobic/hydrophilic sequence patterning, and three addressed the solution conditions of pH, ionic * To whom correspondence should be addressed. E-mail: shell@ engineering.ucsb.edu.

strength, and peptide concentration.22 It was found that hydrophobicity and hydrophobic-hydrophilic patterning were the strongest contributors to prediction. A related model was also developed to study aggregation-prone regions in several neurodegenerative disease-associated proteins.23 More recently, Tartaglia et al.24 also used such an approach to predict aggregation-prone regions in structured proteins by incorporating amino acid hydrophobicity, net charge, and R/β secondary structure propensities with an estimated intrinsic stability of each residue in the native fold. Caflisch and co-workers have also developed several models for predicting aggregation rates based on simple peptide properties.19,25 Using an approximate model, they show that a simple expression captures relative aggregation rates of mutants as a function of changes in solvent-accessible surface areas, dipole moment, β-sheet propensity, and change in charged and aromatic residues.19 With a more intricate phenomenological model incorporating many of these same properties, as well as concentration and temperature, they are able to predict absolute rates and identify aggregation-prone regions of amyloidogenic proteins.25 Several methods have also made greater use of molecular structures and force fields incorporating both physiochemical and statistical terms to predict aggregation propensity, including the TANGO algorithm described by FernandoEscamilla et al.20 and the 3D profile approach of Thompson et al.26 In particular, the former found that hydrophobicity and β-sheet propensity were the largest contributors to prediction success, although detailed electrostatic and hydrogen-bonding interactions still added detectable prediction value. Several more recent web-based aggregation prediction servers now use related strategies.27,28 While these sequence- and template-based approaches have led to the development of robust prediction strategies, it is less clear why they are successful. Are there emergent molecularscale driving forces inherent in the peptide folding and association free-energy landscapes, and can one systematically detect

10.1021/jp104114n  2010 American Chemical Society Published on Web 08/23/2010

11900

J. Phys. Chem. B, Vol. 114, No. 36, 2010

them? Many groups have now studied the formation of amyloid structures using realistic atomic-resolution molecular simulations with physical (i.e., nonstatistical) potentials.29 In one approach, large proteins were broken into short fragments used to form small oligomers, and oligomer formation free energies30 or order parameters31 probed stabilities in β-sheet conformations that successfully signaled aggregation-prone sequences. Interestingly, the study by Khare et al.30 found that the free energies of dimer or tetramer formation of fragments were not particularly wellcorrelated with β-sheet propensity, hydrophobicity, or net charge. A number of additional studies have used brute-force molecular dynamics or advanced sampling techniques like replica exchange32 to examine the structures, interactions, and underlying free-energy surfaces of small oligomers of short peptides.33–46 These efforts have, in particular, emphasized the roles of electrostatics, hydrophobic interactions, and, to some extent, interactions involving aromatic residues.33,34,45 Simulations using lattice and off-lattice coarse-grained models have also supported these findings.47–58 These studies demonstrate the utility of physical simulations of small peptide systems in inferring mechanisms and driving forces for aggregation; however, to our knowledge, the question of whether or not such approaches also bear out the kinds of simple single-peptide predictors of aggregation propensity used in prediction algorithms has not yet been systematically tested. In principle, peptide folding simulations provide much greater information than sequence alone as they have access to the complete equilibrium conformational ensemble and full atomic description, enabling detailed calculation of both geometric (e.g., solvent-accessible surface area) and ensemble properties (e.g., backbone entropy). Here, we perform a very simple test of the ability of peptide folding simulations in predicting aggregation propensity for a large data set of 142 short peptides. For each peptide, we perform extensive canonical sampling with an accurate force field and implicit solvation model to generate the corresponding configurational distribution, which we then use to compute a number of single-peptide metrics. We examine the ability of these properties in distinguishing peptides that readily form fibrils from those that do not and develop predictive models using a logistic regression approach to quantitatively classify the peptides in these two groups. Our method demonstrates a novel approach to coupling molecular thermodynamic properties with experimental data. Moreover, our results provide an interpretation of the driving forces influencing aggregation and, in particular, reinforce the prominent role of the hydrophobic interaction in this context. Methods Systems and Simulation Methods. We investigate a variety of peptide systems from the literature using molecular simulations. We first examine 120 six-residue systems from the AmylHex database assembled by Thompson and co-workers,26 in which peptides were deemed either fibril-forming or -nonforming based on compiled experimental data in the literature; we randomly selected approximately 60 of each type. These hexapeptides stem from point mutants of the amyloidogenic STVIIE peptide, fragments of Islet amyloid protein, tau protein, insulin, and β-microglobulin. We also investigate 22 ten-residue systems from fragments of human prion protein, human lysozyme, and β-microglobulin that were measured for aggregation propensity by Fernandez-Escamilla et al. using circular dichroism spectroscopy.20 From their 71-decapeptide data set, we choose the 11 that were found to aggregate and an additional random 11 that did not. Note that we expressly construct both

Lin and Shell data sets such that half of the entries are fibril-forming. This approach is designed to balance contributions to the models developed for each group but also, in effect, to ensure that the prior aggregation probability for the subsequent statistical models corresponds to a naïve 50% estimate. The specific peptide sequences for the 142 total systems investigated are given in the Supporting Information. To examine the full equilibrium conformational ensembles of each peptide, we simulate them using replica exchange molecular dynamics (REMD).32 REMD facilitates equilibrium convergence in slow-evolving systems by employing a number of system copies (replicas) that are evolved in parallel at different temperatures, spanning the temperature of interest to a heated state where free-energy barriers are more readily crossed. Periodic conformation swaps of neighboring replicas are attempted and accepted according to

Pacc ) min[1, e∆β∆U]

(1)

where ∆β ) β2 - β1 and ∆U ) U(r2) - U(r1). Here, β ) 1/kBT, and the subscripts 1 and 2 indicate the replicas of interest, with r1 and r2 their corresponding configurations prior to the swap. The entire replica cascade forms a Markov chain that maintains an overall Boltzmann-weighted ensemble at each temperature. In these REMD simulations, our replica temperatures span 270-600 K, and swaps are attempted every 20 ps of MD simulation. At each round of replica swaps, five attempts between neighboring pairs are made in order to facilitate convergence; all swaps in each round are randomized in order. The number of replicas is 20, adjusted so that swap attempts are accepted on average with ∼30-50% probability. Other details of the REMD procedure can be found in our previous publications.59–62 All simulations start with peptides in extended conformations, and trajectory frames are recorded every 1 ps for subsequent analysis. Hexapeptides are run for 10 ns of total time per replica, while decapeptides are run for 20 ns. To model the interatomic interactions in our runs, we use the AMBER 9 program63,64 with the ff96 force field65 and the implicit solvation model denoted “igb)5” of Onufriev, Bashford, and Case.66 Python wrappers are used to manage replica exchange events and analysis. Of note, we and our co-workers have found that the ff96/igb)5 force field-solvent combination stabilizes correct folds of a number of small peptides59,61 and, when combined with a mechanism-based conformational search algorithm, is able to fold several small proteins.60,67 In addition to the ff96 force field, we also adjusted the intrinsic hydrogen radii on basic residues as proposed by Simmerling and coworkers68 to better reproduce the stabilities of ion pairs. Metrics of Monomer Physical Properties. For each REMD peptide run, we evaluate a number of equilibrium averages and properties from the final 1 ns of the simulation trajectory at the lowest temperature. For convenience in presenting the results, here, we give these properties abbreviated names as we describe their computation and potential correlation with physical driving forces. • RG denotes the average radius of gyration, which is computed on the basis of backbone alpha carbons. In principle, for polymeric systems, this parameter is one metric of the solvent quality for a given solute; it is also a metric for the extent to which a peptide adopts a compact, potentially “folded” conformation. • HBOND gives the average number of intrabackbone hydrogen bonds that are formed. The DSSP energy criterion69

Peptide Folding and Aggregation is used to determine when a hydrogen bond is present. The number of hydrogen bonds is strongly correlated with the extent to which a peptide forms secondary structure, or is folded, and is anticorrelated with its ability to make interpeptide β sheets without unfolding. • HELIX measures the average number of (i, i+4) helical hydrogen bonds formed, divided by the maximum possible for an ideal helix of the same peptide length. It varies between 0 and 1. This metric is a direct measure of helical propensity specifically and an indirect measure of barriers to interpeptide β-sheet formation. • SASA measures the average solvent-accessible surface area inÅ2,determinedforeachtrajectoryframeusingtheShrake-Rupley algorithm70 with a solvent radius of 1.4 Å. Similar to RG, this metric probes the folding behavior of a peptide but potentially also the solvation free energy and the scaling of van der Waalstype interactions with other peptides. • HSASA1 measures the average hydrophobic solventaccessible surface area in Å2, on a residue basis. It is determined from the contributions to SASA from all atoms involved in hydrophobic residues. The hydrophobic residues are taken as A, C, F, G, I, L, M, P, V, W, Y. For large solutes, this metric gives a coarse-grained measure of potential hydrophobic interactions upon association. • HSASA2 measures the average hydrophobic solventaccessible surface area in Å2, on an atom basis. It is determined from the contributions to SASA from all carbon and carbonbonded hydrogen atoms. • VOL denotes the average molecular volume of the peptide in Å3, determined for each trajectory frame using a spatial grid algorithm. For small solutes, this parameter is often important to the scaling of the solvation free energy. • ENT is a coarse-grained backbone entropy as defined by Ho and Dill.71 It is determined by assigning each trajectory frame a string of R, β, and l ()“loop”) characters representing the state of peptide residues and then computing the entropy from the populations of different string states. Importantly, it is a whole-ensemble metric, characterizing not just average properties, but fluctuations. This metric measures potential entropic penalties for interpeptide association interactions that require formation of an ordered, amyloid structure. • CHARGE gives the net charge of the peptide at neutral pH, determined from the number of acidic and basic groups in the sequence. Importantly, it depends on sequence alone and is independent of the conformational ensemble. This metric is important to the electrostatic component of the solvation free energy and the electrostatic repulsion between associating peptides and plays a role in determining the hydrophobicity of a peptide. • CGROUPS gives the total number of charged groups in the peptide, which may be positive when CHARGE is 0 if there are as many acidic as basic residues. Like CHARGE, this metric is only dependent on sequence. The collection of metric data for a specific peptide sequence can be thought of as a multi-element vector, whose entries we denote as M i with i ) 1, 2, ..., 10. The collection of all such data for the hexapeptide simulation data set will be written as M hexa ik , where k ) 1, 2, ..., 120. Similarly, for the decapeptide data set, we have Mdeca ik , where now k ) 1, 2, ..., 22. Statistical Analysis. To understand the potential role that the computed peptide metrics may play in providing clues to aggregation, we perform a classification analysis to identify those metrics that are most predictive in determining whether or not a peptide will aggregate. This approach is similar to an

J. Phys. Chem. B, Vol. 114, No. 36, 2010 11901 earlier study involving one of us that was designed to predict native residue-residue contacts from protein fragment simulations.62 We propose a logistic regression model in which the probability that a particular sequence will aggregate is given by the following expression

(

ln

)

Pagg ) R0 + 1 - Pagg

∑ RiMi

(2)

i)1

where Pagg gives the classification probability that the sequence belongs to the aggregation class, R0 and the Ri are yet unknown coefficients, and Mi are the simulation-determined metrics for the sequence. Solving this expression for Pagg gives

exp(R0 + ΣRiMi) i

Pagg )

(3)

1 + exp(R0 + ΣRiMi) i

The unknown parameters R in the model are found by likelihood optimization. Separate models are developed for the hexa- and decapeptide data sets. To begin, the set of peptides is first randomly divided evenly into training and test sets; the training set is to be used to determine parameter values, while the test set is used to evaluate their reliability and filter for overfitting. The training set is used to find R values by maximizing the probability that the model correctly classifies fibril-forming and -nonforming peptides as so. This training set likelihood is given by

L)

C (1 - Pagg,k)1-C ∏ Pagg,k k

(4)

k

k

where the product proceeds over all sequences in the training set and Ck is a class identifier that equals 1 if the peptide is a fibril-former and 0 if not. The effect of the Ck in the exponent is to select which probability to use in the likelihood maximization for a peptide, Pagg if it is known to form fibrils readily or (1 - Pagg) otherwise. Substitution of eq 3 into 4 and taking the logarithm yields

ln L )

∑ Ck(R0 + ∑ RiMik) - ln[1 + exp(R0 + k

i

∑ RiMik)]

(5)

i

This is the log-likelihood objective function to be maximized as a function of the R’s. In practice, we use simple NewtonRaphson gradient minimization to optimize this nonlinear equation. After model parametrization, we use the test set to characterize its performance. It is important to note that none of the optimization involves test set data, and therefore, this represents, in some sense, an independent assessment of the transferability of the parametrized model. As a simple metric, we compute the fraction of members of the test set that are correctly classified by the model. That is, we compute Pagg,k for each test system and deem it correctly classified if Pagg,k > 1/2 for Ck ) 1 and Pagg,k < 1/2 for Ck ) 0. To ensure that the specific choice of training and test sets does not bias our results, we repeat the entire optimization procedure 100 times, each one entailing a different random assignment of peptides to the two sets. From each optimization,

11902

J. Phys. Chem. B, Vol. 114, No. 36, 2010

Figure 1. Mean values of metrics computed from REMD simulations. Error bars give the standard deviations of the metrics across all peptide sequences. Some of the metrics have been scaled as indicated in the labels at the bottom of the figure. Units are given in the discussion in the Methods section.

we compute a set of optimal R’s. Rather than average this set of coefficients over the trials, which may generate artificial errors in the subsequent model, we take the parameters that lie at the centroid of the 100 sets in parameter space. To evaluate distances in parameter space, we nondimensionalize each Ri by multiplying by its respective standard deviation of Mi over the data set. Finally, we compute the fraction of correctly classified peptides for each trial using its test set and compute an average fraction for evaluation of the final model. Results Replica Exchange Folding Simulations. The results for the metrics computed from the REMD simulations of the 142 peptides investigated are summarized in Figure 1. Here, the average values and standard deviations of each metric over the entire data set are shown. It is interesting to consider the variation of the computed metrics over the peptide data sets. For the hexapeptides, the CHARGE, atom-based HSASA2, and CGROUPS have the largest relative variance, while the SASA, VOL, and RG have the least. It is of note that the helicity values for the hexapeptides are always 0 as these short systems did not form any helical hydrogen bonds during the course of the REMD runs. For the decapeptides, the metrics with the highest relative variance are CHARGE, HELIX, and CGROUPS, and those with the lowest are again SASA, VOL, and RG. These results may partially reflect the fact that such short peptides have a limited range of conformers and hence a low degree of structural diversity despite substantial differences in hydrophobic and electrostatic groups. A comparison of the average values of the metrics between the fibril-forming and nonforming subsets of the data provides an initial assessment of their potential as classification indicators. Figure 2 plots the difference for each metric between these two groups as a fraction of the standard deviation over the entire data set, such that positive values indicate larger averages for the fibril formers. For the hexapeptides, the metrics varying the most from formers to nonformers are HSASA1 and CGROUPS, while for the decapeptides, it is HSASA1, HELIX, and ENT. A significant trend is that the signs of the metric differences are consistent for both hexa- and decapeptide groups. Those metrics that seem to have, on average, higher values for fibril formers are RG, SASA, HSASA1, HSASA2, and VOL. At the same time, fibril formers tend to have lower values of HBOND, ENT, CHARGE, and CGROUPS. The suggestion of these initial results is that peptides readily aggregating into amyloid

Lin and Shell

Figure 2. Difference in average metric values between fibril-forming and -nonforming sequences. Each value represents the average over all formers in the data set minus that of the nonformers, divided by the corresponding standard deviation over all peptides.

structures tend to be larger (for a given sequence length) and have more exposed hydrophobic surface area but also have fewer intramolecular hydrogen bonds, a lower backbone entropy and net charge, and fewer charged groups. These results at least seem partially consistent with a molecular physical picture of aggregation driven by hydrophobic and hydrogen-bonding interactions. Namely, exposed hydrophobic groups and unsatisfied hydrogen bonds can serve as stabilizing interactions for interpeptide association. Similarly, charged amino acid residues tend to increase solubility and decrease hydrophobic interactions, while peptides that are larger can display an increased interaction surface. The significance of lower backbone entropies in fibril formers is less immediately clear. One possible explanation may be that conformationally flexible and unstructured peptides may have to surrender too unfavorable of an entropy loss upon incorporation into aggregates, while ones that are relatively stiff and in a conformation compatible with the amyloid structure more readily add to growing oligomers. Such a trend makes sense for short peptides, where monomer conformations are likely to be in near-extended form that would not require substantial refolding in order to be incorporated. Indeed, Figure 3 shows a Ramachandran plot of the backbone torsional angles observed in the data set from the dominant cluster structure for each sequence, and these results show a very high degree of extended structure in the hexapeptides and a significant amount in the decapeptides. On the other hand, this line of thinking might not apply to larger fragments, where a high degree of stability in a specific fold may reduce backbone entropy but simultaneously disfavor oligomerization. One unexpected trend in these results is the on-average higher values of HELIX for the fibril formers. This propensity seems at odds with earlier work by a number of groups. Dobson and co-workers,18 for example, found that mutants of several systems that led to increased helical propensity reduced aggregation rates; moreover, work by Kallberg et al.72 suggested that R-helical regions in proteins that were predicted to have β propensity when excised from the full sequence could harbor aggregation “hot spots”. Our apparently conflicting results may be the result of several factors. One possibility is that force field errors may impart erroneous high helicity to the systems studied, although this seems unlikely. Though there has been much debate about R/β balance in protein force fields, the one employed here has not been found to be too helical and in fact has been criticized in the past as favoring β structures.73–76 Our previous work with this force field also does not suggest problems of this sort as we have shown a good degree of balance in folding a variety

Peptide Folding and Aggregation

Figure 3. Ramachandran plots for all sequences and amino acid residues simulated. (top) Results for hexapeptides and (bottom) results for decapeptides. Phi-psi angles are not included for terminal residues in each sequence.

of systems.59,61 Perhaps a better explanation could be that we did not include a broad enough range of sequences with different helical propensities in the relatively small data set to attain an accurate measurement of this effect that overcomes statistical noise. Indeed, the helical values measured are quite small, averaging 5% with a maximum of 34%. That is, none of the peptides exhibited particularly significant populations of helical structures. A third possible explanation is that, for short peptide systems, helicity may play a lesser a role as a deterrent for aggregation than has been assessed in studies of longer peptides. Indeed, short peptide fragments without well-defined folds may explore a range of conformations. To what extent do the measured metrics contain independent information from one another? In Table 1, we evaluate the

J. Phys. Chem. B, Vol. 114, No. 36, 2010 11903 correlation coefficient for metric pairs over the data set. Most of the metrics pairs are uncorrelated; however, there are a few for which a significant correlation (>0.4) exists. HSASA1 is strongly anticorrelated with CGROUPS, which is not surprising since both metrics involve specific classes of amino acids that are mutually exclusive in the metric definitions employed here. Trends in geometric parameters are also evident. VOL is wellcorrelated with the surface area metrics, particularly SASA. HELIX also displays a slight correlation with HBOND, as one might expect. The last major connection is an apparent anticorrelation of RG with HBOND, indicating that structures with a smaller radius of gyration tend to have more intermolecular hydrogen bonds. This is not terribly surprising, except for the fact that HBOND is not particularly well correlated with molecular volume or any of the surface area metrics. Presumably, RG is better than these other variables at indicating peptides in extended, noncollapsed conformations, where hydrogen bonds are unable to form. Detecting Molecular Signatures with Classification Analyses. To better build a connection between simple molecular properties and related driving forces for aggregation propensity, we use the maximum likelihood approach described in the Methods section to develop a number of logistic regression models to quantify the joint relevance of parameters to the picture. We consider systematically all possible models including one, two, and three metrics. These models are developed separately for the hexa- and decapeptide data sets. Each model is scored by its average classification success, that is, the average fraction of peptides in the test sets that are correctly identified as either fibril-forming or not. Figure 4 shows the classification success for each metric, averaged over all models containing that metric. These numbers give an initial assessment of the ability of the metrics to provide predictive information about aggregation propensity. The relatively narrow range of these classification success scores, roughly 56-67%, reflects averaging over all logistic regression models, regardless of quality. It is important to bear in mind that naïve, metric-free prediction based solely on the prior (i.e., on the relative populations of fibril-forming systems) gives a classification success of 56 and 50%, for the hexa- and decapeptide data sets, respectively. Thus, it is apparent that most of the parameter combinations have little predictive capability and thus do not suggest a simple molecular picture. Still,

TABLE 1: Correlation Coefficients for Parameters Measured from REMD Simulations of the Peptide Data Seta

a

The panels above the diagonal give values for the hexapeptides simulated, while those below correspond to the decapeptides. The highlighted regions indicate values with magnitudes above 0.4.

11904

J. Phys. Chem. B, Vol. 114, No. 36, 2010

Figure 4. Average classification success for models as a function of the metric included. The results show the average over all one, two, and three-parameter models including each metric, as shown.

Figure 5. Classification success for regression models as a function of the number of metrics used. Both the best model (highest classification success) and average over the top 10% of all models for a given number of metrics are shown.

examining Figure 4, some metrics do appear to contain more information than others. The residue-based hydrophobic surface area, HSASA1, is among the most predictive metrics for both data sets, according to this analysis. CGROUPS (for the hexapeptides) and HELIX (for the decapeptides) also consistently result in increased classification success when included in models. Importantly, these findings are consistent with the analysis of Chiti et al.18 and later related studies.19,22–25 However, the models involving HELIX infer the opposite trend than expected; the coefficients determined in the regression models tend to favor helicity for aggregation-prone sequences. This result is due to the fact that higher helical content was found, on average, in the fibril-forming members of the decapeptide systems. As the helicity values themselves were quite low, as discussed before, the role of HELIX in the models may be limited to this specific data set. Figure 5 shows the improvement to the regression models as more parameters are added. Relatively little improvement is seen in going from two to three parameters, and thus, we do not consider models with four or more. The best models are able to achieve about 71 and 82% classification success for the hexaand decapeptides, respectively. For the former, this is about 1-3% higher than the average success over the top 10% of models using different metrics. The top hexapeptide models, as given in Table 2, make consistent use of CHARGE and CGROUPS to predict aggregation propensity, with SASA, HSASA1, and VOL emerging as secondary predictors. The regression coefficients for the former two metrics are consistently negative, indicating that the models find peptides with more negative charge and an increased number of charged groups less likely to form fibrils. On the other hand, the latter

Lin and Shell three parameters most frequently have positive regression coefficients, in agreement with the earlier finding that aggregating hexapeptide sequences in the data set, on average, are larger and have greater solvent-accessible surface area. For the surface area metrics, however, a few of the top 10 models manifest negative coefficients; this may be the result of balancing their effect with other metrics to which they are highly correlated, such as VOL. It is interesting to note the presence of two two-parameter models among the most successful hexapeptide regressions. Both of these include CGROUPS, in combination with VOL and SASA. The appearance of CGROUPS in these is consistent with its large relative variation between fibril formers and nonformers in the data set, as shown in Figure 2. However, the absence of HSASA1 in these is surprising given that it had the largest relative variation. One possible explanation for this result is that pairs of metrics not involving HSASA1 were less correlated and provided more independent information for classification. Indeed, Table 1 shows that HSASA1 has a larger correlation with CGROUPS than does VOL or SASA. In contrast to the hexapeptide regressions, the decapeptide models give evidence that a few specific combinations of metrics result in particularly high classification success. The combination of HELIX with HSASA1 provides the most predictive models, with both metrics weighted positively for fibril formers. As discussed before, the somewhat surprising effect of helicity here stems from the higher values found from simulations of aggregation-prone peptides in the data set. Interestingly, this most predictive model includes only two parameters, and the next five top-performing ones include another metric in addition to these two. The decrease in classification success with these additions is likely a result of the small number of decapeptides studied as each training and testing subset includes only 11 parameters. As more metrics are included, the classification success on the test subsets suffers from overfitting to the training data. It is also noteworthy that the decapeptide models have quite large values for the determined coefficients, a likely indicator that the models are sensitive to the small data set. Of the decapeptide models not involving HELIX, the top regressions are able to achieve 71-73% classification success, not significantly higher than the best hexapeptide ones. Unlike the hexapeptide results, however, these models all include HSASA1, underscoring a significant predictive power in the hydrophobic surface area. An analysis of all models involving HSASA1 shows that ENT, CHARGE, and CGROUPS tend to provide the most predictive power in combination with it. In particular, a simple two-metric model involving HSASA1 and CHARGE is the best of these at 73%. Application to Larger Proteins. As a very simple test of the transferability of the analysis developed here, we use them to identify so-called aggregation-prone regions in two larger systems, the Aβ1-42 peptide and the 305-336 fragment of the tau protein, both implicated in the pathogenesis of Alzheimer’s disease. Proline-scanning mutagenesis, NMR, and fragmentation experiments with Aβ1-42 suggest that residues 15-21 and 30-42 are particularly critical to fibril formation.77–80 Similar efforts with tau have shown that the residue range of 306-311 forms a minimal interaction motif for fibril formation.81 These two systems have also been standard tests for algorithms predicting aggregation-prone regions of protein sequences.20,24–27,31 Here, we treat each system by first breaking it into all possible fragments of 6 to 10 residues long, separately performing REMD, and computing metric values for each. We then use the respective hexa- and decapeptide models to predict the fibril-

Peptide Folding and Aggregation

J. Phys. Chem. B, Vol. 114, No. 36, 2010 11905

TABLE 2: Top 10 Regression Models Developed for the Hexapeptide Data Seta model

prior

1 2 3 4 5 6 7 8 9 10

-16.7 -13.7 -7.4 -20.5 -10.8 -11.0 -18.1 -13.6 -18.4 -2.6

RG

HBOND

HELIX

SASA

HSASA1

HSASA2

VOL

ENT

0.9 1.0 0.6 1.2 3.3

-7.2 0.7

1.1 1.1

CHARGE

CGROUPS

classification success

-1.1 -1.2

-2.0 -1.9 -1.4

0.71 0.71 0.68 0.68 0.68 0.68 0.68 0.67 0.67 0.67

-0.6

0.7 6.7 0.5 -0.7 -0.5 1.1

-1.5

1.1 -1.1

-1.2 -1.3 -1.2

a

Only models containing up to three metrics were developed. Except for the prior, the values given indicate the model relevance, that is, the coefficient in the model normalized by the standard deviation of the corresponding metric over the peptide data set. These values are taken from the centroid of the coefficients developed from 100 trial optimizations. Classification success indicates the average fraction of peptides in the test set that were correctly classified using the model.

TABLE 3: Top 10 Regression Models Developed for the Decapeptide Data Seta model 1 2 3 4 5 6 7 8 9 10 a

prior

RG

-143.7 -1279.5 -1472.9 -453.9 -30.9 -331.6 -445.5 -20.3 -948.8 30.1 78.4 -15.2

HBOND HELIX SASA HSASA1 HSASA2 VOL 25.1 277.0 648.9 169.9 54.5 123.0 -431.5

558.6

64.9

49.3 153.8 540.0 246.9 68.2 250.3 10.2 467.6 6.6 13.9

ENT

CHARGE CGROUPS classification success

-477.7 -72.6

9.4 -7.1 -13.8 2.4

-2.6

0.82 0.80 0.79 0.79 0.76 0.74 0.73 0.72 0.71 0.71

See the caption for Table 2 for an explanation of the values.

forming classification probability, averaging over the top 10 models, as reported in Tables 2 and 3. We then consider these average classification probabilities to be aggregation propensities that vary between 0 and 1. Each average propensity is assigned to the mean residue number of the fragment (e.g., fragment 1-6 is assigned to residue 3.5). Figure 6 shows the results of application of the models to these two systems, with the experimentally indicated aggregation-prone regions highlighted for comparison. The hexapeptide models do a good job in predicting the aggregation hot spots in both peptides; a propensity cutoff of at least 0.7-0.8 successfully places all model predictions within the known regions. In both cases, however, there are regions with an average propensity, corresponding to an average classification probability, that is higher than 50% but that is outside of the aggregation-prone regions, near residue 13 for Aβ and 325-332 for tau. Such results may suggest that the priors computed during model optimization may be inaccurate for general predictions, but they might also indicate that the models are only capable of successfully classifying regions with either very small or large propensities, with intermediate values being more indeterminate. The variations of the top 10 models are also significant, as shown by the error bars in Figure 6. Thus, there is likely some gain in statistical robustness in averaging over the top models, which would tend to reduce the effects of bias due to the size of the training data set. The decapeptide predictions for the aggregation-prone regions in Aβ and tau appear to be less accurate. If a cutoff of 0.6 or above is used to indicate aggregation-prone regions, the decapeptide models are able to discriminate the 15-21 hot spot in Aβ but significantly underpredict the full sequence range at the C-terminus of the peptide. At the same time, a high aggregation propensity is erroneously assigned to the first 10residue fragment in the sequence. An analysis of the individual

Figure 6. Prediction of aggregation-prone regions in (top) Aβ peptide and (bottom) tau. The experimentally determined regions are highlighted in gray. The aggregation propensity Pagg for the hexapeptide models is determined by breaking the protein into all six-residue fragments, simulating each with REMD, using computed metrics with the top 10 hexapeptide regression models, and averaging results. Error bars show the standard error over these 10 models. A similar procedure is used for the decapeptide approach.

metrics shows that a high degree of helicity (HELIX ) 0.47) is found for this particular fragment, which increases its propensity in accord with many of the top 10 models. Interestingly, however, high HELIX values are also found for the hot spot range of 15-20 (mean value ) 0.35) but not for any other

11906

J. Phys. Chem. B, Vol. 114, No. 36, 2010

Figure 7. Average hydrophobic surface area per residue (metric HSASA1) for (top) Aβ peptide and (bottom) tau. The experimentally determined regions are highlighted in gray.

residues in the sequence. For tau, the decapeptide models do a poor job of identifying the fibril-forming motif at the beginning of the sequence range; instead, they predict an aggregation hot spot for residues 312-319. Again, these results emerge from high HELIX values for this residue range (mean value ) 0.40). For both systems, the variations in classification predictions among the top 10 models for the decapeptides are more significant than those for the hexapeptides, as shown by the larger error bars in Figure 6. This is compatible with the much smaller data set used in the generation of the former regressions. From where does most of the predictive ability of the models stem? An analysis of the individual metrics measured for fragments along each sequence shows that HSASA1 is noticeably higher in aggregation-prone regions. Figure 7 shows this quantity on a per-residue basis across the sequences. In the experimentally identified parts of the sequence, the hydrophobic surface area measured from hexapeptide fragments is significantly higher for both Aβ and tau. Even though this metric is only positively weighted with aggregation propensity in 3 of the top 10 decapeptide models, its correlation with VOL and anticorrelation with CGROUPS, which appear in all but one of the other top 10, appears to explain the predictive power of this quantity. For the decapeptide fragments, variations in HSASA1 across the sequences are less significant due to averaging larger stretches of sequence, but this metric still shows increased values in the aggregation-prone regions. The sum of these results suggests the possibility of a critical value of exposed hydrophobic surface area, specific to fragment length, after which regions of sequence become highly aggregationprone. The particularly high predictive power of the hydrophobic surface area alone, without other sequence-specific metrics, may also be suggestive of the origins of universal features of aggregation. Discussion and Conclusions In this work, we have mined folding simulations of short peptides for information about aggregation propensity. Our central aim has been to clarify whether or not there are clues to

Lin and Shell aggregation inherent in the equilibrium conformational states of these systems by way of computing from them structural-, sequence-, and ensemble-based metrics for a large number of different amyloid fibril-forming and -nonforming sequences. We find that there are some notable trends in peptides with high aggregation propensity; in particular, these sequences on average tend to display a greater hydrophobic surface area and have a higher number of charged amino acid residues. These findings are mostly consistent with the most predictive factors in a number of sequence-based algorithms for predicting aggregation rates and propensities.18,19,22–25 We have also developed a number of predictive classification models that help identify those combinations of metrics with the most useful information for signaling aggregation. The best such models are able to achieve ∼70-80% success rates in correctly classifying peptides as fibril-forming or not. On the whole, our models are not as predictive as many of the sequenceor template-based approaches in the literature and developed as webservers. This is likely due to the fact that the metrics feeding into these predictive approaches contain only singlepeptide properties and do not speak directly to the aggregate state. For example, we do not consider sequence hydrophobichydrophilic patterning of the kind that would facilitate amyloidlike β-sheet arrangement, we do not use amyloid structures as profiling templates, nor do we consider interpeptide interactions of any sort. Such considerations would improve the models determined here but may require significant additional simulation expense for larger REMD systems. On the other hand, our main goal was to develop a systematic approach for identifying molecular determinants (and related driving forces) for aggregation, and it is nonetheless striking that single-peptide properties alone contain enough information to reach prediction rates at this level. The particularly high predictive power of the hydrophobic surface area in these results warrants some discussion. Specifically, this metric appeared in a number of the best-performing logistic regression models, and it also seemed able to signal, on its own, many of the aggregation-prone regions in Aβ and tau peptides. Such results follow a number of studies that have established the basic importance of the hydrophobic interaction in aggregation and may suggest a robustness of prediction algorithms to model sophistication. At the same time, however, predictions based primarily on this interaction may not be able to correctly identify fibril-forming sequences with few hydrophobic residues, like the GNNQQY peptide from the yeast prion protein Sup35.82 The suggested picture is thus one in which hydrophobic interactions provide a significant driving force that is modified by more specific interactions upon oligomerization and assembly. Interestingly, only the residue-based hydrophobic surface area measurements used in this study provided significant predictive power. The atom-based measure, on the other hand, probably used too fine grained of a decomposition of the solvent-accessible surface to adequately capture which regions effectively appeared as hydrophobic. A potentially more powerful approach to determining which parts of the surface may be hydrophobic, which may in turn facilitate prediction of aggregation propensity, may be the recent proposal by Godawat et al. to use cavity formation probabilities in simulated hydration-layer water.83 Two measures of electrostatic interactions, the net charge and total number of charged groups, also signaled in which class peptides were likely to fall, with a larger absolute charge and number of charged groups disfavored for fibril formers. Naturally, the presence of charged amino acids reduces the

Peptide Folding and Aggregation amount of hydrophobic surface area, both by definition in the current measures and also from a physical standpoint with regard to interfacial water. Thus, to some extent, their predictive power may also be tied to hydrophobic interactions. However, the existence of charged fibril-forming peptides in the data set may require a more nuanced picture of the role of electrostatic interactions. Indeed, Lopez de la Paz et al. suggested a complex picture of electrostatics as driving forces in experiments on hexapeptides with different charge distributions, modulated with pH effects and capping groups. Interestingly, they found that some sequences with (1 net charges more readily formed fibrils than netural ones, and by adding both charged termini to neutral, capped peptides could induce fibrillization in one system. However, peptides with larger net charges, above +2, did not form fibrils. Their results might suggest that electrostatic interactions might actually facilitate the alignment and incorporation of peptides into amyloid structures at moderate to low strengths but hinder it due to a increase in solubility at higher values.84 In our hexapeptide studies, the emergence of both CHARGE and CGROUPS in some of the top regression models may reflect this kind of balance. An alternative approach for future work may be to use as a predictive metric a computed, solvent-mediated electrostatic energy between two peptides or the thermodynamic properties and structure of water solvating a single peptide. Our finding of fibril-forming decapeptides with helical content, both in the original data set and in scanning fragments of Aβ, seems at odds with early results18 and requires close interpretation. The amount of helical content in these was just moderate (maximum ≈ 35%), and no peptides folded to full, highly stable helical structures. Indeed, peptides of this short length are unlikely to be highly stable in any particular fold, instead experiencing a diverse conformational ensemble. Still, that we found fibril-forming peptides with at least some helical propensity questions the utility of this quantity as systematically disfavoring aggregation. Importantly, our predictions of helicity were not based on database-trained secondary structure prediction servers, which might tend to overpredict structure in these small systems, but rather on the simulated conformational ensemble under the action of an accurate force field. The potential irrelevance of helical propensity to aggregation driving forces might also receive support from the number of other prediction methods that use instead only β propensity or β sequence motifs.19,25,85 Finally, this work underscores one of the main challenges in building a basic understanding of aggregation propensity, the difficulty associated with using metrics that are independent. As emphasized in the discussion above, we found significant correlations between different metrics across the data set, in particular, between the hydrophobic surface area and charge metrics and among those involving the geometric quantities like surface area, radius of gyration, and volume. This kind of difficulty is well established,21 but it might also suggest the need for physically relevant variables that more naturally integrate these different quantities. Solvation-based analyses and stabilities and structures of small oligomers of the kind pursued in several atomic-level simulations33–46 may yield metrics for accurate predictive models. Currently, such routes are computationally expensive as tools for quick, webserver-type assessments of aggregation, although coarse-grained models may provide a viable approach.47–58 Acknowledgment. The authors gratefully acknowledge the support of the Camille and Henry Dreyfus Foundation and the National Science Foundation (Award No. CBET-0845074).

J. Phys. Chem. B, Vol. 114, No. 36, 2010 11907 Supporting Information Available: Hexa- and decapeptide sequences used in this work. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Chen, P. Colloids Surf., A 2005, 261, 3–24. (2) Dobson, C. M. Trends Biochem. Sci. 1999, 24, 329–332. (3) Dobson, C. M. Nature 2003, 426, 884–890. (4) Thirumalai, D.; Klimov, D. K.; Dima, R. I. Curr. Opin. Struct. Biol. 2003, 13, 1–14. (5) Chiti, F.; Dobson, C. M. Annu. ReV. Biochem. 2006, 75, 333–366. (6) Eisenberg, D.; Nelson, R.; Sawaya, M. R.; Balbirnie, M.; Sambashivan, S.; Ivanova, M. I.; Madsen, A. O.; Riekel, C. Acc. Chem. Res. 2006, 39, 568–575. (7) Fernandez-Busquets, X.; de Groot, N. S.; Fernandez, D.; Ventura, S. Curr. Med. Chem. 2008, 15, 1336–1349. (8) Wang, W. Int. J. Pharm. 2005, 289, 1–30. (9) Weiss, W. F.; Young, T. M.; Roberts, C. J. J. Pharm. Sci. 2009, 98, 1246–1277. (10) Mihara, H.; Takahashi, Y. Curr. Opin. Struct. Biol. 1997, 7, 501– 508. (11) Zhang, S. Nat. Biotechnol. 2003, 21, 1171–1178. (12) Rajagopal, K.; Schneider, J. P. Curr. Opin. Struct. Biol 2004, 14, 480–486. (13) Fairman, R.; Akerfeldt, K. S. Curr. Opin. Struct. Biol. 2005, 15, 453–463. (14) Chockalingam, K.; Blenner, M.; Banta, S. Protein Eng. Des. Sel. 2007, 20, 155–161. (15) Colombo, G.; Soto, P.; Gazit, E. Trends Biotechnol. 2007, 25, 211– 218. (16) Tycko, R. Curr. Opin. Struct. Biol. 2004, 14, 96–103. (17) Tycko, R. Q. ReV. Biophys. 2006, 39, 1–55. (18) Chiti, F.; Stefani, M.; Taddei, N.; Ramponi, G.; Dobson, C. M. Nature 2003, 424, 805–808. (19) Tartaglia, G. G.; Cavalli, A.; Pellarin, R.; Caflisch, A. Protein Sci. 2004, 13, 1939. (20) Fernandez-Escamilla, A. M.; Rousseau, F.; Schymkowitz, J.; Serrano, L. Nat. Biotechnol. 2004, 22, 1302–1306. (21) Caflisch, A. Curr. Opin. Chem. Biol. 2006, 10, 437–444. (22) DuBay, K. F.; Pawar, A. P.; Chiti, F.; Zurdo, J.; Dobson, C. M.; Vendruscolo, M. J. Mol. Biol. 2004, 341, 1317–1326. (23) Pawar, A. P.; DuBay, K. F.; Zurdo, J.; Chiti, F.; Vendruscolo, M.; Dobson, C. M. J. Mol. Biol. 2005, 350, 379–392. (24) Tartaglia, G. G.; Pawar, A. P.; Campioni, S.; Dobson, C. M.; Chiti, F.; Vendruscolo, M. J. Mol. Biol. 2008, 380, 425–436. (25) Tartaglia, G. G.; Cavalli, A.; Pellarin, R.; Caflisch, A. Protein Sci. 2005, 14, 2723–2734. (26) Thompson, M. J.; Sievers, S. A.; Karanicolas, J.; Ivanova, M. I.; Baker, D.; Eisenberg, D. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 4074– 4078. (27) Conchillo-Sole, O.; de Groot, N.; Aviles, F.; Vendrell, J.; Daura, X.; Ventura, S. BMC Bioinf. 2007, 8, 65. (28) Trovato, A.; Seno, F.; Tosatto, S. C. E. Protein Eng., Des. Sel. 2007, 20, 521–523. (29) Ma, B.; Nussinov, R. Curr. Opin. Chem. Biol. 2006, 10, 445–452. (30) Khare, S. D.; Wilcox, K. C.; Gong, P.; Dokholyan, N. V. Proteins 2005, 61, 617–632. (31) Cecchini, M.; Curcio, R.; Pappalardo, M.; Melki, R.; Caflisch, A. J. Mol. Biol. 2006, 357, 1306–1321. (32) Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 1999, 314, 141–151. (33) Zanuy, D.; Ma, B.; Nussinov, R. Biophys. J. 2003, 84, 1884–1894. (34) Gsponer, J.; Haberthur, U.; Caflisch, A. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 5154–5159. (35) Cecchini, M.; Rao, F.; Seeber, M.; Caflisch, A. J. Chem. Phys. 2004, 121, 10748. (36) Guo, J.-t.; Wetzel, R.; Xu, Y. Proteins 2004, 57, 357–364. (37) Tsai, H. H.; Reches, M.; Tsai, C. J.; Gunasekaran, K.; Gazit, E.; Nussinov, R. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 8174–8179. (38) Baumketner, A.; Shea, J. E. Biophys. J. 2005, 89, 1493–1503. (39) Lo´pez de la Paz, M.; de Mori, G. M. S.; Serrano, L.; Colombo, G. J. Mol. Biol. 2005, 349, 583–596. (40) Ro¨hrig, U. F.; Laio, A.; Tantalo, N.; Parrinello, M.; Petronzio, R. Biophys. J. 2006, 91, 3217–3229. (41) Fogolari, F.; Corazza, A.; Viglino, P.; Zuccato, P.; Pieri, L.; Faccioli, P.; Bellotti, V.; Esposito, G. Biophys. J. 2007, 92, 1673. (42) Priya, A.; Fateh, S. N.; Ulrich, H. E. H. J. Chem. Phys. 2008, 129, 195102. (43) Wolf, M. G.; Jongejan, J. A.; Laman, J. D.; de Leeuw, S. W. J. Phys. Chem. B 2008, 112, 13493–13498. (44) Li, D.-W.; Mohanty, S.; Irback, A.; Huo, S. PLoS Comput. Biol. 2008, 4, e1000238.

11908

J. Phys. Chem. B, Vol. 114, No. 36, 2010

(45) Bellesia, G.; Shea, J.-E. Biophys. J. 2009, 96, 875–886. (46) Takeda, T.; Klimov, D. K. Biophys. J. 2009, 96, 442–452. (47) Giugliarelli, G.; Micheletti, C.; Banavar, J. R.; Maritan, A. J. Chem. Phys. 2000, 113, 5072. (48) Smith, A. V.; Hall, C. K. J. Mol. Biol. 2001, 312, 187–202. (49) Dima, R. I.; Thirumalai, D. Protein Sci. 2002, 11, 1036. (50) Peng, S.; Ding, F.; Urbanc, B.; Buldyrev, S. V.; Cruz, L.; Stanley, H. E.; Dokholyan, N. V. Phys. ReV. E 2004, 69, 41908. (51) Friedel, M.; Shea, J. E. J. Chem. Phys. 2004, 120, 5809. (52) Wei, G.; Mousseau, N.; Derreumaux, P. Biophys. J. 2004, 87, 3648– 3656. (53) Nguyen, H. D.; Hall, C. K. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 16180–16185. (54) Cheung, J. K.; Truskett, T. M. Biophys. J. 2005, 89, 2372–2384. (55) Shen, V. K.; Cheung, J. K.; Errington, J. R.; Truskett, T. M. Biophys. J. 2006, 90, 1949–1960. (56) Cellmer, T.; Bratko, D.; Prausnitz, J. M.; Blanch, H. W. Trends Biotechnol. 2007, 25, 254–261. (57) Bellesia, G.; Shea, J. E. J. Chem. Phys. 2007, 126, 245104. (58) Jianing, Z.; Muthukumar, M. J. Chem. Phys. 2009, 130, 035102. (59) Shell, M. S.; Ritterson, R.; Dill, K. A. J. Phys. Chem. B 2008, 112, 6878–6886. (60) Shell, M. S.; Ozkan, S. B.; Voelz, V. A.; Wu, G. H. A.; Dill, K. Biophys. J. 2009, 96, 917–924. (61) Lin, E.; Shell, M. S. J. Chem. Theory Comput. 2009, 5, 2062– 2073. (62) Voelz, V. A.; Shell, M. S.; Dill, K. A. PLoS Comput. Biol. 2009, 5, e1000281. (63) Pearlman, D. A.; Case, D. A.; Caldwell, J. W.; Ross, W. S.; Cheatham Iii, T. E.; Debolt, S.; Ferguson, D.; Seibel, G.; Kollman, P. Comput. Phys. Commun. 1995, 91, 1–41. (64) Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M., Jr.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. J. Comput. Chem. 2005, 26, 1668–1688. (65) Kollman, P. A.; Dixon, R.; Cornell, W.; Fox, T.; Chipot, C.; Pohorille, A. Comput. Simul. Biomolec. Syst. 1997, 3, 83–96. (66) Onufriev, A.; Bashford, D.; Case, D. A. J. Phys. Chem. B 2000, 104, 3261–3429.

Lin and Shell (67) Ozkan, S. B.; Wu, G. H. A.; Chodera, J. D.; Dill, K. A. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 11987–11992. (68) Geney, R.; Layten, M.; Gomperts, R.; Hornak, V.; Simmerling, C. J. Chem. Theory Comput. 2006, 2, 115–127. (69) Kabsch, W.; Sander, C. Biopolymers 1983, 22, 2577–2637. (70) Shrake, A.; Rupley, J. A. J. Mol. Biol. 1973, 79, 351–364. (71) Ho, B. K.; Dill, K. A. PLoS Comput. Biol. 2006, 2, e27. (72) Kallberg, Y.; Gustafsson, M.; Persson, B.; Thyberg, J.; Johansson, J. J. Biol. Chem. 2001, 276, 12945–12950. (73) Ono, S.; Nakajima, N.; Higo, J.; Nakamura, H. J. Comput. Chem. 2000, 21, 748–762. (74) Garcı´a, A. E.; Sanbonmatsu, K. Y. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 2782–2787. (75) Gnanakaran, S.; Garcia, A. E. Biophys. J. 2003, 84, 1548–1562. (76) Mu, Y.; Kosov, D. S.; Stock, G. J. Phys. Chem. B 2003, 107, 5064– 5073. (77) Wood, S. J.; Wetzel, R.; Martin, J. D.; Hurle, M. R. Biochemistry 1995, 34, 724–730. (78) Tjernberg, L. O.; Callaway, D. J. E.; Tjernberg, A.; Hahne, S.; Lillieho¨o¨k, C.; Terenius, L.; Thyberg, J.; Nordstedt, C. J. Biol. Chem. 1999, 274, 12619–12625. (79) Petkova, A. T.; Ishii, Y.; Balbach, J. J.; Antzutkin, O. N.; Leapman, R. D.; Delaglio, F.; Tycko, R. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 16742– 16747. (80) Williams, A. D.; Portelius, E.; Kheterpal, I.; Guo, J.-t.; Cook, K. D.; Xu, Y.; Wetzel, R. J. Mol. Biol. 2004, 335, 833–842. (81) von Bergen, M.; Friedhoff, P.; Biernat, J.; Heberle, J.; Mandelkow, E. M.; Mandelkow, E. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 5129–5134. (82) Nelson, R.; Sawaya, M. R.; Balbirnie, M.; Madsen, A. O.; Riekel, C.; Grothe, R.; Eisenberg, D. Nature 2005, 435, 773–778. (83) Godawat, R.; Jamadagni, S. N.; Garde, S. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 15119–15124. (84) Lopez de la Paz, M.; Goldie, K.; Zurdo, J.; Lacroix, E.; Dobson, C. M.; Hoenger, A.; Serrano, L. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 16052–16057. (85) Pawar, A. P.; Dubay, K. F.; Zurdo, J.; Chiti, F.; Vendruscolo, M.; Dobson, C. M. J. Mol. Biol. 2005, 350, 379–392.

JP104114N