Scenarios for protein aggregation: Molecular Dynamics simulations and Bioinformatic Analysis

The need to understand the assembly kinetics of fibril formation has become urgent because of the realization that soluble oligomers of amyloidogenic peptides may be even more neurotoxic than the end product, namely, the amyloid fibrils. In order to fully understand the routes to fibril formation one has to characterize the major species in the assembly pathways. The characterization of the energetics and dynamics of oligomers (dimers, trimers etc) is difficult using experiments alone because they undergo large conformational fluctuations. In this context, carefully planned molecular dynamics simulation studies, computations using coarse-grained models, and bioinformatic analysis have given considerable insights into the early events in the route to fibril formation. Here, we describe progress along this direction using examples taken largely from our own work. In this chapter, we focus on aspects of protein aggregation using Abeta-peptides and prion proteins as examples.

antiparallel β-sheet structure in the fibril state, is an ideal system for exploring in detail the mechanism of oligomer formation. Using fairly long and multiple trajectories [6], the assembly pathways for 3 Aβ 16−22 → (Aβ 16−22 ) 3 were probed using all atom simulations in explicit water.
The simulations of Aβ 16−22 and the corresponding mutants allowed us to draw a number of conclusions that may be of general validity.
The simulations of the Aβ 16−22 monomer at room temperature and at neutral pH showed that it is predominantly a random coil. The finite size of the system gives rise to large conformational fluctuations that lead to the population of strand-like structures. There is a very low (∼ 3%) probability of α-helical conformations. The study of this simple system shows that the β-sheet conformation adopted by the monomer must be due to interactions with other peptides.
Oligomerization of three Aβ 16 These studies and other simulations illustrate the complexity in dissecting the assembly of even small amyloidogenic peptides into ordered structures.
The gross features of the fibril structures of a number of proteins and peptides whose monomer sequences and structures are unrelated are similar. This observation might suggest that the interactions that stabilize the oligomers and fibrils must be "universal" involving perhaps only backbone hydrogen bonds. It might appear that side chains, and hence sequence differences, might play a secondary role. Such a conclusion is further supported by repeated observations [21] that any protein or peptide can be made to form cross-β structures under appropriate conditions. However, experiments [22] and simulations [7,23] show that side chain interactions are crucial in directing oligomer formation. Trimers of Aβ 16−22 are stabilized primarily by favorable interpeptide hydrophobic interactions between residues in the central hydrophobic cluster (LVFFA) and secondary by inter-peptide salt bridge formation between K and E.
The importance of side chains can be demonstrated by examining the effect of mutations on the trimer formation in Aβ 16−22 peptides. We showed using simulations that the mutant GLVFFAK, which eliminates the formation of intermolecular salt-bridge entirely destabilizes the trimer. Similarly replacement of L, F, F by S, also destabilizes the trimer. These simulations show that the sequence plays a key role in the tendency of peptides to form amyloid fibrils.
Although no general role has emerged it seems that sequences with enhanced correlation between charges [12] or preponderance of contiguous (> 3) hydrophobic residues might be amenable to amyloid formation on finite time scales.

Dimerization of Aβ 10−35 peptides:
Generation of putative dimer structures: In contrast to Aβ 16−22 fibrils the longer peptide Aβ 10−35 adopts a parallel β-sheet conformation in the amyloid state. It is now suspected that in the fibril state the monomer is stabilized by an intra-molecular salt bridge between Asp23 and Lys28. In order for this salt bridge to form there has to be a bend in the monomeric structures involving the residues VGSN. The importance of a stable turn, which was experimentally determined in the NMR structures, was emphasized in MD simulations as well.
In a recent study [8], we used a number of computational methods to probe dimer formation.
We first generated a putative set of dimer conformations that is based on shape complementarity.
The work on Aβ 16−22 showed that both inter-peptide hydrophobic interactions and the creation of favorable electrostatic contacts are required to produce marginally stable oligomers. In order to dissect their relative importance we generated two homodimer decoy sets by maximizing the number of contacts between the monomer interfaces. The first 2000 dimer structures of each set were selected by minimizing the interaction energy between the monomers. In order to distinguish between desolvation and electrostatic interactions we used two distinct energy functions. The ϕ-dimer (Fig. (2(a))) minimizes the desolvation energies of the dimer at the interface whereas the ε-dimer ( Fig. (2(b))) corresponds to structures that have the highest inter-peptide electrostatic interactions. The structure of the ϕ-dimer is dominated by contacts between hydrophobic segments of the monomers. The hydrophobic core, LVFFA (17)(18)(19)(20)(21), and the hydrophobic C-terminus of both monomers are buried at the dimer interface. The contacts at the interface of the ϕ-dimer are conserved over the lowest energy dimer structures. The ε-dimer interface is characterized by electrostatic inter-monomeric interactions, among which the salt bridge Glu11(A)-Lys28(B) has the largest contribution. Contrary to the ϕ-dimer, the contacts observed at the ε-dimer interface are not conserved across the set of the low energy dimers due to the increased specificity and strength of the electrostatic interaction.
Interior of Aβ oligomers is dry: Insights into the assembly mechanism of the ϕ-dimer and ε-dimer can be obtained from the Potential of Mean Force (PMF). The PMF for the dimerization process was obtained along the center of mass of the two monomers as the Fig.(3) indicates. For each free energy profile, one can distinguish three distinct intervals. In the outer interval, the PMF value is nearly constant, from 6.5Å -7.0Å to maximum separation, which in our case is 9.0Å. At a distance of 6.5Å for the ε-dimer and 7.0Å for the ϕ-dimer, the first solvation shells of the monomers come into contact, and for both dimers the energetics of desolvation of the associating monomers is unfavorable. In the second interval for the ε-dimer, the value of the PMF continues to increase up to 1.2 kcal/mol at a 3.0Å separation; for the ϕ-dimer, the potential energy reaches a value of 0.8 kcal/mol at 5.5Å, and after that the desolvation is favorable, ending in an unstable local minimum at 3.0Å. For the third interval, from 3.0Å to 0.0Å, there is only one solvation shell between the monomers. The water molecules are most strongly ordered near the monomers through electrostatic interactions and hydrogen bonds. As a result, the PMF for the ε-dimer increases sharply between 3.0Å and 1.3Å up to 2.4 kcal/mol. At contact, the van der Waals attraction predominates, making the overall dimerization process energetically favorable. For the ϕ-dimer, the solvation shell between the hydrophobic regions of the monomers is only weakly bound to the solute, and after a small increase in the PMF, corresponding to the van der Waals attraction, the desolvation is entirely favorable.
If the approach along the center of mass of the monomers approximately represents a minimum energy path, then the expulsion of water in the ϕ-dimer must be an early event in the assembly. Explicit simulations for Aβ 16−22 oligomers [6] also show that desolvation occurs early. As a result, the interior of Aβ oligomers is dry.
Hydrophobic interactions between monomers are the driving force in the association of Aβ 10−35 peptides into dimers: Comparing the ϕ-dimer ε-dimer models for monomer association, we find that the former appears to lead to more energetically favorable dimerization than the latter. It appears to be more efficient to remove the entropically unfavorable structured water between the opposing hydrophobic regions of the two monomers than to stabilize the monomer solely through electrostatic interactions. This is in good agreement with the experimental and MD simulations observation that the mutation E22Q -where a charged glutamic acid residue is replaced by a polar glutamine residue -increases the propensity for amyloid formation [24,25].
Molecular dynamics simulations of this increased amyloidogenic activity for the E22Q mutant peptide led to the conclusion that the water-peptide interaction is less favorable for the mutant peptide [26]. Following a more detailed analysis of the structure and dynamics of the WT and E22Q Aβ 10−35 , it has suggested that a change in the charge state of the peptide, due to the E22Q mutation, leads to an increase of the hydrophobicity of the peptide that could be responsible for the increased activity [27].
The time evolution of the ϕ-dimer structure was analyzed and it was observed that the monomers remain in contact during the simulation. It was shown that the hydrophobic interaction between the monomers of the ϕ-dimer acts as a stabilizing force of the dimer. The "extended core" region 15-30 of both monomers in the ϕ-dimer makes the principal contribution to the hydrophobic interaction energy. The ϕ-dimer undergoes internal structural reorganization in the terminal regions of the monomeric peptides. Our simulations indicate that there is substan-tial reorganization of the peptide monomers in the N-and C-terminus regions, as expected for a dimer weakly and relatively non-specifically stabilized by hydrophobic contacts at the dimer interface. Importantly, the structure of the central hydrophobic cluster LVFFA region assumes a conformation similar to that observed for the monomeric peptide in both experiment [28] and simulation [29]. Our simulations suggest that the preservation of the structure of the LVFFA central hydrophobic cluster plays an important role in the stabilization of the ϕ-dimer structure.
The finding that the ϕ-dimer may constitute the ensemble of stable Aβ 10−35 dimer has important implications for fibril formation. The initial event in the dimerization involves, in all likelihood, contacts between the central hydrophobic clusters. In this process, expulsion of water molecules in the interface might be a key event just as in the oligomerization of Aβ 16−22 fragments [6]. Since this process involves cooperative rearrangement of ordered water molecules, it is limited by an effective free energy barrier. Based on our results, we conjecture that events prior to the nucleation process themselves might involve crossing free energy barriers which depend on the peptide-peptide and peptide-water interactions (Fig. (1)).

Initial stages in the PrP C conformational transition:
Prion proteins are extracellular globular proteins which are attached to the plasma membrane by a GPI anchor. They have been linked to various transmissible spongiform encephalopathies (TSEs) including the bovine spongiform encephalopathy, the scrapie disease in sheep, and the Creutzfeldt-Jakob disease in humans. The causative agent in these diseases is believed to be the aggregated form (PrP Sc ) of the cellular prion protein (PrP C ) [30]. The transition to the scrapie form involves a large conformational change from the mainly α-helical PrP C to the PrP Sc state that is rich in β-sheet. According to the "protein-only" hypothesis [2] PrP Sc serves as a template in inducing conformational transitions in PrP C that can subsequently be added to PrP Sc . The "protein-only" hypothesis implies that the conformational change leading to the PrP Sc formation from the normal cellular form PrP C may be spontaneous or might involve interactions with unidentified protein X [31]. Prion proteins, encoded by a single gene, consist of about 250 residues of which the first 22 form a signal sequence. This is followed by unstructured, but likely helical, Cu 2+ binding octarepeats rich in glycine [2]. The NMR [30,32,33] and X-ray [34] structure of PrP C in various species (human, mouse, syrian hamster, bovine, and sheep) shows that the ordered C-terminal part is composed of a short antiparallel β-sheet that contains 8% of the residues in the (90-231) fragment and three helices representing 48% of the secondary structure ( Fig.(4). Fourier transform infrared spectroscopy measurements [35,36] indicate that PrP Sc (90-231) has 47% β-sheet and 24% α-helical content.
We have suggested using structural, bioinformatic, and molecular dynamics simulations that formation of PrP Sc follows scenario II (see Fig. (1)). This implies that, either spontaneously or in the presence of a seed of PrP Sc , the metastable cellular form, PrP C , undergoes a transition to the PrP C * state that is capable of further aggregating or adding to an already present PrP Sc particle. Experiments [37] and scenarios of protein aggregation [14] suggest the proposal that the conformational transition involves the formation of PrP C * that is more stable than PrP C .
The transition from the metastable PrP C → PrP C * , which involves crossing a substantial free energy barrier on the order of 20 kcal/mole [17,38], results in a state that can nucleate and polymerize to the protease resistant form. We also identified the putative regions that are involved in the P rP C → P rP C * transition. Comparison of a number of structural characteristics (such as solvent accessible area, distribution of (Φ, Ψ) angles, mismatches in hydrogen bonds, nature of residues in local and non-local contacts, distribution of regular densities of amino acids, clustering of hydrophobic and hydrophilic residues in helices) between PrP C structures and a databank of "normal" proteins shows that the most unusual features are found in helix 2 (residues 172-194) followed by helix 1 (residues 144-153) [11]. In particular, the C-terminal residues in H2 are frustrated in their helical state. Application of the recently introduced notion of discordance, namely, incompatibility of the predicted and observed secondary structures, also points to the frustration of H2 not only in the wild type but also in mutants of human PrP C .
This suggests that the instability of PrP C proteins may play a role in their being susceptible to the profound conformational change.
We showed [11] that, in addition to the previously proposed role for the segment (90-120) and possibly H1, the C-terminus of H2 and possibly N-terminus may play a role in the α → β transition. Sequence alignments show that helices in avian prion proteins (chicken, duck, crane) are better accommodated in a helical state which might explain the absence of PrP Sc formation over finite time scales in these species. From the analysis it is clear that the conformational fluctuations in the C-terminal end of helix 2 (H2) and in parts of helix 3 (H3) are involved in the transition to PrP C * . Because the stability of PrP C arises from the structures in the C-terminal end, the transition to PrP C * requires global unfolding of PrP C [39] which explains the origin of the high free energy barrier separating PrP C and PrP C * [11]. NMR experiments [37,40] showed that conformational fluctuations that originate in the C-terminal part of H2 are essential in the formation of PrP C * . Structural and mutational studies have also shown that the relatively short helix 1 (H1) is stable over a range of pH values and solvent conditions, and hence is unlikely to undergo conformational change in the transition to PrP C * [41,42,43].
The required conformational fluctuations in PrP C needed to populate PrP C * suggest that the earliest event involves extensive unfolding of the monomeric PrP C . We used results from a database search of sequence patterns in helices of PrP C and extensive all atom molecular dynamics (MD) simulations of helical fragments from the mouse prion protein (mPrP C ) to shed light on the nature of instabilities that drive the PrP C → PrP C * transition [44].
Previously MD simulations have been used to probe other structural aspects of prion proteins including structures of protofibrils [45]. The 10 residue H1, with an unusual sequence pattern Pattern of hydrophobicity in H2 is rare: There are very few sequences that share the pattern of hydrophobicity of H2. In PDBAstral40 [46] (proteins in the PDB having at most 40% sequence similarity) there are only 12 (0.2%) such sequences. In the E Coli genome the number is 46 (1%), while in the Yeast genome it is 122 (1.4%). Inspection of the structures of the 12 proteins from PDBAstral40, shows that the sequence is never entirely helical! For example, in only 13% of these proteins the last 5 residues are found in a helix. A characteristic pattern seen in H2 from mammalian prion proteins is TTTT (positions 190-193). In the PDBAstral40 this pattern occurs in only 18 proteins, including the prion sequence. In an overwhelming number of these cases (15 of the 18 proteins) the TTTT pattern is found in a strand and/or loop conformation (irrespective of the identity of the flanking amino acids). These results add further support to our proposal [11] that the second half of H2 would be better accommodated in non-helical conformations.
"Frustrated" secondary structural elements may be harbinger of tendency to polymerize The ease of aggregation and the morphology of the aggregates depend not only on the protein concentration, but also on other external conditions such as temperature, pH and salt concentration.
Although most proteins can, under suitable conditions, aggregate the observation that several disease causing proteins form amyloid fibrils under physiologically relevant conditions raises the question: Is aggregation or the need to avoid unproductive pathways encoded in the primary sequence itself ? It is clear that sequences that contain a patch of hydrophobic residues are prone to form aggregates [47]. However, it is known that contiguous patches (three or more hydrophobic residues) occur with low probability in globular proteins [48]. For example, sequences with five hydrophobic residues (LVFFA in Aβ peptide) in a row are not well represented. Similarly, it is unusual to find hydrophobic residues concentrated in a specific region of helices such as is found in helix 2 in PrP C [11].
It is natural to wonder if secondary structure elements bear signatures that could reveal amyloidogenic tendencies. Two studies have proposed that the extent of "frustration" in the secondary structure elements (SSE) may be harbinger of amyloid fibril formation [11,49].
Because reliable secondary structure prediction requires knowing the context dependent propensities and multiple sequence alignments (such as used in PHD, Profile network from Heidelberg [50]), it is more likely that assessing the extent of frustration in the SSE rather than analysis of sequence patterns is a better predictor of fibril formation. Frustration in SSE is defined as the incompatibility of the predicted (from PHD, for example) secondary structure and the experimentally determined structure [49]. For example, if a secondary structure is predicted to be in a β-strand with high confidence and if that segment is found (by NMR or X-ray crystallography) to be in a helix, then the structure is frustrated (or discordant or mismatched). The α/β discordance, which can be correlated with amyloid formation, can be assessed using the score S α/β = 1 L L i=1 (R i − 5), where R i is the reliability score predicted by PHD at position i of the query sequence, 5 is the mean score, and L is the sequence length.
The bounds on S α/β are 0 ≤ S α/β ≤ 4 with maximal frustration corresponding to S α/β = 4. Similarly, the measure S β/α gives the extent of frustration of a stretch that is predicted to be helical and is found experimentally to be a strand. Using Proposed structures for PrP C * : Our simulations [44] and recent experiments [42,43] strongly suggest that H1 is unlikely to change conformation in the PrP C → PrP C * transition. The most drastic change occurs in the second half of H2 and parts of H3. Based on the assumption that alterations in the conformation of H2+H3 do not significantly affect the rest of the protein, we have constructed a plausible ensemble of structures for PrP C * (Fig.(6(b)) and 6(c)). In PrP C * (90-231) obtained from the NAMD trajectories ( Fig.(6(b))) the helical content is ∼ 20% (a lower bound), and in PrP C * (90-231) reached during the MOIL simulations (Fig.6(c)) the helical content is ∼ 30% compared to 48% in mPrP C (90-231).
The overall characteristics of these structures are consistent with those proposed by James and coworkers [37]. It remains to be seen if formation of PrP C * , with fluctuating regions in H2+H3, is required for oligomerization of PrP C i.e., if PrP C * is an on-pathway monomeric intermediate on the route to fibrillization. We should emphasize that the conformation of the prion protein in PrP Sc need not coincide with PrP C * .
Comparison of PrP C * structure with the Human prion protein dimer: In an important paper Knaus et al. [61] announced a 2Åcrystal structure of the dimeric form of the human prion protein (residues 90-231). The structure (Fig.(7)) suggest that dimerization occurs by domain swap mechanisms in which H3 from one monomer packs against H2 from another. In fact, Eisenberg and coworkers have suggested that domain swapping mechanism may be a general route of amyloid fibril formation [62]. The electron density map seems to suggest structural fluctuations in the residues 189-198 which coincides with the maximally frustrated region predicted theoretically. The dimer interface is stabilized by residues that are in H2 in the monomeric NMR structures. The header of the PDB file of the monomeric structure of human PrP C indicates that H2 ends at residue 194 and H3 begins at 200. The domain swapped dimer structure shows that residues 190-198 exist largely in a β-strand conformation. The α → β transition minimizes frustration. An implication of the dimer structure is that oligomerization occurs by domain swapping which in PrP C also might implicate the disulfide bond between Cys residues at 179 and 214. The role of the S-S bond in the PrP Sc formation remains controversial. In our full-atom MD simulations of reduced mouse PrP C (data not published), we found structures that closely resemble the monomeric structure in this dimer. For example, H1 remains mostly intact, while H2 breaks into two smaller helices, one running from its normal N-term end to position 187 and the other being formed by the C-term end residues of the original H2 and residues from the loop connecting it to H3. These findings suggest that the dimer structure is a likely route to unfolding and self-assembly of monomeric PrP C .
Sequence pattern matches and long multiple molecular dynamics simulations of helix 1 in mPrP C using two force fields show that the stability of H1 is due to the formation of stabilizing internal salt bridges. In view of the high propensity of α-helix observed in the isolated H1 in conjunction with supporting experimental results [41,42,43] it is clear that alterations in the conformation of H1 are unlikely in the PrP C → PrP C * transition.
The predicted tendency for the second half of H2 to be involved in the formation of

Conclusions
The development of methods to envision the structure of amyloid fibrils has enabled us to obtain molecular insights into the assembly process itself. Computational and experimental studies are beginning to provide detailed information, at the residue level, about the regions in a given protein that harbor amyloidogenic tendencies. We have harnessed these developments to propose tentative ideas on the molecular basis of protein aggregation. These principles (or, more precisely, rules of thumb) may be useful in the interpretation and design of new experiments.
Examination of the stable structures of oligomers and fibrils obtained using experimental restraints and simulations show that these must be stable conformations which maximize the inter-peptide interactions and minimize electrostatic repulsions. Broadly, this is the only amyloid self-organization principle (ASOP) that seems to be obeyed. From the ASOP it follows that the formation of amyloid fibrils should indeed be a generic property of almost all proteins and peptides under suitable conditions. If this were the case then it is remarkable that in the case of normal function aggregation is avoided most of the time. The lack of preponderant protein aggregation may well be due to the efficiency of cleaning mechanisms operating in the cell. This may explain the lack of aggregation of PrP C * under most circumstances. We conjecture that because of efficient degradation processes only mild sequence constraints are needed to prevent oligomer formation during the typical life cycle of newly synthesized proteins.
From a biophysical perspective there are a number of open problems. Are there common pathways involved in the self-assembly of fibrils? Because of the paucity of the structural description of the intermediates involved in an aggregation process a definitive answer cannot be currently provided. The energy landscape perspective, summarized briefly in Fig.(1), suggests that multiple scenarios for assembly must exist. Although the generic nucleation and growth governs fibril formation the details can vary considerably. The microscopic basis for the formation of distinct strains in mammalian prions and in yeast prions remains a mystery. Are these merely associated with the heterogeneous seeds or are there unidentified mechanisms that lead to their growth? What factors may determine the variations in the fibrillization kinetics for the wild type and the mutants? A tentative proposal is that the kinetics of polymerization is determined by the rate of production of N ⋆ (Fig. (1)) [63], which in turn is controlled by barriers separating N and N ⋆ [11,64]. In this scenario the stability of N plays a secondary role.
The generality of this observation has not yet been established. Finally, how can one design better therapeutic agents based on enhanced knowledge of the assembly mechanism? Even in the case of sickle cell disease viable therapies began to emerge only long after the biophysical aspects of gelation were understood [65]. landscape perspective. According to scenario I, the assembly competent state N ⋆ is metastable with respect to the monomeric native state N and is formed due to partial unfolding. In scenario II N ⋆ is formed upon structural conversion either of the native state N (as in prions) or directly from the unfolded state U (as in Aβ-amyloid peptides). In both cases proteins (or peptides) in N ⋆ states must coalescence into larger oligomers capable of growth into fibrils.  computed using the ε-dimer as the starting structure. The red curve is similarly computed using the ϕ-dimer as the starting structure. The difference between the two surfaces suggests that hydrophobic interactions may be more essential to stabilization of the dimer structure than electrostatic interactions.   PrP C is taken from the PDB file 1ag2 (yellow). The conformations for PrP C * contain H1 from 1ag2 while the residues encompassing H2+H3 are shown in a conformation (red) reached towards the end of our MD simulations using the NAMD package (6(b)) or the simulations using the MOIL package (6(c)). The schematic PrP C * structures are representatives from ensembles of fluctuating conformations. In the representative PrP C * structure obtained using NAMD simulations the H1 region, together with the adjacent loops and the β-strands, and residues (205-212) from H3 retain their original conformations and are therefore depicted with same color as in PrP C . In the MOIL representative PrP C * structure the H1 region, together with the adjacent loops and the β-strands, and residues (175-179), (184-188), (193,194) from H2 and residues (203-218) from H3 retain their original conformations and are therefore depicted with same color as in PrP C . The figures are rotated such that the orientation of H1 is the same in all of them. The figures were produced with packages VMD [66] and PovRay (http://www.povray.org/). Fig.7. Cartoon representation of the X-ray crystal structure of the human PrP C dimer (PDB entry 1I4M). For each chain, A and B, the three helices in the 90-231 ordered region of PrP C are colored in red, while the short β-sheet is in yellow. The two cysteine residues (179 and 214) involved in the disulfide bond that connect H2 with H3 are indicated in bond representation and colored in purple. The C-term end of H2 and the N-term end of H3 which we believe to be implicated in the initial stages of the α → β transition are colored in green. We notice that, in contrast to the monomeric PrP C structure from Fig.(4), here this region is no longer entirely helical, but contains a short stretch of β-strand structure and a shorter helix as well. The figure was produced with packages VMD [66] and PovRay (http://www.povray.org/).