Thermodynamics of Amyloid-β Fibril Elongation: Atomistic Details of the Transition State.

Amyloid-β (Aβ) fibrils and plaques are one of the hallmarks of Alzheimer's disease. While the kinetics of fibrillar growth of Aβ have been extensively studied, several vital questions remain. In particular, the atomistic origins of the Arrhenius barrier observed in experiments have not been elucidated. Employing the familiar thermodynamic integration method, we have directly simulated the dissociation of an Aβ(15-40) (D23N mutant) peptide from the surface of a filament along its most probable path (MPP) using all-atom molecular dynamics. This allows for a direct calculation of the free energy profile along the MPP, revealing a multipeak energetic barrier between the free peptide state and the aggregated state. By definition of the MPP, this simulated unbinding process represents the reverse of the physical elongation pathway, allowing us to draw biophysically relevant conclusions from the simulation data. Analyzing the detailed atomistic interactions along the MPP, we identify the atomistic origins of these peaks as resulting from the dock-lock mechanism of filament elongation. Careful analysis of the dynamics of filament elongation could prove key to the development of novel therapeutic strategies for amyloid-related diseases.

Introduction

The presence of abnormal fibrillar aggregates of amyloid species in the brain is associated with the development of neurodegenerative diseases such as Alzheimer’s (AD), Parkinson’s, and Huntington’s.[1−3] Amyloid-beta (Aβ) peptides in Alzheimer’s are produced from the cleavage of amyloid-β precursor protein (AβPP) via sequential cleavage by β and γ-secretases. While the exact pathogenic role of amyloid-β aggregating into senile plaques is still not well-understood, fibrillar associations of Aβ have emerged as key neurotoxic species in Alzheimer’s disease affected brain areas such as the hippocampus, cerebral cortex, and amygdala.

The aggregation and elongation processes of amyloid-β fibrils have been previously studied in kinetics experiments.[4−11] On the basis of such, a two-step process for fibril elongation has been proposed, commonly known as the dock-lock mechanism. In this paradigm, first a free Aβ peptide rapidly adheres (docking) to a mature preformed fibril, followed by a slower second step where a rearrangement of the incoming peptide occurs into the ordered fibrillar shape (locking).[5,6,10] From kinetics experiments, this process is known to follow the Arrhenius equation with an associated energetic barrier. Detailed molecular analysis of this binding process of Aβ fibrils is crucial for the better understanding of protein aggregation processes into neurotoxic species in AD, identification of key residues/domains involved in fibrillation of Aβ, and the design of potential inhibitors of fibrillar growth.

Molecular dynamics (MD) simulations have proven a valuable and precise tool to elucidate details in complement to experimental studies on Aβ kinetics. The dynamics of fibril aggregation have been simulated at various levels of resolution by coarse-grained,[12−17] implicit water,[18,19] and hybrid models.[20] All-atom simulations have typically been restricted to short peptides,[21−23] although some longer peptides have been studied recently on this level.[24−26] Recent work has also studied the free energy landscape as a function of peptide number.[14] The long time scales, elucidated by kinetics experiments, involved in the aggregation process present a major hurdle for atomistic simulations of fibril elongation. In addition, the atomistic origins of the Arrhenius barrier observed in experiments[4−11] have not been explained in the current literature. Such all-atom simulations are necessary for a detailed molecular picture of the elongation process and for the determination of the relevant thermodynamic parameters.

In the present study, we have simulated the elongation process of Aβ using all-atom molecular dynamics. By conducting (slow) steered molecular dynamics (SMD) experiments coupled with the thermodynamic integration method, we simulated the process of dissociating a monomeric peptide from a mature amyloid-β filament. By relaxing the dissociated peptide at each step of the unbinding pathway (similar to the gentlest ascent algorithm[27]), we have directly sampled the most probable path (MPP) of fibril growth under the imposed SMD constraints. Evaluating the Gibbs free energy (the potential of mean force,[28] PMF) along the MPP yielded a multipeak activation barrier consistent with the kinetics experiments of the current literature. Analysis of the detailed atomistic interactions along the MPP allowed us to identify the origins of the Arrhenius barriers as related to the dock-lock mechanism of Aβ fibril elongation.

Results and Discussion

We have simulated an Aβ peptide unbinding event from an amyloid filament by steering a monomer along its MPP. We computed the PMF along the MPP, revealing a multipeak energetic barrier separating the dissociated state from the aggregated state. In what follows, we elucidate the origins of this barrier via atomistic analyses of the Aβ peptide-filament system along the unbinding pathway.

Atomistic Origins of the Energetic Barrier

The key findings of the present work are shown in Figure . Our calculations indicate the presence of a multipeak energetic barrier along the unbinding pathway, showing at least three characteristic peaks along the reaction coordinate (hereafter referred to as RC or as P0s displacement interchangeably). The estimated Arrhenius barrier was 8.7 ± 1.4 kcal/mol above the dissociated state (average of three sets of simulations, uncertainty is the largest standard deviation). The overall energetic barrier obtained from the simulations of the system along the Aβ filament elongation pathway is in agreement with the experimental results of Hasegawa and co-workers, who estimated the energetic barrier to be approximately 10 kcal/mol for wild-type Aβ40.[29] Recently, Young et al.[30] measured the elongation barrier for the wild-type Aβ42 to be 11.2 kcal/mol. Thus, our result of 8.7 kcal/mol for the D23N Iowa mutant isoform of Aβ40 agrees with the existent experimental measurements.

Figure 1

Top: Free energy of elongation along the most probable path (MPP) with dock-lock steps identified. Bottom: Snapshots of the dissociation process corresponding to the three different peaks in the free energy profile. Note how during the unlocking steps the N- and C-termini of P0 and P1 are still tightly interacting, whereas during the undocking step they break away from each other. All molecular graphics in this work were rendered with VMD.[38]

Free energy profiles along the elongation pathway have been computed previously. Schwierz et al.[24] computed a monotonic profile involving only a one-dimensional PMF, in contrast to our 18-dimensional profile, which could help to explain why their PMF does not exhibit a barrier. Our computational strategy allows us to directly probe the most probable path, which yields the Arrhenius barriers along the elongation pathway, in agreement with experimental reports. The free energy computed by Gurry and Stultz[19] exhibits a barrier of approximately 4 kcal/mol in the locking regime, most likely due to their choice to restrain the coordinates of the topmost layer (P1) of the aggregate. Our simulations indicate that the dynamically fluctuating hydrophobic interactions between P0 and P1 are critical for the elongation process of Aβ fibrils.

An additional parameter of interest is the elongation rate (k+) for monomer addition. Our computed barrier allows us to estimate k+ aswhere the diffusion rate for monomer addition has been estimated[24] to be 3 × 109 M–1 s–1, Ea is the Arrhenius barrier, kB is the Boltzmann constant, and T is the temperature (298 K in our simulations). This yields an estimate for k+ of 1.2 × 103 M–1 s–1, in contrast with Young et al.’s[30] experimental value of 9.3 × 105 M–1 s–1 for Aβ42. The discrepancy between these values is due to our neglect of the entropic contributions to the elongation rate (see the computational details in the Methods and in the Supporting Information).

To elucidate the atomistic basis for the barrier, we analyzed the main hydrophobic contacts between the steered peptide (P0) and the nearest strand in the aggregate (P1) (Figure ). Figure shows a snapshot along the steering path. In what follows, we use filament/aggregate interchangeably to refer to an organized layer of Aβ peptides (P1–P8 in Figure ). As evidenced by the proximity between P0 and P1 during the first two displacement peaks, these states may be assigned to the locking step of filament elongation kinetics.[5,6,10] In contrast, the definite separation during the third displacement peak is consistent with the docking step. Thus, our simulations are in agreement with the experimental results as well as with the template-dependent dock-lock mechanism of Aβ elongation.

Figure 2

Snapshot of the steering pathway with the steered alpha carbons shown as magenta-colored spheres. The numbering used to refer to the peptides constituting the filament is indicated.

Hydrophobic Interactions Promote Assembly of Aβ Filaments

Which residues of Aβ are essential for the transition from a free peptide to the well-packed elongated filament? To answer this, we studied the hydrophobic interactions between P0 and P1. We monitored the α carbon (Cα) distances between the corresponding hydrophobic residues in P0 and P1, (e.g., the distance between P0-Val36-Cα and P1-Val36-Cα) as a function of P0s displacement. Clusters of hydrophobic residues have been suggested to stabilize Aβ filament structure.[31] In addition, mutations of specific residues such as Glu22 into hydrophobic residues have been found to increase the aggregation of Aβ as a function of the mutation’s hydrophobicity.[14,32] We identified three distinct elongation stages involving three sets of hydrophobic residues. The first set of residues included Ala30, Ile31, and Ile32, which initially showed stable Cα–Cα distances until RC reached a value of 5 Å (Figure ), consistent with the first peak in the free energy profile. The second peak in the free energy profile can be identified with residues Ala21, Ala30, and Ile32 (Figure ). Ala21 presented a particularly abrupt increment in the Cα–Cα distance at RC ∼ 9 Å, in conjunction with transitions in the Cα–Cα distances of Ala30 and Ile32, which had already contributed to the first peak (Figure ). The third peak (Figure ), corresponding to the undock step resulted from the overall contribution of hydrophobic amino acids in the C-terminus of Aβ: Leu34, Met35, Val36, Val39, and Val40. Our results indicated that the Cα distances of this third group transitioned abruptly at RC around 12–13 Å, which we attribute to a breakage of the associated hydrophobic contacts among the residues that provide stability and promote the well-organized assembly of Aβ filaments. Of this last set, Leu34, Met35, and Val36 had been identified as being highly flexible by solid state 2H NMR line shape experiments.[33] In particular, this NMR study identified Val36 as having the highest conformational variability for the Iowa mutant, which agrees with our results. Snapshots of the three elongation stages are shown in Figure , bottom. Trajectory movies showing snapshots centered around these peaks and highlighting the residues involved are available in the Supporting Information as Movies 1, 2, and 3.

Figure 3

Residues involved in the first peak of the PMF. Top left: Alpha carbons of Ala30, Ile31, and Ile32 represented as magenta-colored spheres. P0 and P1 are shown as ribbons colored by residue type (white for hydrophobic, green for hydrophilic, blue for positively charged, and red for negatively charged), while the rest of the peptides are shown in cartoon representation colored by residue type. Top right, bottom left and right: Distance between corresponding alpha carbons of P0 and P1 as a function of the reaction coordinate. Arrows highlight the transition, occurring at displacements consistent with the first peak.

Figure 4

Residues involved in the second peak of the PMF. Top left: Alpha carbons of Ala21, Ala30, Ile32 represented as magenta-colored spheres. P0 and P1 are shown as ribbons colored by residue type (white for hydrophobic, green for hydrophilic, blue for positively charged, and red for negatively charged) while the rest of the peptides are shown in cartoon representation colored by residue type. Top right, bottom left and right: Distance between corresponding alpha carbons of P0 and P1 as a function of the reaction coordinate. Arrows highlight the transition, occurring at displacements consistent with the second peak.

Figure 5

Residues involved in the third peak of the PMF. Top left: Alpha carbons of Leu34, Met35, Val36, Val39, and Val40 represented as magenta-colored spheres. P0 and P1 are shown as ribbons colored by residue type (white for hydrophobic, green for hydrophilic, blue for positively charged, and red for negatively charged), while the rest of the peptides are shown in cartoon representation colored by residue type. Top right, bottom left and right: Distance between corresponding alpha carbons of P0 and P1 as a function of the reaction coordinate. Arrows highlight the transition, occurring at displacements consistent with the third peak.

Recently, Mason et al.[34] studied the self-assembly of phenylalanine dimers and found that they exhibit remarkably similar activation parameters compared to full-length Aβ42. The importance of these residues in the assembly of Aβ filaments has been recognized in multiple studies.[35,36] Our metric of measuring the corresponding Cα distances between hydrophobic residues of P0 and P1 did not reveal the same abrupt transitions for Phe19, Phe20 as it did for the residues in Figures –5. However, the aromatic interactions between these residues are still significant for the peptide–aggregate interactions. In Figure S4 in the Supporting Information, we show the Cα–Cα distances between P0 and P1 for every pair of hydrophobic residues of the Aβ15–40 atomic model as a function of the reaction coordinate.

Our simulations indicate that the C-terminal hydrophobic residues are responsible for the initial docking stage. In contrast, previous studies[19,20,37] have identified the central hydrophobic core as the initial point of attachment. In the study by Gurry et al.[19] and the work by Han et al.,[20] both groups positionally restrained the atoms of the core filament. We, on the other hand, only fixed the alpha carbons on P4 (Figure ), which allows us to capture the dynamics between the incoming monomer (P0) and the topmost layer of the aggregate (P1). Recently, Bacci et al.[37] used advanced sampling techniques to probe the monomer-filament interactions that give rise to the slow locking pathways. Our findings are in complement to theirs since our study focused on the Arrhenius barrier for elongation rather than the complex conformational rearrangement that accompanies the slower aggregation steps of Aβ. In addition, our findings on the importance of the dynamics of P1 for elongation are consistent with the study by Bacci et al., which also illustrated the role of the penultimate peptide in the elongation process.[37]

Release of Water Molecules from Filaments

Aβ is generated by endoproteolysis of AβPP, a Type I transmembrane protein of 695–770 amino acids. The fragment corresponding to Aβ is part of the transmembrane domain, originally located in a hydrophobic environment of the cell membrane. Upon release of Aβ, it can exist in multiple assembly states (monomer, oligomers, protofilaments and filaments). The filament assembly process is complex, nucleation-dependent and the mechanism driving this process is not well-understood.[39] In the aggregate state, the P0–P1 interface is essentially devoid of water molecules.[23,24] Our simulations along the unbinding pathway showed that a hydration layer begins to form around P0 and over the exposed surface of the aggregate. Figure shows the number of water molecules (N) within 5 Å of P0 or the aggregate as a function of the reaction coordinate. A trajectory movie showing the formation of a hydration layer as the peptide transverses the binding pathway is available in the Supporting Information as Movie 4. We highlight that these water molecules are located outside the inner filament. Interestingly, a recent report[37] suggests that protofibril-internal waters do not contribute to the elongation process. In the neurodegenerative context, the opposite phenomenon occurs, where the incoming or newly formed Aβ monomers or dimers are fully hydrated (in the cerebrospinal fluid) and will achieve a dry interface state directly related to the approach distance to the filament template. This assembly process is repeated until the formation of mature fibers, diffuse plaques and neuritic plaques in the affected areas of the brain. Hydrophobic residues of Aβ mediate specific interactions that direct the self-assembly and are sufficient to promote Aβ aggregation as confirmed through site specific mutation variants of Aβ.[40] Taken together, these results confirm the importance of hydrophobic residues in the aggregation mechanism.

Figure 6

Waters surrounding the Aβ-peptide–filament system as a function of the reaction coordinate. Waters were considered to surround the system if their oxygen atoms were within 5 Å of either P0 or the filament. Error bars are standard deviations from the three sets used to calculate the free energy profiles in Figure .

Thermodynamics of Elongation

Additional insight into the origin of the energetic barriers of Aβ elongation is possible from analyses of the system’s enthalpy. For a better understanding of the individual thermodynamic contributions, we decomposed the total enthalpy into various components: solvent–P0 (SP), aggregate–P0 (AP), P0–P0 (PP), aggregate–aggregate (AA), solvent–solvent (SS), and solvent–aggregate (SA). These individual contributions can be further dissected into bonded and nonbonded interactions, the latter of which consist of electrostatics and van der Waals (vdW) interactions. The relevant hydrophobic contacts identified in Figures –5 are primarily based on vdW interactions, which are distance-dependent and are weakened as P0 undergoes the corresponding unlocking/undocking steps. As P0 is progressively exposed to the solvent, we expect the vdW interactions between the solvent and P0 to increase. Figure shows the total, the vdW, and the electrostatic interactions between the solvent and P0 as well as between the aggregate and P0 as a function of P0s displacement. Note how the displacement at which the vdW plot rapidly increases (for SP) is consistent with the transition states of Figure . Similar decompositions for the peptide–peptide, aggregate–aggregate, solvent–aggregate, and solvent–solvent subsystems are shown in Figures S2 and S3 in the Supporting Information. It should be noted that, while the computation of the free energy profile is robust, the decomposed interactions are much less so. This is, in part, a consequence of the limited sample space molecular dynamics simulations can probe. Figure thus reveals qualitatively consistent trends with our robust PMF results in Figure , but have limited quantitatively predictive power.

Figure 7

Energetics between the solvent and the peptide (SP) and the aggregate and the peptide (AP) along the binding pathway. Left column: Total, vdW, and electrostatic interaction along the binding pathway for SP. Right column: Total, vdW, and electrostatic interaction along the binding pathway for AP. The rapid transitions and subsequent relaxations in the vdW plots are consistent with the free energy profile of Figure . In particular, the subtle jumps at displacements ∼9 and ∼12 Å in the AP vdW plot are consistent with those of Figures and 4.

The elongation steps present in our results involve dynamical fluctuations between the incoming monomer (P0) and the topmost layer of the filament (P1). In particular, P1 undergoes significant structural changes in order to “welcome” P0 onto the preformed filament. Dynamics similar to these have been described previously, e.g., Bacci et al.[37] identified multiple putative docking states associated with increased disorder at the filament tip. The consistency between our decomposed enthalpies and the fluctuations between P0 and P1, quantified by their hydrophobic Cα-Cα distances (Figures –5), serves to corroborate our proposed atomistic mechanism of filament elongation.

In summation, we simulated an elongation event for the D23N “Iowa” mutant of Aβ15–40 by steering six α carbons on the odd end of the filament. By equilibrating the system under the steering constraints, we determined the overall free energy of the system along the 18-dimensional (6 × 3 degrees of freedom) most probable binding path. Our results show a free energy barrier (8.7 kcal/mol) consistent with experimental results showing that the elongation of Aβ is an activated process.

Along the elongation pathway, we found three distinguishable peaks in the free energy profile. Our analyses indicate that the origins of these peaks follow from the relevant hydrophobic interactions between the steered Aβ peptide and the template filament. The aggregation mechanism showed significant similarities to the dock-lock mechanism of filament elongation. Initially, a set of hydrophobic contacts from the C-terminus of the Aβ filament “docks” the incoming peptide to the preformed filament template, after which additional interactions “lock” it in place. This process also releases water molecules that were part of the hydration shell of the Aβ peptide and the exposed filament surface.

Analyzing the energetics of our system reveals consistent trends with the hydrophobic interactions we identified along the elongation pathway. By decomposing the total enthalpy of the system into pair contributions and further dissecting these into bonded and nonbonded interactions, it was possible to identify distinct transitions in the van der Waals interaction energies between the steered Aβ peptide and the solvent, as well as between the steered Aβ peptide and the Aβ filament at reaction coordinate values of ∼5, 9, and 12 Å, consistent with the transition states in the free energy profile. The specific identification and classification of individual contributions of hydrophobic amino acids to the self-assembly process of Aβ peptides provides valuable information on the most critical and reactive sites triggering the formation/stabilization of Aβ filaments. This information may open opportunities for the design of novel diagnostic or therapeutic compounds that specifically target these active subdomains in Aβ.

Methods

System Preparation and Equilibrium Dynamics

The starting point for our simulations was the NMR structure for the in-register, parallel, D23N Iowa mutant of Aβ15–40 filaments (PDB ID: 2MPZ) determined by Sgourakis et al.[41] We extracted a single filament/aggregate (Figure ) from the trimeric structure (chains A, D, G, J, M, P, S, V, Y) and fully solvated it by adding water molecules (TIP3P[42]) in a box with dimensions of 80 × 80 × 100 Å3. Then, we added counterions to neutralize the system and set the NaCl concentration to 150 mM, mimicking the physiological conditions of the cerebrospinal fluid (ionic strength and pH). We used NAMD 2.12[43] in conjunction with the CHARMM36[44] force field for the equilibrium simulations. We optimized the structure in the isobaric–isothermal (NPT) ensemble for 10 ps after which we equilibrated the system for 22 ns. We used the Nosé–Hoover barostat with a piston target of 1 bar. We also used Langevin dynamics with a friction coefficient γ of 5/ps for the thermal control. For electrostatics calculations, we used periodic boundary conditions with the particle-mesh Ewald method (with grid points of 128 × 128 × 128). We used a time step of 1 fs for bonded interactions, 2 fs for short-range nonbonded interactions, and 4 fs for long-range nonbonded interactions.

Sampling the Transition States of Aβ

Starting from the final state of the equilibrium simulations, we steered the alpha carbons of residues Gln15, Phe20, Gly25, Asn27, Gly33, and Val40 (the collective positions of which will be denoted as R = {r15,r20,r25,r27,r33,r40}, where r denotes the position vector of residue j’s α carbon) of the odd tip (peptide P0) of the Aβ filament as shown in Figure . We steered the alpha carbons (R) along the filament axis at a rate of 2.5 Å/ns. We chose as our reaction coordinate the center of mass displacement of R along the filament axis (that is, if we define the center of mass position rcm as , then the reaction coordinate is the component of rcm along the filament axis) with respect to the initial configuration. Each steering segment was followed by 4 ns of relaxation in the degrees of freedom orthogonal to the reaction coordinate (by disallowing fluctuations of the steered alpha carbons). This is similar to the gentlest ascent algorithm of Crippen and Scheraga.[27] Restricting the movement of the alpha carbons is also necessary to ensure each subsequent steering section is continuous with the previous one, that is, the configuration of the steered carbons at the end of section j is identical to the one at the beginning of section j + 1. We divided the dissociation pathway into 115 sections, the first 85 of which were 0.2 Å in length, while the remaining covered 0.4 Å each. The choice of atoms to steer is guided by how well they represent the overall position and orientation of the peptide. Different choices may lead to convergence issues for the free energy profile. Additional details regarding the PMF computation can be found in the Supporting Information.

To calculate the Gibbs free energy (ΔG or the PMF) as a function of P0s displacement, we employed thermodynamic integration. First, we equilibrated the resulting states at the end of each steering section for an additional 2 ns (for a total initial equilibration of 6 ns because of the initial 4 ns between each steering segment). Second, we conducted three independent 1 ns long samplings for the mean force acting on the steered atoms, where each data set was used as the input for the next one, i.e., the first set was sampled for 6 + 1 ns (discarding 6 ns), the second one for 6 + 1 + 1 ns (discarding 7 ns) and the third set was sampled for 6 + 1 + 1 + 1 ns (discarding 8 ns). This was done to improve the statistics of our results as well as to test the convergence of the method used. The positions of the steered alpha carbons were not allowed to fluctuate during this sampling, but the forces acting on them were recorded during the course of the simulations. By integrating over the mean force on the steered atoms, we are effectively calculating the thermodynamic integral along the 18-dimensional curve defined by the displacement of the 3 × 6 degrees of freedom of the steered atoms in the phase space. If the steering speed is low enough and/or the system has been equilibrated enough, the steering pathway will be near the MPP of the physically feasible unbinding pathway (the reverse of which represents the physical binding pathway). In addition, the free energy profile calculated from each set will converge to a common value, as seen in Figure . Two nonconverged trials, with initial equilibrations at each steering section of 4 and 5 ns, respectively, are shown in Figure S1 in the Supporting Information.