Insight into Cross-Amyloid Interactions and Morphologies: Molecular Dynamics Simulations of Model Peptide Fragments of Amyloid-β (Aβ16-22) and Islet Amyloid Polypeptide (IAPP20-29)

Amyloid-beta (Aβ) and islet amyloid polypeptide (IAPP) are small peptides, classified as amyloids, that have the potential to self-assemble and form cytotoxic species, such as small soluble oligomers and large insoluble fibrils. The formation of Aβ aggregates facilitates the progression of Alzheimer’s disease (AD), while IAPP aggregates induce pancreatic β-cell apoptosis, leading to exacerbation of Type 2 diabetes (T2D). Cross-amyloid interactions between Aβ and IAPP have been described both in vivo and in vitro, implying the role of Aβ or IAPP as modulators of cytotoxic self-aggregation of each peptide, and suggesting that Aβ-IAPP interactions are a potential molecular link between AD and T2D. Using molecular dynamics simulations, “hot spot” regions of the two peptides were studied to understand the formation of hexamers in a heterogenous and homogenous peptide-containing environment. Systems of only Aβ(16-22) peptides formed antiparallel, β-barrel-like structures, while systems of only IAPP(20-29) peptides formed stacked, parallel beta strands and had relatively unstable aggregation structures after 2 μs of simulation time. Systems containing both Aβ and IAPP (1:1 ratio) hexamers showed antiparallel, β-barrel-like structures, with an interdigitated arrangement of Aβ(16-22) and IAPP(20-29). These β-barrel structures have features of cytotoxic amyloid species identified in previous literature. Ultimately, this work seeks to provide atomistic insight into both the mechanism behind cross-amyloid interactions and structural morphologies of these toxic amyloid species. Statement of Significance Molecular knowledge, biophysical characterization, structural morphologies, and formation pathways of amyloid oligomers - specifically low-molecular weight, cross-amyloid oligomers - remain preliminary and undefined. Characterizing interactions between homogenous and heterogenous amyloid oligomers is of great interest given that certain oligomer morphologies contribute to cytotoxicity, eventually resulting in comorbid diseases such as Alzheimer’s disease (AD) and Type 2 Diabetes Mellitus (T2DM). Utilizing model systems (e.g., fragments of full-length peptides) and molecular dynamics (MD) simulations to probe the biophysical underpinnings of cross-amyloid oligomer structures is the first step in understanding the dynamics, stability, and potential modes of cytotoxicity of these species, providing important insights into targetable biomolecular structures.


Introduction
The accumulation of amyloidogenic peptides into morphologically variable aggregates and highly-ordered plaques are hallmarks of amyloid diseases [1][2][3] . Amyloid diseases constitute a variety of degenerative disorders, including Alzheimer's Disease (AD) and Type 2 Diabetes (T2D). In 2020, T2D and AD affected an estimated 6 million 4 and 30 million 5 Americans, and are the sixth and seventh leading causes of death in the US, respectively 6 . AD, the most prevalent form of dementia, is characterized by neuronal cell degradation and death 7 . T2D is characterized by insulin insensitivity and impaired pancreatic β-cell function, and is exacerbated by poor diet, lack of exercise, and a host of other genetic and environmental factors 8,9 . While both AD and T2D are characterized as individual pathologies, these diseases are highly comorbid. Evidence implicates AD and T2D in mutual pathogenic pathways via overlapping molecular mechanisms [10][11][12] , and each increases the risk for the other. [13][14][15] . Furthermore, AD pathology alters insulin signaling pathways in the brain 16 , suggesting AD is implicated in Type III Diabetes 17 .
AD and T2D pathologies are associated with the amyloidogenic peptides amyloidβ (Aβ) [18][19][20][21] and the islet amyloid polypeptide (IAPP or amylin) [22][23][24][25] , respectively. Both Aβ and IAPP are capable of self-assembling [26][27][28][29] into cytotoxic species of variable morphologies, specifically fibrils [30][31][32][33] and oligomers 20,34,35 . These oligomers, particularly low-molecular weight oligomers, have been identified as key contributors to the onset of each disease and are capable of cross-amyloid interactions. [36][37][38] Elucidating the mechanisms by which these amyloids aggregate and their resultant structural morphologies can provide insight into the pathogenic events of each disease, allow assessment of the extent of their comorbidity at the molecular level, and provide insight into therapeutic design to disrupt cross-amyloid formation. Fibrillar amyloid morphologies generally exhibit a cross-β structure, composed of repeating β-strand monomers arranged in a steric zipper 39,40 . Oligomers formed onpathway to fibrillar aggregate formation are thought to be the primary cytotoxic species in amyloid diseases 41,42 . These have been found to non-specifically destabilize cellular membranes and cause ionic leakage and impaired Ca 2+ regulation by the formation of βbarrel structures [43][44][45] . Amyloid β-barrels are thought to be composed of a low number of peptide subunits, each assuming a U-shaped 'β-strand-turn-β-strand' motif common to most amyloids, and stacked annularly, parallel, and usually in-register to form the barrel. 46,47 Cytotoxic β-barrels comprised of out-of-register subunits have also been identified 48 .
Common structural features of amyloid oligomers and their mechanisms of membrane perturbation have yet to be described at an atomistic level, likely due to their transient nature and widely varying structural morphologies. Molecular dynamics (MD) simulations are ideal for investigating oligomerization events, allowing us to characterize their formation, structure, and stability. To date, MD simulations have been utilized to describe oligomerization of full-length Aβ 49, 50 and IAPP 51,52 , and to determine which portions of their sequences are necessary and sufficient for inducing aggregation and providing stability. In particular, the core hydrophobic regions, composed of residues 16-22 of Aβ (KLVFFAE) and residues 20-29 of IAPP (SNNFGAILSS), have been shown to exhibit similar aggregation potential as full-length peptides 53,54 , and thus are of primary interest as model peptide fragments for studying aggregation.
Previous work characterizing fragments of Aβ and IAPP has provided foundational evidence for modeling these particular fragments. Residues 20-29 of IAPP (IAPP (20-29)) form the spine of IAPP fibrils under physiological conditions 40 , indicating the ability of this fragment to participate in β-strand formation and its importance for aggregate stability.
Furthermore, the primary sequence of this fragment is poorly conserved in IAPP analogues from other species that do not exhibit amyloidogenic properties, 55 indicating that this fragment is influential in the human IAPP aggregation pathway. Likewise, residues 16-22 of Aβ (Aβ (16)(17)(18)(19)(20)(21)(22)) include the central hydrophobic core of full-length Aβ, which is capable of forming aggregates itself 56 and is essential for the aggregation of the full-length peptide. 57 Furthermore, residues 16-22 of Aβ exist in a natively helical state; a conformational transition of this region to a coiled state capable of integrating into βstrands is considered a major rate-limiting step in Aβ aggregation 36 , and interactions between IAPP and this region of Aβ may facilitate this conformational change and subsequently promote co-aggregation. 36,58 Full-length Aβ and IAPP share approximately 25% sequence identity and 50% sequence similarity, with the greatest conservation occurring within these "hot spot" regions 59 ; as such, these fragments likely play a central role in cross-amyloid aggregation.
Here, we utilized MD simulations to observe Aβ (16)(17)(18)(19)(20)(21)(22) and IAPP(20-29) hexamer formation, using both homogeneous and mixed systems, to provide detailed insight into the structural transitions and architecture of these amyloid fragments. This work describes the underlying structural organization and inter-peptide interactions that these model peptides sample in route to an organized, β-strand rich oligomer structure.
The steepest-descent method was employed for energy minimization, after which all systems were subjected to an equilibration by a canonical (NVT) ensemble to maintain temperature of the system at 310K using the Berendsen weak coupling method 64 , for 100 ps. Each replicate system utilized different, randomly selected starting velocities assigned at the onset of NVT. An isothermal-isobaric (NPT) ensemble was then employed for a subsequent 100 ps to equilibrate the systems to a pressure of 1 bar using the V-rescale modified Berendsen temperature coupling method 64, 65 and Parrinello-Rahman barostat. 66,67 Positional restraints were imposed on peptide atoms during equilibration and released afterward. Bond lengths were constrained using Parallel Linear Constraint Solver (P-LINCS) 68 with an integration time step of 2 fs. The smooth Particle-Mesh Ewald method (PME) 69,70 was used for calculating long-range electrostatics using cubic interpolation. A Fourier grid spacing of 0.16 nm, three-dimensional periodic boundary conditions, and a real-space coulomb cutoff of 1.4 nm were also utilized. Input parameter files, coordinate files, and scripts utilized in this study are located on our Open Science Framework page (https://osf.io/82n73/). MD production was performed for 2 µs for all replicates in each molecular system for a total simulation time of 24 µs. Root-mean-square-deviation (RMSD) analysis was used to determine simulation convergence for each system; subsequent, aggregate analysis was then conducted over the last 500 ns of each simulation to ensure the characterization of stable hexamers and provide adequate sampling time (6 µs shape was assessed by calculating eccentricity (e) from the three moments of inertia (I1, I2, and I3) and semiaxes a, b, and c, as previously described. 72,75,76 Other analysis, including hydrogen bonding, radius of gyration, and residue-residue contacts were performed using GROMACS analysis tools and in-house scripts. All averages presented were calculated over four replicates for each simulation system and are presented with the corresponding standard deviation. A two-tailed t-test was used for statistical analysis, with statistical significance determined as p < 0.05.

Results and Discussion
Detailed descriptions of amyloid morphologies and stepwise aggregation events are necessary for understanding and preventing the progression of amyloid diseases like AD and T2D. However, the stochastic nature of these events makes atomistic characterization of aggregate structure particularly difficult. Thus, MD simulations are a valuable tool for providing foundational insight into these mechanisms and identifying potentially targetable structures or co-interaction capabilities. Here, we utilized unitedatom MD simulations to model the aggregation of Aβ (16)(17)(18)(19)(20)(21)(22) and IAPP (20)(21)(22)(23)(24)(25)(26)(27)(28)(29)  Adoption of β-strand content is indicative of amyloid aggregation pathways. 77,78 As such, monitoring the progression of secondary structure allows us to identify canonical on-pathway aggregation. To properly model secondary structure of Aβ (16)(17)(18)(19)(20)(21)(22) and IAPP (20)(21)(22)(23)(24)(25)(26)(27)(28)(29), force field selection is critical and must be taken into account in system construction and in analysis outcomes. The GROMOS96 53a6 forcefield has been shown to model βstrand formation of Aβ in agreement with experimental data, with other force fields tending to over-stabilize α-helices 71 . For the purposes of this work and based on our previous force field comparison study, the GROMOS96 53a6 was chosen for modelling Aβ and IAPP fragments. A recent review highlights that there are three major problems to consider in simulating amyloid aggregation, including insufficiently accurate force fields, protein concentration, and time-scale limits. 79 More work is necessary to further refine force fields for both the field of IDPs, for model systems, and for aggregation processes, with this work contributing how one force field models similarities and differences between model systems of homogenous and mixed amyloid formation. This work focuses on homogenous and mixed amyloid structure formation on that facet and does not seek to be a comparision between force fields. Furthermore, recent work from discrete MD (DMD) simulations of full-length Aβ42 hexamer formation indicate a β-barrel formation in line with the results discussed below, as well as in vivo experimental validation, that highlighted the formation the β-barrel oligomers of Aβ42 80 . To describe aggregate structure with respect to secondary structure content, secondary structure over time for each system was assigned using the DSSP algorithm 81 . Aβ control hexamers showed a rapid timewise transition from 100% random coil at MD onset to 51% ± 4 β-strand structure over the last 500 ns of simulation (Figure 2A, Table 1). Likewise, IAPP control hexamers showed a similarly rapid progression from predominantly random coil at the onset of MD, stabilizing at 48% ± 8 β-strand structures in the last 500 ns of simulation ( Figure 2B, Table 1). These data match previous analyses of secondary structure content of amyloid oligomers, which adopt β-strand structure on-pathway to aggregation 82,83 . Fragments adopting experimentally determined secondary structure content supports the use of GROMOS 53a6 for modelling these fragments.
IAPP (20)(21)(22)(23)(24)(25)(26)(27)(28)(29) hexamers were the most polymorphic of any of the systems, exhibiting the most numerous cluster structures sampled over the last 500 ns of simulation (Table   2). Furthermore, the free energy of IAPP systems was on average the most favorable with the lowest free energy calculation of the final formed morphologies, with the lowest amount of hydrophobic SASA over the last 500 ns of simulation (Table 2). IAPP, like Aβ, benefits from strong interactions between hydrophobic residues (Figure 3). We predict that IAPP systems were more energetically favorable given the decrease in order and barrel formation observed as compared to Aβ only systems. Interestingly, IAPP is modestly more compact than the Aβ and mixed systems as assessed by radius of gyration, suggesting that hydrophobic packing is strongest in IAPP systems ( Table 2).
The reduced stability in IAPP systems is likely due to weaker interactions from flanking polar residues relative to Aβ; average minimum distance between the polar termini was longer than distance between charged terminal residues in Aβ ( Table 2). The polar termini are thus more flexible, and likely contribute to the structural polymorphism.
Combating amyloid diseases like AD and T2D also requires translating these differences in oligomeric structure and stability to relative toxicity. Out-of-register β-strand alignment, leaving staggered hydrogen bonds, has been implicated in the formation of toxic, cylindrin-like amyloids 91 . Out-of-register β-strands were observed in mixed Aβ (16)(17)(18)(19)(20)(21)(22) + IAPP (20)(21)(22)(23)(24)(25)(26)(27)(28)(29) replicates. We calculated the average number of hydrogen bonds between residues and water for each system to determine if out-of-register hydrogen bonds are of higher frequency in heterogeneous systems. Hydrogen bonding with water was increased in mixed systems, particularly with respect to N-terminal IAPP (20)(21)(22)(23)(24)(25)(26)(27)(28)(29) residues Asn21 and Asn22 ( Figure S23). If hydrophobic core contacts are preserved in full-length crossaggregates, out-of-register hydrogen bonds could potentially be present on the N-terminal end of IAPP. This supports the formation of β-barrels for full length cross-aggregates and may provide insight into their toxicity. Additionally, Aβ Glu22 residues were more likely to have an increase in number of hydrogen bonds with water, in both the mixed and homogenous system.
The formation of several distinct pore-like amyloid oligomers have been reported in the literature. It is hypothesized that these structures may confer toxicity by permeabilizing the membrane to ions, leading to cellular dysfunction. While the β-barrel structures observed appear pore-like in shape, the internal cavity of the barrel is generally packed with hydrophobic side chains. As such, barrel-like structures are accordingly unlikely to accommodate passage of significant quantities of water, corroborated by water shell analysis (Figure 4, Figures S24-S26). Similar hydrophobic packing has been observed in other studies examining the formation of pore-like Aβ structures 90 . In general, more apparent barrel-like structures appear to accommodate notably less water in the internal cavity due to increased hydrophobic packing, and therefore have less internal surface area for water contact. However, it is worth mentioning that permeability to solvent and ions may be dependent on oligomer size; Aβ 20-mers appear to have characteristics of ion channels 92 . Furthermore, disordered oligomers with increased hydrophobic SASA have been identified as cytotoxic 86,87 , suggesting that water permeability may be a factor in cytotoxicity.