Computational Simulations of the Early Steps of Protein Aggregation

Abstract

There is strong evidence that the oligomers of key proteins, formed during the early steps of aggregation, could be the primary toxic species associated with human neurodegenerative diseases, such as Alzheimer’s and prion diseases. Here, we review recent progress in the development of computational approaches in order to understand the structures, dynamics and free energy surfaces of oligomers. We also discuss possible research directions for the coming years.

Introduction

More than 20 human diseases are associated with the aggregation of key proteins. For instance, the 40- and 42-residue amyloid β-protein (Aβ),1 99-residue microglobulin (β2m) protein,2 and 210-residue PrP prion protein,3 are linked to Alzheimer's disease, dialysis-related amyloidosis, and prion diseases, respectively. All these proteins have different amino acid sequences and biophysical properties in solution: Aβ is random coil,4 β2m has an all-β topology while PrP is essentially α-helical,5 however, they all misfold and aggregate into amyloid fibrils with a cross-β structure, characterized by β-sheets perpendicular to the fibril axis.6–10

Although the amyloid plaques are one of the hallmarks of the diseases, there is strong evidence that early oligomers during the aggregation process may be the primary toxic species.11–15 The initial assemblies of oligomers are difficult to characterize at an atomic level of detail using biophysical methods, because they are transient and in dynamic equilibrium between dimers, trimers, tetramers etc. The experimental signature of the rate-limiting step is the presence of a lag phase in polymer growth, which varies with protein concentration, salt and agitation. Once a nucleus is present, maturation into fibrils is rapid.16–18 It follows that we have very limited information on both the secondary and tertiary structures of these early oligomers, although an α-helix signal has been detected by circular dichroism,19,20 and based on antibody experiments,21 they may share a common topology, different from the fibrils, independently of the amino acid sequence.

In this context, computer simulations of protein aggregation have provided insights into the early steps of amyloid formation. Such theoretical investigation is made possible by the fact that peptides of four to seven amino acids can form fibrils indistinguishable from those formed by proteins. For example, a reptation mechanism first identified numerically,22,23 seems to be dominant in the rearrangement of amyloid chains.24 These results, and others presented below, show that simulations can complement and even guide experimental measurements in the study of amyloid aggregation. In this review, we limit ourselves to coarse-grained protein aggregation simulations taken largely from our work. Excellent reviews on the use of all-atom molecular dynamics25 or discontinuous molecular dynamics simulations26 to understand protein aggregation or bioinformatics27 to predict regions promoting amyloid fibril formation can be found elsewhere.

In the first part, we report the progress in the development of computational approaches to monitor amyloid fibril formation with their strengths and weaknesses. Then, we present the main results that have been extracted from coarse-grained protein simulations in order to understand the structures, dynamics and free energy surfaces of oligomers. Finally, we list three directions that others and we will probably follow in the field.

In Silico Approaches to Simulate Amyloid Aggregation

Molecular dynamics (MD), which integrates directly Newton's equation of motion, offers the most detailed picture at the atomic level, providing both dynamical and thermodynamic information. Because of the integration 2 fs timestep, all-atom MD in explicit solvent is limited to trajectories of 1 µs—typically about 100 ns. Such a time scale might be sufficient to study the stability of preformed structures28–34 such as the cross-beta-spine steric zipper of the Sup-35 prion fragment,7,35 the very early events in the dynamics of Aβ36–40 or the docking of unstructured monomer on preformed structured oligomers,41 but other methods of various degrees of efficiency and accuracy are needed to span the aggregation regime from monomers to fibrils, which requires several days in vitro.

A first approach is the replica exchange procedure, initially proposed by Swendsen and Wang42 and then reformulated with a molecular dynamical scheme by Sugita and Okamoto.43 The replica exchange molecular dynamics (REMD) mixes a series of MD trajectories (or replica) run in parallel at different temperatures through Monte Carlo accept/reject moves. The probability of exchanging two conformations i and j run at T_h and T_l, respectively, is given by min {1,exp[(β_h - β_l)(E_j - E_i)]} where β is 1/kT. Because REMD allows conformations to move between various T, one expects a more extensive sampling of the low-energy structures than with standard MD, providing a better description of the thermodynamic properties of the system, at the cost of the loss of dynamical information. However, even at high temperature, the decorrelation time can be longer than the simulations in slow systems. This is the case, for example, when explicit solvent is used or when large assemblies of full-length proteins are studied.44 It follows that this limitation directly affects the quality of the results. Another severe limitation of REMD is that the number of replicas increases rapidly with the size of the system, limiting its applicability to small all-atom systems in solvent such as the monomer of full-length Aβ, and dimers and trimers of short fragments.45–48

The use of implicit solvent and coarse-grained models lifts some of the limitations of both MD and REMD on two counts. First, by removing some of the dampening due to the collision with the solvent, they accelerate the sampling of the phase space by as much as two orders of magnitude.49 Second, the simplification of the potential decreases the computational cost allowing longer simulations.50–52 For example, Jang and Shin could generate multiple 150-ns MD trajectories for a trimer of Aβ(10–35) using a Born-generalized implicit solvent.50 Similarly, Paci et al could reach a total time of 0.5 µs on tetramers of TTR(105–115).53 The use of implicit solvent and coarse-grained models do not affect all motions uniformly, however. This is not important for thermodynamics, but these modifications can affect the details of the aggregation dynamics.54

Implicit solvent protein models also allow the use of a wider range of methods. Activated approaches, such as the activation-relaxation technique (ART nouveau) for example,55,56 explores the space of conformations, through well-defined transition states. Coupled with the implicit solvent coarse-grained protein OPEP force field,57,58 this method has provided the best fit with the NMR data of the fragment Aβ(21–30) in solution59 and has allowed to monitor the aggregation of various amyloid peptides with 4 to 15 amino-acids, in settings ranging from dimer to dodecamer.22,23,60–66 Because it lacks detailed balance ART cannot provide solid thermodynamic information. Trajectories are physically-based, however, and are qualitatively correct, based on comparisons with other approaches using a 16-residue b-hairpin model.67

Other approaches have also been used to characterize the first steps of aggregation. Irbäck et al uses a simplified off-lattice potential with a Monte Carlo approach based on two elementary backbone moves. This reduction in motion speeds up the simulations, making it possible to follow aggregation of hexamers of Aβ(16–22).68 This is also the case of discrete molecular dynamics (DMD), which requires a much-simplified force field with square-well interactions in order to evolve the time based on collisions. Such an approach would be unworkable in the presence of explicit solvent, but with implicit solvent DMD generates trajectories corresponding to seconds or more, treating up to 100 chains or so.69–71 It remains to be determined whether the very simple force field used provides the correct dynamics or even the proper thermodynamics. More characterization on a wide range of sequences and structures is clearly necessary at this moment.

Free Energy Surfaces of Dimers

It is well established that the early steps of amyloid-fibril formation are characterized by the formation of low molecular weight oligomers consisting of a mixture of dimers, trimers, tetramers, and more in rapid equilibrium.17,18 Because Aβ exists as a stable dimer in solution at low concentrations,72,73 and dimers may act as seeds for larger oligomers, there has been considerable theoretical efforts to characterize the free energy surface and the dynamics of dimer formation using various protein models.22,44,74–78

Here, we probe the free energy surface of the Aβ(16–22) dimer, resulting from a 50-ns REMD-OPEP79 simulation in a sphere of 70 Å diameter, starting from two chains in random orientation. The integration timestep is 1 fs and T is controlled by the Berendsen's bath80 with a coupling constant of 0.1 ps. We use eight replicas with T varying between 287 and 500 K with exponential distribution and an exchange time between neighboring replicas of 20 ps, leading to an acceptance ratio between 30–40%.

Figure 1 shows the free energy surface of Aβ(16–22) dimer at 310 K. The two reaction coordinates are the cosine of the angle between the two chains and the extended status of the chains. The extended status is the product of the end-to-end distance of the chains divided by the product of the end-to-end distance for two ideal β-strands. We see multiple free energy minima. These are in-register and out-of-register parallel strands (structures A, B and D), parallel chains (structure C), cross chains (structure E), antiparallel loops (structure F), and antiparallel strands (structures G and H). The Boltzmann probabilities of structures A–H at 310 K are: 1, 1, 6, 5, 3, 21, 10 and 12%. The calculated percentage of β-strand content is 36% at 310 K. Taken together, these results indicate that the antiparallel arrangement of the chains is preferred over the parallel arrangement, in agreement with previous reports.22,77 The in-register and out-of-register antiparallel β-structures (structures G and H) are similar in free energies, helping explain the experimental dependency of β-sheet registry on pH conditions.81 It also follows that this free energy surface generated by REMD-OPEP is very similar in character to that generated by all-atom explicit solvent REMD simulations.44 This finding is interesting for two reasons. First, it shows that cross conformations (structure E, Fig. 1) are also populated using a coarse-grained protein model with implicit solvent if thermal fluctuations are considered. Second, it opens the door to the study of the free energy surfaces for many peptides in dimers or higher-order species using reasonable computer resources. We find, for example, that the free energy surface of the human β2m(83–89) dimer is very similar to that shown in Figure 1, although the 7-residue sequences differ: KLVFFAE for Aβ(16–22) vs. NHVTLSQ for β2m(83–89).

Structures and Dynamics of Transient Oligomeric Species

While simulations of dimers and trimers in solution provide insights into the populated structures,23,41,49,48 larger systems are needed to better understand the early events of aggregation and the dynamics leading to ordered structures. Experimentally, we often observe amorphous and annular aggregates prior to the formation of protofibrils.82,83

Using ART and MD simulations coupled to OPEP, we study the aggregation of a number of peptides of increasing length (from 4 to 15 amino acids) and oligomer size (from trimer to dodecamer), starting from random orientations of the chains and random coil conformations of each chain. While we observe a high percentage of ordered β-sheet structures for trimers and tetramers of Aβ(16–22) and KFFE,23,63 tetramers and heptamers of β2m(83–89),84,85 and tetramers to octamers of KFFE,60,61,63,64 self-assembly to fibrillar structures is more problematic for a tetramer of Aβ(11–25)62 and a dodecamer of NFGAIL.66 Figure 2 shows a generic aggregation picture leading to fibrillar states, emerging from all these simulations. Starting from an initial state with all chains randomly placed, the peptides first come together to form amorphous aggregates with two-stranded or three-stranded β-sheet rapidly in place. Then, the outcome of both ART and MD simulations using OPEP varies according to the chain length and the oligomer size.

Small oligomer size and small chain length.

For a system containing less than 8 chains of small lengths, the amorphous aggregates evolve either directly to fully ordered structures (orthogonal or parallel β-sheets) or indirectly through closed or open β-barrels. The generated parallel ordered structures display the cross-β characteristics observed experimentally,7 with Cα..Cα distances of 5.0 Å between the strands and around 10.0 Å between the layers. In contrast, orthogonal β-sheets only display the meridional 5.0 Å reflection. Both structures have been observed by MC simulations of six Aβ(16–22) chains68 and DMD simulations of polyalanines69 and prion fragments.70 Similarly, orthogonal layers have been observed by all-atom MD simulations starting from parallel or antiparallel layers.31,63,64 The amorphous aggregates, the bilayer β-sheet and the β-barrel structures are all in dynamic equilibrium.64,65 Using MD-OPEP on seven β2m(83–89) chains, we find that the transition is more rapid from β-barrel to β-sheets than from β-sheets to β-barrel, and the estimated time scale for both reactions is on the order of the ms range in explicit solvent.85 Of course, this time will increase with the oligomeric size, and the shape of the β-barrel will also vary with the oligomer size and the length of the peptides. Seven or eight chains of KFFE can only form open β-barrel,61,64 while seven chains of β2m(83–89) are sufficient to stabilize into a closed β-barrel with hydrogen-bonding interactions between all chains.

We identify two important mechanisms in all ART simulations, and more recently in MD simulations of four Aβ(16–22) chains.52 The first one is the reptation move of the chains. This motion allows the chains to rearrange their H-bond networks without having to fully detach, decreasing significantly the free-energy cost of realignment.22,52,62 This move has recently been observed experimentally in oligomers formed by Aβ(16–22) using isotope-edited IR spectroscopy.24 The second mechanism we find refers to as the two-stage dock-lock mechanism for adding a monomer on a preformed oligomer.41,66 Both mechanisms are likely general processes for error-correction in the early steps of protein aggregation, but also in amyloid growth by monomer addition.60,66,86

Large oligomer size or long chain length.

When the number of chains is twelve or more, or when the chain length increases to 15 amino acids, aggregation to ordered β-sheet structures is very rare, and the peptides form amorphous aggregates that are in dynamic equilibrium. Among a total of 12 ART simulations of twelve NFGAIL chains, one run locates highly ordered β-sheet structure within 12,000 events.66 Similarly, multiple simulations of four Aβ(11–25) chains fail to explore a four-stranded β-sheet within 30,000 events.62 In both cases, the elongation process is difficult, as the amorphous oligomer is dense, even when detachment/reattachment of the peptides and reptations moves of the chains are considered. This implies that, just like in the case of monomeric protein folding, there are many kinetic traps and the acquisition of a β-sheet oligomeric structure requires minimal frustration and a low conformational entropy.66

Future Directions

Simulating protein aggregation is a challenging problem that pushes the limits of current methods. Progress will likely be achieved by combining these various methods to take advantages of their particular strengths while minimizing their weaknesses. Thus far, most simulations of oligomers have been performed in solution at various pH and concentration conditions. While they have already provided insights into self-assembly pathways and structures of oligomers, three open questions remain to be addressed.

First, it is essential that these simulations are repeated in more realistic cellular environment, with full treatment of metallic ions (especially copper) and membrane. Interactions with copper87 and membrane37 have already been investigated by short MD simulations, but they only explore local fluctuations around the starting structures. In the case of Aβ, we may go one step beyond and incorporate the effect of cholesterol, and even apolipoprotein E and its allele A4, but the latter system is problematic because there is no structure available.88

Second, there have been many reports suggesting constrained pathways leading to fibrils, but very few solid verifications. Several theoretical studies have suggested intermediate states with α-helix89 and even α-sheet character.90 Similarly, an α-helical signal has been detected by circular dichroism in the late steps of Aβ aggregation, but there is no evidence that these intermediates are on-pathways.17 Rather, our simulations on Aβ(16–22), NFGAIL, Aβ(11–25) and β2m(83–89)23,62,66,85 suggest several pathways for self-assembly, in agreement with experimental data, and transient sampling of species including amorphous aggregates, α-helical intermediates with a very low population, and β-barrels. This β-barrel differs from the disorganized annular prototype seen by DMD91 and microscopy measurements.83 While we do not know if the β-barrel is accessible to full-length Aβ, its structural characteristics makes it an ideal system to create pores within the membranes, which in turn could contribute to the enhanced toxicity of the oligomeric intermediates.83

Third, there has been no theoretical attempt to target the interaction sites between current inhibitors and oligomers at an atomic level of detail. These compounds reported to block amyloid aggregation are based either on D-peptides, N-methylated peptides, peptides containing proline substitutions.92–94 For instance, Meredith et al find that the N-methylated Aβ(16–22) peptide at positions 17, 19 and 21 inhibits the fibrillogenesis of full-length Aβ.95 Similarly, Yang et al showed that curcumin inhibits formation of Aβ oligomers and fibrils, and reduces amyloid in vivo.96 This is a very difficult numerical problem, but free energy surfaces of oligomers-inhibitors using coarse-grained models should tell us where inhibitors bind, as this would help rational drug design.

Figures and Tables

Figure 1 Free energy surface of Aβ(16–22) dimer at 310 K obtained from REMD-OPEP simulation. The two reaction coordinates used are the cosine of the angle between the two chains and the extended status of the two chains. The structures A–H of low free energy (in kcal/mol) are shown.

Figure 2 A generic aggregation picture derived from ART- and MD-OPEP simulations. Starting from a randomly chosen state, the peptides form amorphous aggregates. From there, the outcome changes with the oligomer size (OS) and chain length (L). For OS < 9 and L < 8, rapid aggregation proceeds directly to ordered β-sheets or indirectly through β-barrels. The double arrows indicate reversibility. For larger OS or L, aggregation into β-barrels and ordered β-sheets is very slow and rare.

Acknowledgements

G.W. thanks the National Natural Science Foundation of China under Grant No.10674029, the Young Foundation of Fudan University, the Senior Visiting Scholar Grant of Fudan University. N.M. acknowledges partial funding from the NSERC, FQRNT and the Canada Research Chair Program. P.D. acknowledges financial support from CNRS, University Paris 7 and the 6^th European PCRD. N.M. and P.D. also thank financial support from the Alzheimer Society of Canada.