Aggregation of Aβ40/42 chains in the presence of cyclic neuropeptides investigated by molecular dynamics simulations

Alzheimer’s disease is associated with the formation of toxic aggregates of amyloid beta (Aβ) peptides. Despite tremendous efforts, our understanding of the molecular mechanisms of aggregation, as well as cofactors that might influence it, remains incomplete. The small cyclic neuropeptide somatostatin-14 (SST14) was recently found to be the most selectively enriched protein in human frontal lobe extracts that binds Aβ42 aggregates. Furthermore, SST14’s presence was also found to promote the formation of toxic Aβ42 oligomers in vitro. In order to elucidate how SST14 influences the onset of Aβ oligomerization, we performed all-atom molecular dynamics simulations of model mixtures of Aβ42 or Aβ40 peptides with SST14 molecules and analyzed the structure and dynamics of early-stage aggregates. For comparison we also analyzed the aggregation of Aβ42 in the presence of arginine vasopressin (AVP), a different cyclic neuropeptide. We observed the formation of self-assembled aggregates containing the Aβ chains and small cyclic peptides in all mixtures of Aβ42–SST14, Aβ42–AVP, and Aβ40–SST14. The Aβ42–SST14 mixtures were found to develop compact, dynamically stable, but small aggregates with the highest exposure of hydrophobic residues to the solvent. Differences in the morphology and dynamics of aggregates that comprise SST14 or AVP appear to reflect distinct (1) regions of the Aβ chains they interact with; (2) propensities to engage in hydrogen bonds with Aβ peptides; and (3) solvent exposures of hydrophilic and hydrophobic groups. The presence of SST14 was found to impede aggregation in the Aβ42–SST14 system despite a high hydrophobicity, producing a stronger “sticky surface” effect in the aggregates at the onset of Aβ42–SST14 oligomerization.


Introduction
Alzheimer's disease (AD) is one of the most devastating neurodegenerative disorders of our time due to its high prevalence in the aging population, the challenges of early diagnosis and the lack of efficient therapeutics. AD is associated with the misfolding and formation of toxic aggregates of the amyloid β (Aβ) peptide [1,2]. It is believed that misfolding in spontaneously formed aggregates of Aβ chains eventually leads to a cascade of self-replicating misfolding events producing β-sheet rich amyloid deposits in the brain. However, mounting evidence suggests that relatively small prefibrillary oligomers rather than mature amyloid fibrils may be the primary toxic assemblies underlying the pathogenesis of the disease [2,3]. Insights into the molecular mechanisms of misfolding and aggregation, as well as co-factors that might influence these processes, remain incomplete.
Contributing to this status quo is a high level of heterogeneity in regards to both the building blocks and the architecture of early oligomeric assemblies. First, the Aβ peptide itself exists in multiple alloforms. Aβ 42 comprises two C-terminal residues, I41 and A42, not present in Aβ 40 . Although the Aβ 40 variant is more abundant than Aβ 42 , the latter is believed to be more pathogenic, primarily on the grounds of a relative increase in Aβ 42 in the brain of AD patients, and higher rates of fibrillization in-vitro [4][5][6]. Second, Aβ oligomers are an ensemble of highly dynamic assemblies [7,8]. Experiments indicate that the initial small aggregates tend to adopt a largely unstructured morphology [9][10][11]. Quaternary structures without pronounced alignments of peptide chains were also predicted by molecular dynamics (MD) simulations for dimers of Aβ 40 and Aβ 42 [12][13][14] and larger multimeric Aβ aggregates [15,16]. Characterizations by a variety of biophysical methods in vitro suggest that initially unstructured oligomers may either undergo a transition into more ordered β-sheet-rich conformations, or mediate formation of β-sheet-rich assemblies [9][10][11][12]. However, the majority of early non-fibrillary aggregates dissociate into monomers rather than convert into fibrils [7,8], albeit they are potentially toxic [8,11]. Neither the detailed molecular mechanisms of aggregation, nor the accompanying misfolding or the specific structures and morphologies of the various nonfibrillary and pre-fibrillary assemblies are sufficiently understood. On theoretical grounds, it has been inferred that the aggregation process is driven by a competition of hydrophobic collapse and hydrogen bonding, with the former favoring unstructured conformations, and the latter giving rise to β-sheets [17][18][19]. The toxicity of non-fibrillary and pre-fibrillary oligomeric assemblies has been hypothetically linked to the exposure of hydrophobic groups and unpaired β-strands at their surfaces [9,20], often referred to as a "sticky surface" effect. Recent modeling studies [21,22] indicate that accumulation of subtle structural perturbations of early aggregates at the onset of the oligomerization process may result in distinct aggregation pathways at later, more advanced stages of Aβ oligomerization. Altogether, both experimental and theoretical evidence suggests that molecular events occurring early in the process of aggregation play a key role in determining both the structure and toxicity of Aβ oligomers [14].
The potential role of extrinsic factors in the aggregation and amyloidogenic conversion is a relatively under-explored aspect of immense importance for both the basic understanding of the aggregation process and the rational design of therapeutic strategies [23,24]. In particular, toxic Aβ � 56 complexes have been hypothesized to require unknown cofactors for their assembly [25]. The cyclic neuropeptide somatostatin-14 (SST 14 ) was reported to promote the proteolytic degradation of Aβ through neprilysin induction [26]. The level of SST 14 in the brain decreases with ageing [26] and an accelerated decline is found in AD patients [27], potentially leading to an increase in steady-state Aβ levels. Recent experiments indicate that the cyclic neuropeptide somatostatin-14 (SST 14 ) is the most selectively enriched peptide in human frontal lobe extracts that binds oligomeric Aβ 42 aggregates [28]. Moreover, SST 14 's presence was found to inhibit fibrillization of Aβ 42 in vitro, while promoting formation of smaller oligomers, reminiscent of toxic Aβ � 56 complexes [28,29]. Interestingly, Aβ 42 but not Aβ 40 peptides were prone to delayed fibrillization under the influence of SST 14 . This effect was also not observed upon replacement of SST 14 with other cyclic neuropeptides, such as arginine vasopressin (AVP) [29]. Another recent biochemical study [30] investigated binding of SST 14 to a specific membrane-associated Aβ 42 tetramer, termed βPFO Aβ  . This report validated the ability of SST 14 to selectively interact only with oligomeric assemblies of Aβ 42 . Consistent with previous data gathered with soluble Aβ 42 oligomers [29], the authors reported that the binding interface between SST 14 with βPFO Aβ  may involve a central tryptophan (W8) within SST 14 . Whereas prior experiments [29] with soluble Aβ 42 oligomers had pointed toward a contribution of the N-terminal half to the SST 14 binding, the interactions with the membrane associated βPFO Aβ(1-42) oligomer seems to primarily rely on residues 18-20 within the C-terminal half of Aβ 42 . Although at first glance these results may seem contradictory, the N-terminal half of Aβ 42 is required for the formation of βPFO Aβ  , an observation that would have precluded the ability to detect binding of SST 14 in assays that relied on truncated Aβ  in prior binding studies. Both studies agreed that binding occurs predominantly to Aβ 42 oligomers and is not observed with oligomers formed from Aβ 40 . Because analyses of [30] were restricted to a highly purified membrane-associated βPFO Aβ  , it remains unclear whether an alternative binding pose or binding stoichiometry would be available in soluble Aβ 42 oligomers. Consistent with such a scenario, the binding constant for binding of SST 14 to βPFO Aβ  was approximately threefold higher than the previously observed Kd for binding of SST 14 to soluble Aβ 42 oligomers.
In order to elucidate how SST 14 might influence the onset of Aβ oligomerization, we have performed all-atom molecular dynamics (MD) simulations of model systems containing mixtures of Aβ 42 or Aβ 40 monomeric peptides and SST 14 molecules in explicit water. We also performed similar simulations for mixtures of Aβ 42 peptides and AVP molecules as negative controls. We investigated early stages of aggregation in each system, and analyzed the structure and dynamics of the early self-assembled aggregates. For this purpose, we combined wellestablished structural analysis tools with a novel essential collective dynamics (ECD) method developed in our group [14,16], which allows to accurately identify persistent correlations of motion (dynamics) in a molecule or supramolecular system without exhaustive conformational sampling, based on an original fundamental concept [31,32]. The method allows analysing dynamics correlations between selected pairs of atoms; characterising main-chain flexibilities; and identifying domains of correlated motion within the same framework. The presence of SST 14 was found to impede the aggregation in the Aβ 42 -SST 14 mixtures despite having a high hydrophobicity. We attribute differences in the structures and dynamics of the Aβ 42 -SST 14 , Aβ 42 -AVP and Aβ 40 -SST 14 systems to distinct tendencies of SST 14 and AVP to (1) interact with specific regions of the Aβ chains; (2) develop hydrogen bonds with Aβ peptides; and (3) expose hydrophilic or hydrophobic groups to the solvent.

Aggregation of Aβ peptides in the presence of small cyclic peptides
Initially, eight Aβ 42 or Aβ 40 chains and eight SST 14 or AVP molecules were placed in the simulation box in random positions, as outlined in the Methods section. In each system, the Aβ peptide chains are labelled A, B, C,. . .H, and the small cyclic peptides are labelled I, J, K,. . .P. Three independent 500 ns long MD simulations were conducted for each of the Aβ 42 -SST 14 , Aβ 42 -AVP and Aβ 40 -SST 14 systems in explicit OPC water. To distinguish between the three trajectories, each of them were numbered by I, II, and III. Unless otherwise stated, the examples displayed below originate from trajectories Aβ 42 -SST 14 -I, Aβ 42 -AVP-II, and Aβ 40 -SST 14 -II. These trajectories were selected as representative examples of leading trends for each of the systems. Fig 1 shows the configurations for the three representative systems after the equilibration (see Methods), and after 500 ns production MD simulations. Additional details on the structures obtained in these trajectories can be found in S1 and S2 Figs in the Supporting Information. Table 1 summarizes the aggregation status in all nine systems after 500 ns MD simulations.
After equilibration of the Aβ 42 -SST 14 system, Aβ chains were still relatively sparsely positioned, as shown in Fig 1A. After the 500 ns long production MD simulations, all three trajectories for this Aβ 42 -SST 14 system developed self-assembled aggregates. Most aggregates contained two or three Aβ 42 chains and one or two SST 14 molecules, as illustrated in Figs 1B and S2A for system Aβ 42 -SST 14 -I. Somewhat larger aggregates were observed occasionally along with dimers and trimers. In system Aβ 42 -SST 14 -II four Aβ 42 chains and five SST 14 molecules formed an aggregate, and in system Aβ 42 -SST 14 -III five Aβ 42 chains and three SST 14 molecules formed an aggregate ( Table 1). The SST molecules tend to attach to the surface of the Aβ 42 oligomers such that their N-and C-terminal arms face the solvent.
When the eight SST 14 molecules were replaced by AVP molecules, the starting configurations for the Aβ 42 -AVP system after equilibration also exhibited sparse positioning of the chains ( Fig 1C). However, the final structures from the three trajectories comprise large aggregates composed of mixed Aβ 42 and AVP ( Table 1). The final configuration for system Aβ 42 -AVP-II is shown in Figs 1D and S2B. In this case, the largest aggregate is an octamer containing all eight Aβ 42 chains and eight AVP molecules. In each of the other two Aβ 42 -AVP trajectories Aβ 42 peptides formed a hexamer and a dimer. The hexamers contained 6 or 7 AVP molecules and at least one AVP was bound to each of the Aβ 42 dimers.
For the Aβ 40 -SST 14 system containing eight Aβ 40 chains and eight SST 14 molecules, all three trajectories developed aggregates composed of seven Aβ peptides and six SST 14 chains. One Aβ 40 chain remained free in each of the trajectories (Table 1). Fig 1E shows a representative configuration in system Aβ 40 -SST 14 -II after the equilibration, and Figs 1F and S2C depict this system after 500 ns of MD simulation. In this case, seven Aβ 40 and six SST 14 chains formed a curved semi-hollow structured oligomer. Six SST 14 molecules are inserted inside the aggregate, whereas two SST 14 are attached to the remaining free Aβ 40 chain preventing its attachment to the larger aggregate.
Due to the putative importance of the C-terminal residues in defining the aggregation behaviour of Aβ peptides [4,19], the locations of the Aβ 42 and Aβ 40 C-termini are marked, respectively, by blue and green spheres in Fig 1. From Fig 1B and 1D it appears that the C-termini of Aβ 42 peptides (blue spheres) tend to remain at the surface of the self-assembled aggregates in both the Aβ 42 -SST 14 -I and Aβ 42 -AVP-II systems. In the Aβ 40 -SST 14 -II system, the Ctermini of several Aβ 40 chains were buried inside the aggregate (Fig 1F). The N-termini of both Aβ 42 and Aβ 40 chains tend to remain at the surface of the oligomers in all systems.
The aggregation process is accompanied by a build-up of hydrogen bonds (HBs) in each system. then the build-up slows down although the number of HBs fluctuates considerably during subsequent simulations. In the system Aβ 42 -AVP-II a persistent increase in the total number of HBs continued throughout the entire simulation. By 500 ns this system developed more HBs than the other two systems, presumably due to the large size of the aggregate. As can be seen in S1 Table, the average number of HBs over three Aβ 42 -AVP trajectories at the end of the simulations was greater than similar averages in the other two systems as well.
In order to provide more detailed information on the hydrogen bonding, we analyzed the number of HBs between and across Aβ chains and the small cyclic peptides separately. In Fig  2, red lines represent the evolution of the number of HBs between Aβ chains in the three representative systems. As is evident from the figure, Aβ-Aβ bonds contribute a significant portion (60-80%) of the total number of HBs, and their time dependencies resemble those of the total number of HBs. According to S1 Table, at the beginning of the simulations the average numbers of Aβ-Aβ bonds were close in Aβ 42 -SST 14 and Aβ 42 -AVP systems (49 and 52, respectively), and slightly above than in the Aβ 40 -SST 14 system (44). In the course of the simulations, the average numbers of Aβ-Aβ bonds increased approximately twice in each of the systems.
The most dramatic differences were observed for hydrogen bonding across Aβ chains and small cyclic peptides (blue lines in Fig 2), and between small cyclic peptides (green lines in Fig  2). After equilibration and early into the production runs Aβ-Aβ, SST 14 -SST 14 , and AVP-AVP bonds were almost exclusively intra-molecular, due to distances between the chains in the initial structures. Examples illustrating this can be found in S1 Fig. During the first 20 ns after equilibrations, approximately 20 SST 14 -SST 14 bonds in average were formed for each of two somatostatin-containing systems (S1 Table). In striking contrast, only~9 AVP-AVP bonds have developed. In the course of the simulations the number of SST 14 -SST 14 bonds changed only slightly adopting average values of 23 (for Aβ 42 -SST 14 system) and 17 (for Aβ 40 -SST 14 system) over the last 20 ns. At the same time, the average number of AVP-AVP bonds decreased barely remaining above~7 over the last 20 ns. Such a large difference cannot be explained by different numbers of HB acceptor sites in SST 14 and AVP (which are 24 and 16, respectively) and rather indicates differences in bonding behaviours of SST 14 and AVP. As can be seen in S1 Table, during the first 20 ns of production simulations when the small cyclic peptides just started binding to Aβ chains, the average values of cross-species hydrogen bonds were in the

PLOS COMPUTATIONAL BIOLOGY
proportion of Aβ 42 -AVP > Aβ 40 -SST 14 > Aβ 42 -SST 14 . More specifically, the greatest average number of 12 cross-species bonds was found in Aβ 42 -AVP system, and the smallest average number of 5 cross-species bonds was observed in Aβ 42 -SST 14 system. During the simulations the number of cross-species bonds increased by 3-5 times in each simulation, and exhibited a significant variability across individual trajectories (S1 Table and S3A Fig). However, the average number of Aβ 42 -AVP bonds remained approximately 2.4 times greater than the number of Aβ 42 -SST 14 bonds. Similarly to the other two systems, the number of Aβ 40 -SST 14 bonds varied across trajectories. In system Aβ 40 -SST 14 -II, for instance, the number of cross-species bonds is close to that in system Aβ 42 -AVP-II over significant parts of the two trajectories (S3B Fig). However, the average number of Aβ 40 -SST 14 bonds is less than that of Aβ 42 -AVP (S1 Table  and S3A Fig). Overall, regardless of the size and morphology of these aggregates, SST and Aβ tend forming less HBs than AVP and Aβ. This explains the large total number of hydrogen bonds in AVP-containing systems.
In order to better understand how fluctuations in the number of hydrogen bonds relate to the morphologies, we examined the structures adopted by the three representative systems at several local maxima and minima of the cross-species HBs' dependencies on time.  14 chains, which appears to contribute to the observed fluctuations of the hydrogen bonding. For example, the structure at local minimum of 98 ns exhibits a monomeric ("free") Aβ 42 chain, three free SST 14 molecules, and one SST 14 molecule semi-detached from a small aggregate. Clearly there are more free chains than for 293 ns, where a maximum of HBs is observed. However, the local maximum of hydrogen bonding at 79 ns seems to exhibit a sparser morphology than the local minimum at 454 ns. We attribute this to the dynamic nature of the small aggregates and their propensity to restructure, which emerged as another factor contributing to the fluctuations. In the Aβ 42 -AVP system ( Fig 3B) small aggregates self-assembled early in the simulations. The local maximum of hydrogen bonding at 86 ns already exhibits four aggregates and only one free AVP molecule. This is followed by a local minimum at 157 ns, where we observe a large but sparse aggregate, two smaller ones, and one free AVP. A significant decrease in the number of HBs can be attributed to the ongoing change in the aggregation status. Two aggregates observed at 346 ns are significantly more compact, and they exhibit a pronounced increase in hydrogen bonding. By 442 ns the two aggregates have merged into one C-shaped oligomer; however, the accompanying restructuring resulted in a loss in hydrogen bonding. The system Aβ 40 -SST 14 -II ( Fig 3C) quickly developed two aggregates without any free chains resulting in a maximum in HB at 113 ns. Although the gross aggregation status was retained throughout the course of the simulation, the subsequent ups and downs of the hydrogen bonding suggest continuing restructuring of the aggregates.
After the 500 ns long MD simulations, the secondary structure of the self-assembled aggregates contains predominantly random coils and turns, as depicted in Fig 1. However, stable βsheets are also observed. S4 Fig shows the secondary structure evolution over the 500 ns of three representative trajectories, and Table 2 (Table 2). In most cases we observe anti-parallel, intramolecular β-sheets in a hairpin conformation of the chain, often involving Aβ C-terminal residues 37-40 (in Aβ 42 systems) or 36-39 (in Aβ 40 system) as part of the β-sheets. Regardless of the system, most inter-chain β-sheets involve Aβ residues from the central region between positions 23 and 35, although specific locations of βstrands may vary. In several simulations (Aβ 42 -SST 14 -II and III, and Aβ 40 -SST 14 -II and III) β-

PLOS COMPUTATIONAL BIOLOGY
strands from this region formed β-sheets with N-terminal residues of a different Aβ chain. The smaller SST 14 and AVP molecules appeared less prone to contribute to the stable β-sheet content (see

ECD correlations of motion in self-assembled aggregates
To probe collective motions of the Aβ 42 -SST 14 , Aβ 42 -AVP and Aβ 40 -SST 14 systems, essential collective dynamics analyses (ECD) [16,[31][32][33][34][35] have been performed. In the ECD method, principal eigenvectors of the covariance matrix calculated from MD trajectories are employed to characterize correlations of motion (dynamics) between atoms, as outlined in the Methods. Importantly, the method allows to reliably identify stable dynamical correlations without requiring exhaustive conformational sampling or convergence to equilibrium [31,32]. The dynamics correlations identified in the ECD framework can also be visualized in the form of domains of correlated motion [16,31,32], which represent relatively rigid parts of the system consisting of atoms moving coherently. Fig 5 shows the largest domains of correlated motion color-mapped onto the secondary structure of self-assembled aggregates in the three systems. The domains are colored blue, red, green, yellow, etc. in the order of decreasing size. In the Aβ 42 -SST 14 -I system shown in Fig 5A, the three largest domains of correlated motion colored blue, red, and green are located in central parts of the three aggregates. These domains include both Aβ and SST 14 chains that are strongly inter-correlated in Fig 4A, and contain most of the β-structure. Peripheral regions, which often include termini, tend to exhibit relatively independent dynamics. In the case of the Aβ 42 -AVP-II system shown in Fig 5B we observe three large domains as well; however, all three are included in a single dumbbellshaped aggregate. The domain colored blue contains Aβ chains A and G and AVP chains N and O, whereas the domain colored green consists mainly of Aβ chain F and AVP chain I. These two domains are located at the core of two sides of the "dumbbell", which are held together through the mediation of the third domain (colored red), which includes parts of Aβ chains C and E along with AVP molecule J. In the Aβ 40 -SST 14 -II system depicted in Fig

Main-chain flexibility profiles
Within the same ECD framework, main-chain flexibility profiles [16,[33][34][35]  In the Aβ 42 -SST 14 system, most Aβ chains show a decreased flexibility at the last 20 ns of production MD simulation compared to the first 20 ns (Fig 6A), as one could expect since oligomerization limits mobility of the chains. The flexibility profiles of individual Aβ chains tend to adopt oscillating shapes with several maxima and minima. Although by the end of the simulations N-and C-terminal ends of the chains remain flexible, deep minima developed within 10-12 residues from one or both of the termini in most chains. This is especially the case for chains B, D, and H, which have formed a trimeric self-assembled aggregate (   simulation. An analysis of the trajectory indicated that initially these molecules were surrounded by neighbouring Aβ 42 chains (S1A Fig), which somewhat limited their motion. By the end of the simulations the Aβ 42 chains and other SST 14 molecules formed three aggregates, whereas the molecules K, O, and P remained free. Eventually, remaining in solution resulted in a less constrained motion for these SST molecules by the end of the simulation. Fig 6B depicts the main-chain flexibility profile for the Aβ 42 -AVP system. After the equilibration most chains exhibit relatively uniform flexibility patterns, whereas by the end of the simulation pronounced differences emerged across the various chains. Both termini of chains B and C, as well as the N-terminus of chain F, and the C-terminus of chain G adopted higher flexibility levels in the course of the simulation, which can be explained by their location at the periphery of the self-assembled aggregate. In contrast, extensive regions of chains A, E, F, and G, which are located at the core of the three domains of correlated motion (Figs 5 and S6B), exhibited low flexibility. Flexibilities of AVP molecules tended to follow those of the Aβ chains with which they interacted over the last 20 ns of the simulation.
The main-chain flexibility profile of the Aβ 40 -SST 14 system is shown in Fig 6C. Extensive regions of Aβ chains A, B, D, G and H, which are located in the interior regions of the aggregate, developed a decrease in flexibility during MD simulation. In contrast, chain F and most of chain C adopted relatively high flexibilities due to their less constrained positions. Central regions of SST 14 molecules K, M, O, and P exhibited relatively low main-chain flexibility over the last 20 ns of the MD simulation. These chains were actively involved in the oligomerization process, which resulted in their insertion inside the self-assembled aggregate where they became a part of the largest domain of correlated motion (see also Figs 5C and S6C). SST 14 molecules I and M, which attached to the free Aβ chain F, exhibited high flexibility during last 20 ns of the simulations. Table 3 lists the main-chain flexibility of the Aβ 42 or Aβ 40 C-termini averaged over eight Aβ chains in each of the three MD trajectories in each system. The table also lists the results of averaging over the three MD trajectories for each system. Further details of C-terminal flexibilities in individual Aβ chains can be found in S4 Table in the Supporting Information. As  Table 3 illustrates, the average flexibility of the C-termini is greater in the AVP-containing systems than in the SST-containing systems. AVP-containing systems also exhibit the greatest variation of C-terminal flexibility across individual trajectories. Interestingly, the average flexibility of the C-termini in the Aβ 40 -SST 14 system is very close to that in the Aβ 42 -SST 14 system (Table 3), despite the absence of two hydrophobic residues at the C-terminus of Aβ 40

Intermolecular interactions in self-assembled aggregates
In order to elucidate the molecular mechanisms behind the influence of the small cyclic peptides onto the self-assembled aggregates, we identified the SST 14 and AVP molecules that exhibited the strongest inter-molecular dynamics correlations by averaging over the last 20 ns, in each MD trajectory and for all the Aβ 42 -SST 14 , Aβ 42 -AVP and Aβ 40 -SST 14 systems. Fig 7  depicts close-ups of the identified highest correlated regions involving small cyclic peptides and Aβ chains. Fig 7A and 7B illustrate, respectively, the interactions between an Aβ dimer (A+C) and the SST chains L and M attached to it, and between Aβ chains F+G and the SST chain N in system Aβ 42 -SST 14 -I. Fig 7C shows the interactions between Aβ chains G and B and SST chain P, and Fig 7D illustrates the interactions between Aβ chain F and SST chain J in system Aβ 42 -SST 14 -II. In these examples, the strongest inter-molecular correlations involve pairs of residues: K16-F7, G33-T10, and V40-W8, (Fig 7A), G33-T12 and M35-T10 (Fig 7B), I32-N5 and V40-K4 (Fig 7C), and K28-N5 (Fig 7D). Remarkably, the SST residues involved in strong correlations with Aβ 42 are often hydrophobic (F7, W8), or carry a long hydrocarbon side chain (K4), and/or have hydrophobic neighbors such as F6 or F11. Similarly, the pairing Aβ residues often carry hydrophobic groups (M35, V40) and/or have hydrophobic neighbors (L17, I31, I32, L34, I41). These hydrophobicity-driven interactions of side chains are accompanied by backbone-backbone hydrogen bonding, for example L16-F7, G33-T10 (Fig 7A), L34-T12 and V36-T10 (Fig 7B), or I41-K4 (Fig 7C). Complementary to hydrophobic interactions and hydrogen bonding, we also detected π-π interactions between residue F20 of Aβ 42 and F7 of SST 14 (Fig 7D).
Next, Fig 7E and 7F depict the strongest interactions between Aβ 42 and AVP molecules identified in systems Aβ 42 -AVP-I and Aβ 42 -AVP-II. In the first system illustrated by Fig 7E, AVP residue R8 from chain N developed strong dynamics correlations with D23 of the Aβ 42 chain E. Moreover, residue C1 developed contacts with residue F20 located nearby, although this interaction is somewhat weaker. In this case, both AVP residues have formed hydrogen bonds with their Aβ 42 counterparts resulting in formation of a cross-species β-sheet ( Fig 7E). Fig 7F illustrates the strongly correlated pair consisting of Aβ 42 residue K16 (chain G) and AVP residue Y2 along with a weaker interaction of closely positioned pair of E22 (chain A) and Q4. Despite an involvement of hydrophobic contacts between Y2 and a hydrocarbon group of K16, electrostatic forces rather than hydrophobicity appear to play a role in this interaction.
Interactions between Aβ 40 chains and SST 14 molecules are illustrated with examples from system Aβ 40 -SST 14 -I ( Fig 7G) and system Aβ 40 -SST 14 -II (Fig 7H). In system Aβ 40 -SST 14 -I, the strongest dynamics correlations involve residues M35, V36, and G37 from Aβ chain F and residues T10 and T11 from SST chain J. Backbones of residues M35 and G37 developed hydrogen bonds with F11 and T10, respectively, whereas residue V38 formed a hydrogen bond with T10 as well as hydrophobic interaction with F11. In system Aβ 40 -SST 14 -II, Aβ residues V40 (chain A) and K28 (chain B) exhibit strong dynamic correlations with residues T10 and F11 from SST chain P, respectively. All these interactions involve hydrogen bonding (Fig 7H). Strong correlations that involve hydrophobic contacts are less frequent in the Aβ 40 -SST 14 system than in the Aβ 42 -SST 14 system, evidently due to the lack of two C-terminal residues I41 and A42 in Aβ 40 .

Solvent accessibility of the aggregates
The solvent accessible surface area (SASA) is an important quantitative indicator of the oligomerization process. Fig 8 shows the evolution of total SASAs for Aβ chains and small cyclic peptides for three representative systems Aβ 42 -SST 14 -I, Aβ 42 -AVP-II, and Aβ 40 -SST 14 -II during 500 ns production MD simulations. Consistent with the ongoing aggregation in these systems, total SASAs of Aβ chains decrease significantly in the course of the simulations (Fig 8A). Although the greatest decrease occurs over the first 200-250 ns, a tendency of continuing collapse persists over the entire 500 ns interval. Total SASAs of small cyclic peptides also exhibit a steady decrease (Fig 8B). At the beginning of the simulations SASAs of SST 14 are close across the Aβ 42 -SST 14 -I and Aβ 40 -SST 14 -II systems; however subsequently SST 14 tends to exhibit a greater SASA in the Aβ 42 -SST 14 -I system than in the Aβ 40 -SST 14 -II system. There are two reasons for that. First, one or more SST 14 molecules always remain free in the Aβ 42 -SST 14 -I system. Second, from our analysis of the trajectories it follows that the SST 14 molecules tend to occupy peripheral regions of aggregates in Aβ 42 -SST 14 system, whereas in Aβ 40 -SST 14 system, they are often found in the interior of aggregates. The total SASA of AVP is less than that of SST 14 in either system due to smaller size of the AVP molecules. Table 4 lists average hydrophobic and hydrophilic SASAs for all the peptides in each of the three systems Aβ 42 -SST 14 , Aβ 42 -AVP and Aβ 40 -SST 14 during the first and the last 5 ns of 500 ns long MD simulations, and S5 Table provides similar data for SST 14 and AVP molecules separately. Additional details of the SASAs in individual trajectories are given in S6 Table. As one could expect, the hydrophobic SASA is less than the hydrophilic counterpart in all systems at all stages of the production MD simulations, since the aggregation process facilitates burying of hydrophobic groups. As Table 4 indicates, exposures of both hydrophobic and hydrophilic groups to the solvent have decreased considerably in the course of the simulations. The relative decrease amounts to approximately 28-37% of the initial SASA's values in each of the three systems. Upon averaging over three trajectories for each of the systems Aβ 42 -SST 14 , Aβ 42 -AVP  and Aβ 40 -SST 14 , the total hydrophilic SASAs are quite uniform across all three systems during the first 5 ns of MD simulations, with a slightly increased contribution of hydrophilic residues in Aβ 42 -AVP (Table 4). In contrast, the initial hydrophobic SASAs of the three systems exhibit diverging trends across the three systems. At the beginning of MD simulations, total hydrophobic SASAs in all Aβ 42 -SST 14 trajectories are higher than in any of Aβ 40 -SST 14 trajectories (S6 Table), which can be explained by the presence of the two hydrophobic C-terminal residues in Aβ 42 . On the other hand, the initial hydrophobic SASAs in both SST 14 containing systems are pronouncedly higher than in all Aβ 42 -AVP trajectories (Tables 4 and S6), which we attribute to the larger hydrophobic SASA of SST 14 molecules compared to AVP molecules (S5 Table).
In the course of the MD simulations, the Aβ 42 -SST 14 system exhibited the smallest relative decrease in hydrophobic SASA (33%), consistent with a smaller size of aggregates and greater solvent exposure of Aβ chains and SST 14 molecules in these systems ( Table 4). The reduction in hydrophobic SASA in both Aβ 42 -AVP and Aβ 40 -SST 14 systems is approximately 37% despite a relatively weak hydrophobicity of AVP. Together with a pronounced reduction in exposure of hydrophilic residues in Aβ 42 -AVP, this may suggest more efficient collapse in Aβ 42 -AVP system in comparison to Aβ 40 -SST 14 system.
Average hydrophobic and hydrophilic SASAs calculated for SST 14 alone (S5 Table) exhibit consistent trends with those in Table 4. Initially, SST 14 's hydrophobic and hydrophilic SASAs are close in both systems. In the course of the subsequent simulations, the decrease in both hydrophilic and hydrophobic SASA of SST 14 is less pronounced in Aβ 42 -SST 14 system than in Aβ 40 -SST 14 system. In the case of AVP, the relative decrease in hydrophobic and hydrophilic SASA of 49% and 43%, respectively, is the greatest of all systems. In addition, the high initial hydrophilicity in AVP results in a three times greater absolute reduction of hydrophilic SASA in comparison to the hydrophobic one.
Summarising, the greater initial hydrophobicity of the Aβ 42 -SST 14 system might result in substantial hydrophobic collapse; however, the observed decrease in solvent exposure of hydrophobic residues in the Aβ 42 -SST 14 system was relatively weak. Contrary to the expectations, a higher decrease of hydrophobic SASA during the MD simulations was observed in Aβ 40 -SST 14 system instead. The Aβ 42 -AVP system exhibited the strongest reduction in solvent exposure of hydrophilic residues, suggesting a more important role of electrostatic interactions in the course of aggregation in these systems.

Discussion
We have conducted all-atom MD simulations of early aggregation events at the onset of oligomerization in mixtures of Aβ peptides and small cyclic peptides, Aβ 42 -SST 14 , Aβ 42 -AVP, and Aβ 40 -SST 14 in water. Initially, the Aβ peptides and small cyclic peptides were sparsely positioned with respect to each other. Aggregation of the Aβ peptides, accompanied by a pronounced decrease in solvent exposure of both hydrophobic and hydrophilic groups, went on throughout the 500 ns long MD runs. In the course of the simulations, all systems developed self-assembled aggregates containing Aβ chains and small cyclic peptides. Consistent with published experiments [10,11] and MD simulations [15,16] of Aβ aggregation in the absence of small cyclic peptides, the self-assembled aggregates observed in this study exhibit largely unstructured morphologies without pronounced alignment of the chains. Although such unstructured morphologies indicate a considerable influence of hydrophobic collapse onto the oligomerization, the observed build-up of hydrogen bonding accompanied by development of stable β-sheet content also suggests a potential tendency toward β-conversion.
Overall, self-assembled Aβ 42 -SST 14 aggregates are compact and well-correlated dynamically. However, they tend to be smaller in size than aggregates formed in the other systems; their hydrophobic groups exhibit the highest solvent accessibilities in average by the end of the 500 ns simulation; and they develop less hydrogen bonds across Aβ 42 peptides and SST 14 molecules in comparison to the other systems. The replacement of SST 14 with AVP molecules produced larger aggregates. In addition, the Aβ 42 -AVP system tended to develop more hydrogen bonds than the other systems, which is especially true for cross-species hydrogen bonding. In the course of the 500 ns simulations, the AVP molecules developed almost 3 times more HBs with Aβ 42 than the SST 14 molecules did. Both hydrophobic and hydrophilic SASA are the lowest in the Aβ 42 -AVP systems in average after 500 ns simulations. The two Aβ 42 -containing systems exhibited similar tendencies of Aβ's C-termini to remain at the aggregate's surface, although C-terminal flexibility tended to be greater in AVP-containing systems than in SST 14containing systems. The most pronounced differences between the Aβ 42 -SST 14 and Aβ 40 -SST 14 systems concern the morphologies of the self-assembled aggregates. Unlike the Aβ 42 -SST 14 system, the Aβ 40 -SST 14 system formed large aggregates that contained a domain of correlated motion surrounded by relatively weakly correlated parts. In the Aβ 40 -SST 14 system the C-termini are buried more often than in the other two systems, and they exhibit the lowest average C-terminal flexibility. According to a recent report [21], such morphological differences may lead to different pathways at later stages of oligomerization.
Our examination of inter-molecular contacts that exhibited the strongest dynamics correlations in average over the last 20 ns of the simulations revealed differences in the interactions of SST 14 and AVP molecules with Aβ peptides. In Aβ 42 -SST 14 aggregates, highly correlated regions tended to involve residues K4, N5, F7, W8, and T10 of SST 14 , in agreement with biochemical observations [28]. These residues were found to interact primarily with the residues 28-35 of Aβ 42 , although contacts with its C-terminus were also detected (K4-V40 and W8-V40). The interactions tended to involve hydrophobic contacts, although backbone-backbone hydrogen bonds were also present. In the Aβ 42 -AVP aggregates, strong correlations involved predominantly contacts between N-terminal AVP resides C1, Y2, and Q4, and residues 20-23 of Aβ 42 . In the Aβ 40 -SST 14 aggregates, T10 and F11 of SST 14 were found to interact with C-terminal residues 35-40 of Aβ 40 . Hydrogen bonds are found frequently between strongly correlated residues in systems Aβ 42 -AVP and Aβ 40 -SST 14 .
At a coarser level, we observe a tendency of SST 14 molecules to develop hydrophobic contacts with surfaces of Aβ 42 aggregates such that both N-and C-terminal arms of SST 14 face the solvent. This may prevent hydrophobicity-driven coalescence the Aβ 42 aggregates, explaining the relatively small size of aggregates in the Aβ 42 -SST 14 system. Distinct from SST 14 , AVP molecules tend to penetrate inside the Aβ 42 aggregates, where they develop extensive hydrogen bonding with central regions of Aβ 42 peptides. When SST 14 interacts with small Aβ 40 aggregates, the morphologies resemble those of Aβ 42 -SST 14 , whereas large Aβ 40 aggregates tend to incorporate SST 14 molecules resembling the interaction of AVP with Aβ 42 aggregates.
Based on the results obtained in this study, the experimentally observed differences in the aggregation behaviours of Aβ 42 -SST 14 , Aβ 42 -AVP, and Aβ 40 -SST 14 mixtures [28,29] may be attributed to the different interaction of the small cyclic peptides with each of the Aβ peptides and the solvent.
First, the tendency of SST 14 to develop predominantly hydrophobic contacts with residues 28-35 and the C-terminus of Aβ 42 favors morphologies where several Aβ 42 chains form small aggregates compacted by Aβ-Aβ hydrogen bonding, and carrying SST 14 molecules attached to their surfaces. This surface attachment prevents the coalescence of small Aβ 42 aggregates into bigger structures. Within the time regimes of the present MD study, this results in a prevalence of multiple small Aβ 42 aggregates. The resulting delay in hydrophobicity-driven coalescence may explain the extended lag phase of Aβ 42 aggregation observed experimentally [28]. In contrast, AVP molecules preferentially bind to central regions of Aβ 42 peptides, and are prone to develop extensive inter-species hydrogen bonding. As a result, large Aβ 42 aggregates selfassemble quickly and contain inserted AVP molecules that strengthen connections within the aggregates. The binding of SST 14 to Aβ 40 tended to rely on strong connections with the C-terminal region of the Aβ 40 chains. Combined with the propensity of the Aβ 40 -SST 14 system to develop inter-species hydrogen bonding, this facilitated insertion of SST 14 into the interior of Aβ 40 aggregates.
Second, differences in hydrophobicity of the three systems result in different surface reactivities of the aggregates. Clearly, the presence of SST 14 imparts the aggregates with a greater hydrophobicity. The resulting "sticky surface" effect is enhanced even more due to the location of both SST 14 molecules and hydrophobic C-termini, preferentially at the surface of the aggregates. In the Aβ 42 -AVP system hydrophobic C-termini are also exposed and flexible. However, the "sticky surface" effect is lessened by a lower hydrophobicity of AVP and a lesser surface area due to the larger size of the aggregates. In the Aβ 40 -SST 14 system, the insertion of SST 14 molecules into the interior of the aggregates combined with the reduced number of exposed hydrophobic C-terminal residues plays a role.
Summarising, our observations suggest that mixtures of Aβ 42 and SST 14 , although pronouncedly more hydrophobic than their Aβ 42 -AVP and Aβ 40 -SST 14 counterparts, tend to aggregate slower than the other two mixtures retaining extensive "sticky surface" areas exposed for longer times. This may explain the increased pathogenicity of Aβ 42 , but not Aβ 40 peptides, in the presence of SST 14 , but not AVP molecules.

Modeling structures
Each system investigated in this study contained eight randomly positioned monomeric Aβ 42 or Aβ 40 peptides and eight randomly positioned small cyclic peptides, somatostatin-14 (SST 14 ) or arginine vasopressin (AVP), as listed in Table 5. For Aβ 42 peptides, the sequence 1 DAEFRHDSGYEVHHQKLVFFAEDVGSNKGAIIGLMVGGVVIA 42 was used, and for Aβ 40 peptides the two C-terminal residues were excluded. The sequences 1 AGCKNFFWKTFTSC 14 and 1 CYFENCPRG 9 -NH 2 were used for SST 14 and AVP, respectively, to match the conditions of the experiments [28,29]. As starting coordinates for these monomeric units, we used the coordinates of chains A-F of the NMR-derived Aβ 1-42 hexamer from PDB ID 2NAO [36] and the coordinates of chains A-H of the Aβ 1-40 nanomer from PDB ID 2M4J [37], which were equilibrated before modeling. For small cyclic peptides, we used the coordinates of chain A of somatostatin-14 with S-S bond from PDB ID 2MI1 [38], and coordinates of chain B of vasopressin from PDB ID 1YF4 [39]. Using the Accelrys VS [40] and VMD [41] packages, we built the starting systems to contain eight monomeric Aβ 1-40/42 units from 2M4J/2NAO.pdb, and Table 5. List of the systems studied.

System
Details of preparation Number of MD trajectories eight SST14/AVP molecules by randomly rotating and shifting the units against each other. All peptides were randomly placed in the simulation box such that distances between their surfaces were longer than 4 Å. This resulted in sparsely positioned units with a concentration of approximately 4 mM for Aβ, SST 14 , and AVP each. Minimizations, equilibrations, and production MD simulations were carried out using the GROMACS v5.0.7 package [42]. The AMBER99SB-ILDN forcefield [43] was used for protein atoms. Then, we minimized each set of 16 randomly positioned units using a steepest descent algorithm for energy in vacuum, and added the solvent using the optimized general-purpose 4-point OPC water model [44], as well as counterions Na+/Cl− to electro-neutralize the systems. Subsequent solvent minimizations involved decreasing position restraints on nonhydrogen protein atoms, as well as heating with the Berendsen thermostats. NVT-equilibrations were performed on all the systems. To collect sufficient statistics on the systems' dynamics, we ran independent MD trajectories for each of the systems, as listed in Table 5. The same sets of starting coordinates of atoms were used in each of the trajectories, whereas the seed numbers were varied when generating the initial Maxwell distributions of atom velocities. The three MD trajectories were labeled as "I", "II" and "III".

Molecular dynamics simulations
The production MD simulations, as well as the last equilibration step, were conducted at a temperature of 310 K and a pressure of 1 atm using isotropic pressure coupling (NPT ensemble), the Verlet cut-off scheme for neighbour searching, the Particle-Mesh Ewald treatment of electrostatics, a twin-range cut-off for van der Waals interactions, and bond lengths restrained through the linear constant solver (LINCS) algorithm with a fourth order of expansion. A Vrescale scheme was used for temperature coupling, and the Parrinello-Rahman algorithm was employed for pressure coupling applying the scaling to the center of mass of reference coordinates. Production MD simulations were performed for 100 ns for 9 systems containing 16 chains each (see Table 5). 2 fs time steps were used, and snapshots were saved every 20 fs in order to analyze the essential collective dynamics of the systems.
Structural analysis of the trajectories, including assessment of the secondary structure, the number of hydrogen bonds, and solvent accessible areas, has been done using the GROMACS scripts [42] and the VMD package [41]. For graphical representation, the VMD and Accelrys VS packages were utilized.

Essential collective dynamics analysis
To analyze the dynamics of the aggregates in greater depth we employed the novel essential collective dynamics (ECD) method [16,[31][32][33][34][35], which stemmed from a recently developed statistical-mechanical framework [31,32]. In this framework, a macromolecule or supramolecular complex is described by generalized Langevin equations with a set of essential collective coordinates identified as principal eigenvectors of a covariance matrix, which in turn is calculated by applying the principal component analysis (PCA) on MD trajectories. It has been demonstrated [31,32] that persistent correlations between atoms' motion (dynamics) can be calculated from a projected all-atom image of the system in a multi-dimensional space of the essential collective coordinates, as described in details elsewhere [16,[31][32][33]. A suite of dynamics descriptors has been derived within this framework, including pair correlation maps [16,33,35], dynamics domains of correlated motion [16,31,33,34], and main-chain flexibilities [16,[33][34][35]. The method has been validated extensively against NMR [31][32][33][34]45] and X-ray [33,46] structural data, and the corresponding descriptors were demonstrated to predict accurately the persistent dynamics trends from short fragments of MD trajectories without exhaustive conformational sampling [31,32]. In this study, we applied the ECD method with a set of 10 principal eigenvectors, which were sufficient to sample more than 95% of the total displacement in the systems considered. The ECD descriptors were calculated from multiple 0.2 ns long segments of the MD trajectories, and averaged over pertinent segments of production MD simulations, as specified in Results and Discussion. The pair correlation descriptors displayed in Fig 4 have been derived from the projected images of Cα atoms in the respective aggregates as specified earlier [16,33,34,45]. The dynamics domains of correlation motion, which represent relatively rigid parts of the aggregate composed of atoms moving coherently (Fig 5) were identified by applying the clustering technique described in earlier reports [31,33,34,45] with a threshold of 0.008 to the projected all-atom images of all systems. The main-chain flexibility profiles (Fig 6) illustrate the level of dynamic coupling of main-chain atoms with the entire aggregate, which in turn is represented by the centroid of the projected all-atom images of the aggregate.  14 -II (C) during the 500 ns long production MD simulations. β-sheets/bridges are shown with yellow/dark yellow color, α/3-π-helices-with purple/blue, turns-with green, and random coils-with white. The Y-axis represents residues in chains A-P. chains C, D and F are colored cyan; and chains G and H are colored purple. In (B) and (C), chains A to H are colored yellow, cyan, purple, lime, mauve, ochre, iceblue and black, respectively. All SST 14 and AVP molecules are colored red. (TIF) S1 Table. Number of hydrogen bonds between and across Aβ chains and small cyclic peptides (SCP), and total number of hydrogen bonds in each system in average over the first and the last 20 nanoseconds of 500 ns long MD runs. (XLSX) S2 Table. Population of β-bridges in all trajectories averaged over the last 20 ns of the 500 ns long MD simulations, their averages over three trajectories for each system, the averages over the entire 500 ns simulations, and standard deviations of the data. The data are given in % of occupied Aβ residues. (XLSX) S3 Table. Population of α-helices and 3-π-helices in all trajectories averaged over the last 20 ns of the 500 ns long MD simulations, their averages over three trajectories for each system, the averages over the entire 500 ns simulations, and standard deviations of the data. The data are given in % of occupied Aβ residues. (XLSX) S4 Table. ECD main-chain flexibility of C-terminus in each Aβ 42