Amyloid aggregates of human islet amyloid polypeptide (hIAPP or human amylin) have long been implicated in the development of type II diabetes. While hIAPP is known to aggregate into amyloid fibrils, it is the early-stage prefibrillar species that have been proposed to be cytotoxic. A detailed picture of the early-stage aggregation process and relevant intermediates would be valuable in the development of effective therapeutics. Here, we use atomistic molecular dynamics simulations with a combination of enhanced sampling methods to examine the formation of the hIAPP dimer in water. Bias-exchange metadynamics calculations reveal relative conformational stabilities of the hIAPP dimer. Finite temperature string method calculations identify pathways for dimer formation, along with relevant free energy barriers and intermediate structures. We show that the initial stages of dimerization involve crossing a substantial free energy barrier to form an intermediate structure exhibiting transient β-sheet character, before proceeding to form an entropically stabilized dimer structure.
I. INTRODUCTION
Amyloidogenic proteins have long been implicated in a host of human diseases. These proteins exhibit a shared tendency to self-assemble into aggregates known as amyloid, which share morphological properties, including a fibrillar shape with heavily β-sheet secondary structure.1 One such amyloidogenic protein is human islet amyloid polypeptide (hIAPP or human amylin); this 37-residue polypeptide hormone is co-secreted with insulin in the pancreas by β-cells and plays a role in regulating blood glucose levels.2,3 Heavily β-sheet amyloid aggregates of hIAPP have been pathologically linked to the loss of pancreatic β-cells and the onset of type II diabetes,4 prompting studies of hIAPP fibrils and their formation.
Extensive characterization efforts have probed the structure of the mature hIAPP fibril, including the use of solid-state NMR (ssNMR) experiments to propose the arrangement of individual hIAPP chains within a mature fibril. In this proposed model, hIAPP monomers are stacked along the fibril axis and each monomer is arranged in a U-shaped conformation consisting of two β-strands (residues 8-17 and 28-37), connected via a loop region. Parallel β-sheets are formed as a monomeric hIAPP stack to form the mature fibril structure. This model is consistent with two-dimensional infrared spectroscopy (2D-IR) experiments,5 as well as electron paramagnetic resonance (EPR) measurements.6 Furthermore, recent experiments and simulations have studied amylin in the presence of inhibitors,7–9 binding of amylin to metals,10 interactions between amylin and its mutants or other amyloidogenic proteins,11,12 as well as amylin behavior at a membrane8,13–17 and structural rearrangements during aggregation.18
Mature fibrils are relatively biologically inert, and the formation of prefibrillar species, or protofilaments, has been linked to cytotoxicity,19–22 prompting a shift toward the study of early-stage amyloid aggregates. Dimers, trimers, and larger oligomers have been proposed to be responsible for disrupting cell membranes, inhibiting metabolic functions, inducing oxidative stress, and triggering apoptosis.5,20,23–25 Experimental evidence includes observations of membrane leakage prior to mature fibril formation26 and of disrupted cell membranes isolated from areas of fibril growth.27 In order to better understand the role of hIAPP in disease onset, as well as to design effective therapeutics, it is crucial to uncover the mechanisms of early-stage fibril formation, as well as identify any relevant intermediate structures in the aggregation process.
In particular, residues 20-29 in hIAPP are suspected to play a key role in amyloid formation.28 This 10-residue segment is itself amyloidogenic, and mutations to this sequence suppress amylin aggregation.28–31 However, residues 20-29 are not located within the heavily β-sheet fibril core in the ssNMR-based fibril model.1 2D-IR experiments and molecular dynamics (MD) simulations have previously been combined to identify a key intermediate in the early-stage aggregation of hIAPP, which featured the FGAIL region in residues 23-27 in a transient parallel β-sheet prior to the formation of the final U-shaped configuration.32 In addition to this FGAIL region, residues L12A13 have recently been suggested to form a stacked turn or disordered β-sheet during the aggregation through 2D-IR experiments with dihedral indexing.33 Furthermore, Chiu and de Pablo have utilized MD simulations combined with bias-exchange metadynamics (BE-MetaD) to study the mechanism by which two disordered hIAPP monomers assemble into a U-shaped dimer,32,34 and more recent work has extended this approach to investigate dimer formation in the presence of a lipid membrane.17 Two potential dimerization mechanisms were proposed, both of which exhibit an intermediate parallel β-sheet structure in residues 23-27 (FGAIL region). The disordered dimer was found to be less thermodynamically stable than both the intermediate and the final U-shaped dimer; however, the work concludes with a cautionary note that these results are highly dependent on the force field used (GROMOS96 53a6).34,35 Subsequent force field reviews suggested that GROMOS96 53a6 tends to over-stabilize the formation of β-sheet in both human and rat islet amyloid polypeptide (rIAPP),36 as well as the amyloidogenic polypeptide polyglutamine.37 Newer force fields were shown to more accurately capture experimentally determined properties, for example Cα secondary NMR chemical shifts: one such force field is AMBERff99SB*-ILDN.38,39
Although GROMOS96 53a6 does predict an intermediate structure consistent with previous experimental and computational studies, given the findings of Hoffmann et al.,36 further investigation of amylin aggregation warrants the use of a force field that better captures the conformational behavior of the hIAPP molecule, especially if studies are to be extended to higher order aggregates. In this work, we seek a more accurate understanding of early-stage hIAPP aggregation in water using the AMBERff99SB*-ILDN force field. We use bias-exchange metadynamics to reveal the free energy landscape for hIAPP dimerization under the new, more accurate force field. Furthermore, we go beyond past studies of the dimer and employ finite temperature string method calculations to unveil specific aggregation mechanisms for the hIAPP dimer, as well as the associated changes in free energy, revealing a substantial free energy barrier associated with the formation of an intermediate β-sheet structure prior to the formation of a locally stable structured amylin dimer.
II. RESULTS AND DISCUSSION
A. Free energy landscape for the hIAPP dimer
Bias-exchange metadynamics (as described in Sec. IV) using the AMBERff99SB*-ILDN force field was used to produce the free energy landscape of the hIAPP dimer in solution [Fig. 1(a)], plotted as a function of two collective variables (CVs), Qres 8-16 and Qres 27-35. The two Q parameters measure similarity of residues 8-16 and 27-35 in the simulated structure to the ideal U-shaped dimer, which is extracted from the ssNMR mature fibril structure. Further details on the chosen collective variables are found in Sec. IV. The global free energy minimum is found at (Qres 8-16, Qres 27-35) = (0.43, 0.38). Two local minima corresponding to fully dimerized structures are found at the upper right of the plot at (0.88, 0.88) and (0.88, 0.68), with free energies of 7.4 kJ/mol and 6.8 kJ/mol, respectively. Additional local minima are found at (0.83, 0.33), with a free energy of 6.5 kJ/mol, and (0.23, 0.68), with a free energy of 7.3 kJ/mol; these minima correspond to partially dimerized structures.
Several features of the free energy landscape can be pieced together to understand the dimerization process. First, the global minimum is found in a single, wide basin in the lower left quadrant of the free energy landscape. In Fig. 1(a), conformations within 5 kJ/mol of the global minimum are bounded by the second contour line out from the global minimum, shown in yellow. Structures in this lower left quadrant correspond to low values of Qres 8-16 and Qres 27-35; these conformations show little similarity to the ideal U-shaped dimer. The free energy minimum that most closely matches the U-shaped dimer is found in the upper right quadrant at (0.88, 0.88), elevated over the global minimum by 7.4 kJ/mol, or approximately 3 kBT at room temperature. The relative depths, widths, and locations in CV space of these two minima suggest that the conformational ensemble of the hIAPP dimer is largely dominated by disordered structures dissimilar from the U-shaped dimer. Additionally, the free energy minima in the upper right quadrant at (0.88, 0.88) and (0.88, 0.68) are more elongated along the Qres 27-35 axis versus the Qres 8-16 axis, suggesting that C-termini β-strands are more flexible than the N-termini β-strands. Furthermore, the conformational terrain is relatively smooth outside of the global minimum, with no local free energy minima deeper than 2.5 kJ/mol below their immediate surroundings; the major barrier to the dimerization process is escaping the wide well surrounding the global minimum centered at (0.43, 0.38).
A number of qualitative differences arise between the AMBERff99SB*-ILDN and GROMOS96 53a6 models. The GROMOS model previously suggested that the U-shaped dimer is found at the global free energy minimum, at 1.1 kJ/mol more stable than the disordered dimer. The GROMOS model also identified multiple metastable states, with free energies within kBT of the global minimum; each of these metastable states featured a significant amount of β-sheet structure. In comparison, the free energy landscape obtained using AMBERff99SB*-ILDN features a U-shaped dimer at 7.4 kJ/mol above the disordered dimer in free energy, with a single wide basin centered around the global minimum containing disordered conformations. This discrepancy is consistent with the findings of Hoffmann et al. that found that GROMOS96 53a6 overstabilizes the β-sheet secondary structure in hIAPP and that the AMBER model can more accurately reproduce the true structure of hIAPP in solution;36 this prompts reexamination of the true mechanistic details behind the hIAPP dimerization process.
While Fig. 1(a) highlights individual conformational clusters and their relative stabilities, little mechanistic information can be drawn from the free energy surface. In order to extract information regarding the role of a possible transient β-sheet intermediate formed in residues 20-29, we overlay the average amount of parallel β-sheet in residues 20-29 onto the free energy contours from Fig. 1(a). The overlay pinpoints two clusters of parallel β-sheet character in residues 20-29: one tight cluster at (0.33, 0.58) approximately 5.0 kJ/mol above the global minimum and a second diffuse cluster at (0.73, 0.73) approximately 10 kJ/mol above the global minimum. If a β-sheet intermediate is indeed formed in residues 20-29 during the dimerization process, the transition pathway must pass through these regions of the elevated parallel β-sheet mapped in Fig. 1(b).
However, there are two such clusters of high β-sheet character in residues 20-29 and little information to distinguish whether these two separate regions in CV space share degenerate hIAPP conformations. It is difficult to pinpoint a specific aggregation pathway on this free energy landscape, even more so considering that the two Q parameters are only dependent on a subset of all the residues in the system. In order to obtain details on the mechanism of aggregation, we perform finite temperature string method calculations using collective variables that may better capture conformational changes across the entire polypeptide. While performing yet another round of metadynamics with these more general collective variables might be prohibitively expensive, string methods inherently focus on the transition region, making them an attractive alternative for studying the dimerization process.
We use the finite temperature string method to examine the hIAPP dimerization process along two collective variables: the amount of parallel β-sheet formed between hIAPP chains and the distance between the centers of mass of each hIAPP monomer.
B. hIAPP dimerization pathway
Figure 2 shows the final transition pathway for the hIAPP dimerization process, obtained from finite temperature string method calculations (see Sec. IV). One end of the string corresponds to the disordered state of the hIAPP dimer (labeled Structure I, located at dCOM = 1.31 nm and ) and the opposite end of the string corresponds to the fully formed dimer (Structure V, at dCOM = 0.40 nm and ). Representative snapshots of the system are shown alongside the pathway in Fig. 2, highlighting key structural changes that occur during the dimerization process.
In the disordered configuration shown in Panel I, no secondary structure is formed between the two amylin chains. However, any individual disordered chain may independently form α-helix or β-strand motifs. Panel II highlights an intermediate structure, consistent with previous studies, exhibiting a parallel β-sheet structure that has been suspected to play a role in amyloid aggregation.28 Increasing amounts of parallel β-sheets are formed between the two amylin monomers as dimerization progresses, as witnessed in Panels III and IV, which show parallel β-sheet forming in the C-termini before advancing to the N-termini to form the full fibrillar dimer in Panel V.
The free energy profile along this dimerization pathway is shown in Fig. 3. An initial free energy barrier of approximately 7 kBT is discovered, and the images provided in Fig. 3 indicate that this barrier corresponds to the formation of the intermediate β-sheet structure. As the system climbs the 7 kBT energy barrier, the two amylin monomers approach each other and align; as the system traverses past the peak of the barrier, the transient intermediate β-sheet forms, accompanied by a drop in free energy.
Figure 4 shows the average secondary structure (calculated via the DSSP algorithm40) per residue for the dimerization process, with barriers corresponding to those in Fig. 3 marked. The intermediate β-sheet structure formed in the Barrier I region is found to be localized not only in the previously proposed region in residues 20-29 but also in the L12A13 region recently proposed by Maj and co-workers.33 Note that this is slightly shifted toward the N-termini compared with the GROMOS model.
A small free energy barrier is observed at the opposite end of the pathway, on the order of 2 kBT. By examining the corresponding snapshots, it becomes clear that this barrier is associated with a conformational rearrangement of the two amylin monomers, transitioning from two extended chains stacked side-by-side in parallel to a more compact “ribbon-like” structure exhibiting a slight twist. This rearrangement is coupled with additional formation of the N-termini β-sheet structure, as reflected in the plot of secondary structures in Fig. 4. Figure 4 also reveals that the additional formation of β-sheet in the N-termini, and its associated free energy stabilization, comes at the cost of relinquishing a fraction of β-sheet in the C-termini. Furthermore, formation of the fully dimerized structure corresponds to a shift in localization of β-sheet character away from the center of the hIAPP chain toward the termini.
Additional mechanistic details can be uncovered by tracking the number of protein-protein and protein-water contacts (within 3 Å) over the course of the dimerization process (Fig. 5), as well as the entropic and enthalpic contributions to the changes in free energy (Fig. 6). The latter is accomplished by calculating a change in potential energy ΔU and a change in free energy ΔA for each sampled Voronoi cell along the reaction coordinate with respect to the disordered bin, followed by calculating the entropic changes using ΔA = ΔU − ΔTS. Figure 5 shows the number of protein-protein contacts steadily increasing over the course of dimerization while protein-water contacts fluctuate around a value slightly lower than found in the disordered state. The corresponding enthalpic and entropic profiles are shown in Fig. 6. Taken together, Figs. 5 and 6 suggest that enthalpic and entropic changes fluctuate as more protein-protein contacts are formed, before a sharp jump in entropic contributions occurs, corresponding to Panel E in Fig. 3. By overlaying multiple root-mean-square deviation (RMSD) aligned snapshots taken from the same trajectory as Panel E (Fig. 7), we see the jump in entropy reflected in the conformational degeneracy of the dimer at both termini, while the core of the dimer remains compact, the ends of each chain explore many configurations. The entropic contribution then drops as we move to Panel F in Fig. 3, stabilizing the final dimer structure by approximately 2 kBT.
III. CONCLUSIONS
Bias-exchange metadynamics simulations have been used to study the thermodynamics of hIAPP dimerization using the AMBERff99SB*-ILDN force field. A global free energy minimum corresponding to the disordered dimer has been identified, as well as a metastable free energy minimum corresponding to the fully formed dimer has been identified, whose value is 7.4 kJ/mol or approximately 3 kBT higher.
The dimerization pathway for hIAPP, determined from finite temperature string method calculations, reveals an energy barrier of approximately 7 kBT associated with the formation of an intermediate β-sheet structure. Interestingly, this structure is localized in the previously proposed region in residues 20-29, as well as in the more recently proposed L12A13 region. Consistent with the results of metadynamics simulations, the fully formed dimer corresponds to a local free energy minimum whose energy is approximately 4.5 kBT higher than that of the disordered dimer. The fully formed dimer lies in a shallow well of approximately 2 kBT; importantly, further free energy decomposition suggests that the final dimer structure is entropically stabilized.
While a moderate (7 kBT) free energy barrier is associated with the formation of the intermediate dimer structure, the question that now arises is whether aggregation will remain uphill in free energy as the oligomer grows larger, or whether the aggregation process will become favorable after a certain size oligomer has been formed, thereby initiating exponential fibril growth. Other remaining questions include whether higher order aggregates nucleate from the globally stable disordered state, an intermediate structure formed during dimerization, or the locally stabilized structured dimer, and how these higher order aggregation events proceed. Such simulations are considerably larger in magnitude and are being pursued in our laboratory—the results will be presented in a future publication. Further investigations of that nature may shed light on the formation of higher order aggregates and the effects of different physiologically relevant environments on the aggregation process.
IV. METHODS
A. Human amylin dimer
We base the design of our dimer system on that of Chiu and de Pablo,34 with the major difference being the change in force field. The simulated system consists of two hIAPP molecules, 26 626 water molecules, and six chloride ions. The amino acid sequence for hIAPP is KCNTATCATQRLANFLVHSSNNFGAILSSTNVGSNTY. The C-termini of each polypeptide are amidated, and the side chains of Cys2 and Cys7 are linked by a disulfide bond. Protonation states of all ionizable functional groups are assigned on the basis of their pKa values in aqueous solution at a pH of 7.0. Each hIAPP molecule carries a formal charge of +3, and chloride counterions are included to ensure zero net charge in the system. Polypeptides and ions were modeled by the AMBERff99SB*-ILDN force field,39,43 and water was modeled by the TIP3P model.44 The system was placed in a periodic cubic box of side length 9.4 nm. Coulombic forces were calculated using the particle mesh Ewald algorithm,45,46 temperature coupling at 298 K was achieved using the Nosé-Hoover thermostat,47 volume was held constant, and simulations were performed with a time step of 2 fs. The LINCS (Linear Constraint Solver) algorithm was used to constrain hydrogen bond lengths to equilibrium values.48
B. Bias-exchange metadynamics
Bias-exchange metadynamics (BE-MetaD) simulations were performed to construct a free energy landscape for the hIAPP dimer. BE-MetaD performs conventional metadynamics in parallel on multiple replicas of the system, with atomic coordinates periodically exchanged between replicas à la parallel tempering.49 As in standard metadynamics,50,51 the final bias potential consists of a sum of small, Gaussian potentials deposited periodically along the system trajectory in collective variable (CV) space. Free energies are calculated from the combined statistics using the weighted histogram analysis method (WHAM).52 In this study, and in the following equations, each replica is biased along one CV at most. The cumulative bias potential of replica i at time t, denoted as , is50
where ξi denotes the CV, ξ(t) is the atomic coordinate vector, and the times at which Gaussian potentials were previously deposited are t0. The amplitude and width of the Gaussian potentials are W and σi, respectively. Exchanges of the atomic coordinates and velocities are attempted every 50 ps between randomly selected pairs of replicas. Coordinates of replicas i and j are exchanged with probability pij,
where kB is the Boltzmann constant and T is the temperature. Since all replicas are held at the same, constant temperature, the exchange probability is independent of the energy from the true, non-biasing MD interactions. Beyond some elapsed simulation time tfill, the system diffuses freely in CV space and a time average of the inverse cumulative bias potential is an estimator of the underlying free energy surface. Using WHAM, statistics from all replicas may be combined and the free energy A may be calculated as a function of any CV ξ′41,52
where k is the number of replicas, is the average biasing potential acting along ξ′ in replica j, and the fj is the normalization constant for replica j, calculated iteratively through the WHAM algorithm. represents the number of times that replica i visits CV value ξ′.
Two types of CVs were chosen for the BE-MetaD calculations: Qres u-v, a measure of similarity of to the proposed amylin fibril structure derived from ssNMR experiments and , a measure of parallel β-sheet content. In each case, a range of amino acid sequences is bounded by residues u and . Specifically, Qres u-v is the RMSD structural similarity of residues u through to the same set of residues from the ssNMR structure,34,53
where rij and are the distances, in Å, between backbone atoms i and j in the sampled and reference configuration, respectively. The angle brackets indicate an average over all pairs of backbone atoms belonging to residues in the sequence bounded by residues u and . The value of Qu-v ranges from 0, which indicates no similarity, to 1, which indicates an identical conformation. In this study, we adopt Qres 8-16 and Qres 27-35 as CVs. The reference structures are residues 8-16 and residues 27-35, respectively, of the ssNMR structure,1 i.e., the N- and C-terminal β-strands in the fibril.
The parallel β-sheet content of a particular sequence bounded by residues u and is defined as54
where the sum runs over all pairs of three-residue segments bounded by residues u and in both molecules. RMSD measures the root mean square deviation, in Å, of the positions of the N, Cα, Cβ, C, and O backbone atoms in those (3 + 3)-residue blocks from those in an ideal, parallel β-sheet. Put another way, counts the number of pairs of three-residue segments similar to the ideal parallel β-sheet, scaled by a switching function.
We adopt , , and as CVs. Biasing along and accelerates sampling of all parallel β-sheets involving those segments, not just the ones in the ssNMR structure. Biasing along accelerates sampling of parallel β-sheets in the central regions of the molecules, which have been identified as a key intermediate in the formation of fibrils.28
BE-MetaD simulations were carried out for a system consisting of six replicas. Five replicas were subject to the metadynamics biasing potential acting along one of the five CVs defined above: Qres 8-16, Qres 27-35, , , or . The sixth replica evolved with zero biasing potential but was allowed to exchange coordinates with the other replicas according to Eq. (2). For all replicas, W = 2.0 kJ mol−1 and Gaussian potentials were deposited every 5 ps. Widths of the deposited Gaussian potentials were σi = 0.02 for Qres 8-16 and Qres 27-35; widths of the deposited Gaussians were σi = 0.2 for , , and . The total simulation time was 500 ns. The filling time tfill was 400 ns. Simulations were conducted using the GROMACS 4.5.5 simulation package55,56 and a modified version of the PLUMED 1.3 plugin.57 WHAM calculations were conducted using the METAGUI plugin for VMD (Visual Molecular Dynamics).41,42
C. Finite temperature string method
In order to identify a pathway for hIAPP dimerization, we utilize the finite-temperature string method,58 which calculates a transition pathway as a series of local points (or “nodes”) connected by a smooth curve (or “string”) in collective variable space. Here, we investigate the free energy landscape described in terms of two intuitive collective variables that depend on all residues in the system: (1) parallel beta sheet character , following Eq. (5) and using all residues in the system and (2) a measure of the spatial separation between individual hIAPP chains. We use the distance between the centers of mass for each hIAPP monomer.
The string is discretized into 16 nodes, located in collective variable space at zα, where α indicates the index along the string (α = 0, …, 15). The nodes split collective variable space into a Voronoi tessellation, where each node has an associated Voronoi cell consisting of the points in CV space closer to itself than any other node. We assume Euclidean geometry. At every string method iteration, each Voronoi cell is sampled such that no bias is applied while the system is within the boundaries of its own Voronoi cell and a harmonic restraining potential is applied when the system departs from its own cell
Each string method iteration samples every Voronoi cell for 100 ps. We track the running average of each node’s location in collective variable space since the first iteration . At the nth iteration, the string is updated according to
where Δτ is chosen to be 0.1 and the smoothing parameter rα is equal to 0 for nodes on either end of the string (α = 0 or 15), otherwise
where smoothing parameter κ is 0.1 and the number of nodes along the string N is 16. After each update to the string, a cubic spline interpolation is drawn through the 16 nodes and the nodes are redistributed along the string to maintain equal arc-length between adjacent nodes.
We iterate through these steps until the string converges, after which we run a secondary string method calculation with only N = 4 nodes. This additional string method calculation aims to ensure discovery of the true fibrillar state, which BE-MetaD suggests may lie in a relatively narrow, isolated free energy basin. This “miniature” string is initialized from an idealized dimer structure (extracted from the ssNMR mature fibril structure1) to the fibrillar end of the original, converged 16-node string, which is held pinned in CV space throughout the evolution of the 4-node string.
Upon convergence of the secondary string, the free energy is computed along the final composite string, consisting of the 16-node string and 4-node string stitched together. This is done by calculating πα, the equilibrium probability of the system to be found in Voronoi cell α, which is then used to calculate the corresponding free energy Aα,59
To improve resolution of the resulting free energy profile, we further discretize the original 16 node string to 32 nodes, resulting in the complete pathway being described by a total of N = 35 Voronoi cells. Each of the Voronoi cells is sampled using the same soft wall restraints described in Eq. (6), for 50 ns for multiple runs. For each system sampling cell α, we collect Tα, the total simulation time spent within cell α, as well as Nαγ, the number of times the system escapes into a neighboring cell γ. The equilibrium probabilities πα are calculated with the following system of equations, where is the rate of escape from cell α into γ:
ACKNOWLEDGMENTS
We thank Dr. Kyle Hoffmann and Lucas Antony for helpful discussions. The results presented in this work were obtained using the computational resources of the University of Chicago Research Computing Center and the Laboratory Computing Resource Center at Argonne National Laboratory. This work is supported by the National Science Foundation, Grant No. DMR-1710357.