Experimentally Consistent Simulation of Aβ <sub>21–30</sub> Peptides with a Minimal NMR Bias

■ INTRODUCTION

Amyloid fibrils are self-assembled fibrillar peptides 1 which are linked to Alzheimer's, 2 Parkinson's, 3 Huntington's, 4 prion diseases 5 and type II diabetes. 6 In Alzheimer's disease (AD), soluble amyloid β (Aβ) peptides misfold and self-assemble into neurotoxic insoluble fibrils. 6 Aggregation of amyloid monomers may fit classical nucleation theory: 7 aggregation is initially dictated by the rate of critical nucleus formation. The formation of a critical nucleus from Aβ monomers is called a primary nucleation pathway 8 and depends only on monomers. 9 In the secondary nucleation pathway, the fibrils elongate as monomers attach to the nucleation site, which depends on both monomer and fibril concentrations. 9 A fibril-catalyzed secondary pathway seems to be a driving force in further aggregation of both Aβ40 and Aβ42, the two forms of the amyloid peptide. 8,9 In order to understand how monomers form nuclei during the primary nucleation pathway, one needs to understand the structural and dynamical properties of the monomers. Molecular understanding of both the primary nucleation and secondary nucleation pathways of aggregation is essential to finding the routes that halt aggregation and prevent the progression of fibrillogenesis.

Although some recent work on the mechanism of AD has focused on neuroinflammation, 10 the propensity of Aβ is still correlated to neurotoxicity, and its understanding is relevant for Alzheimer's and other prion diseases. Familial mutants of the Aβ42 peptide that increase its aggregation were also found to be more neurotoxic. 11 When monomers interact with fibrils and start forming small oligomers, they become more toxic. 8,9 Oligomers of various sizes have been found to be more toxic than mature fibrils. 5,[12][13][14] This shows that early stages of aggregation, difficult to study with experiments but accessible with simulation, may be important for understanding the disease.

The aggregation process of Aβ has been carefully studied in the past via both simulations and experiments. Lazo et al. 15 used proteolytic cleavage (at nondenaturing conditions) to study both Aβ40 and Aβ42 and to show that there is a protease-resistant β turn involving residues between V24 and K28. They hypothesized it as a nucleation site for aggregation and subsequent fibrillogenesis. Experimental groups [15][16][17] and theory groups 18,19 have investigated the Aβ 21-30 fragment to further support this region as a nucleation site. The salt bridges, E22-K28 and D23-K28, in the protease-resistant region can withstand denaturing conditions, such as high temperature. 20 It is not clear whether the salt bridges (D23-K28 and E22-K28) or the hydrophobic interactions of V24 with hydrophobic parts of K28 lead to the resistance to denaturing conditions in this fragment. Based on the importance of the Aβ 21-30 region in aggregation behavior, this fragment was chosen for study.

The most commonly observed conformations of the Aβ 21-30 fragment prior to fibrillization are random coil and β turn conformations, according to experiments and simulations. 15,[20][21][22][23] Simulations that predict the β turn to be a prominent structure also suggest the importance of the D23 and K28 salt bridges for stability and formation of metastable conformations. 19,[24][25][26][27][28][29] Understanding how a random coil conformation of intrinsically disordered proteins/peptides (IDPs) changes to β hairpins is important in halting aggregation. Mutagenesis simulation of Aβ 21-30 fragments found that D23 is important in stabilizing the β turn by adding Ca 2+ , which destabilized and disrupted the salt bridges in the β turn. 26 Phosphorylation of S26 in the 21-30 fragment decreased the ability to form bends/turns by removing the salt-bridge between D23 and K28 and decreasing the plasticity of the fragment. 21 Studying peptide interactions is important in possibly using Aβ-derived peptides as inhibitors of aggregation. Molecular dynamics (MD) simulations of Aβ42 helped to understand how peptide oligomers bind to a surface (self-assembled monolayers) and understand their behavior at the interfaces. 30 Other peptides have been investigated as a mode of reducing neurotoxicity and fibrillogenesis of amyloid peptides. 31,32 The heptapeptide (Aβ 16-22 ) component of the Aβ42 peptide inhibited fibrillogenesis and reduced cytotoxicity. 32 When conjugated with larger cyclo-β-dextrin, the heptapeptide's inhibitory and anticytotoxicity effects were improved. A different set of short peptides derived from Aβ42, conjugated in tandem on the surface of gold nanoparticles, inhibited fibrillogenesis and reduced the toxicity of the gold nanoparticles. 31 Therefore, peptide interactions are potential routes to halting fibrillogenesis.

In this work, we simulate Aβ 21-30 fragments using a method that ensures MD simulations' agreement with experiments by incorporating information from experimental measurements into the simulation. 33 NMR experiments and simulations were performed to investigate the structure and dynamics of Aβ 21-30 . Experiment directed simulation (EDS) was used to simulate Aβ 21-30 with high accuracy by using experimentally measured chemical shifts (δ H ) as inputs to the simulation. 34 EDS is a minimal biasing technique that matches the ensemble average of simulations to experimental biases. 34 Bias can be added along any position-dependent collective variable (CV), such as δ H , radius of gyration (R g ), and so forth. Understanding Aβ peptide aggregation requires understanding its dynamics, and EDS has been shown to work well in improving thermodynamic accuracy while also maintaining or improving accuracy of dynamics. 35 To address the large free energy barriers of peptide-peptide interactions, the enhanced sampling technique parallel tempering well-tempered ensemble (PT-WTE) was used. 36,37 PT-WTE was chosen because it requires an order of magnitude fewer replicas than parallel tempering simulations to achieve a comparable replica exchange efficiency. 36 Incorporating experimental results into simulations of peptides and proteins has been performed in the past to improve the accuracy of simulations. [38][39][40] For rigid proteins, that do not change their secondary structure easily, chemical shift MD (CS-MD) method that uses δ H as restraints has been used to improve the accuracy of simulations and improve agreement with experimental results. 39 CS-MD was able to reproduce the X-ray results of proteins that are rigid and have a secondary structure. Aβ 21-30 is a dynamic IDP because it does not have a permanent secondary structure, and thus, EDS was chosen as it provides a maximum entropy guarantee.

■ MATERIALS AND METHODS Experimental Details. Aβ 21-30 peptide was purchased from GenScript USA Inc. The peptide purity was 96.2%, as specified by the supplier. D 2 O was purchased from ACROS Organics. NMR samples were prepared by dissolving 0.5 mM peptide in 90% H 2 O/10% D 2 O containing 25 mM sodium phosphate (pH 6.05) and 25 mM NaCl. The watergate pulse sequence was used to suppress the water signal. All NMR experiments were conducted in a Varian Inova 600 MHz spectrometer. One-dimensional and two-dimensional (2D) 1 H homo-nuclear NMR spectra were obtained at 278 K. The 2D total correlation spectroscopy (TOCSY) spectrum was collected using 57 ms mixing time. The 2D NOESY spectrum was obtained with 200 ms mixing time to investigate spatial proton connectivities. Because NOESY intensities of Aβ 21-30 were low because of the intermediate tumbling time, rotating frame Overhauser effect spectroscopy (ROESY) was conducted with a same mixing time of 200 ms. The TOCSY and NOESY/ROESY spectra were used to complete the proton resonance assignment. NMRPipe and NMRFAM-SPARKY 41 were respectively used for data processing and data analysis.

A pulse gradient spin echo (PGSE) diffusion study was conducted to measure the diffusion coefficient of the Aβ 21-30 monomer. In PGSE, [42][43][44] a series of spectra is recorded with a progressively increasing pulse gradient amplitude. The diffusion coefficient can then be calculated by fitting the resulting peak decay to the Stejskal-Tanner equation. 45 Simulation Details. MD simulations were used to study the Aβ42's2 1 -30 (AEDVGSNKGA) region. One, two, and four monomers of the fragment were simulated with 0.013 g/ cm 3 peptide density. This peptide density was used to better replicate high peptide concentrations used in NMR experiments (0.5 mM in our experiments). We closely followed the method developed for using EDS with enhanced sampling described in Amirkulova and White. 37 EDS minimally biases the potential energy such that a desired experimentally observable CVs converge to the experimental values. 34 EDS matches multiple ensemble averaged experimental observables in one iteration of the simulation. EDS is a type of maximum entropy technique which is used to incorporate experimental data into simulations. It is one of the many recent approaches to minimal biasing simulations to better match experimental

The Journal of Physical Chemistry B data. 34,37,[46][47][48][49][50] A recent review of EDS can be found in the work by Amirkulova and White. 33 Given a CV, f(r⃗ ), EDS computes a potential energy bias U′(r⃗ ,f(r⃗ )) and adds it to the potential energy of the whole system, U(r), such that the ensemble average of f(r) matches a chosen value: ⟨f(r)⟩ = f̅ . Here, f̅ is measured from an experiment. 34 Using the stochastic gradient descent algorithm, EDS matches experiments by computing coupling constants, α i , that control how strongly each biased ith CV affects the potential energy. 34 δ H ij is a proton chemical shift where i refers to one of the ten experimental proton chemical shifts and where j iterates over to the number of monomers, K.

(1)

All simulations were performed using the Amber ff99SB 51 force field and the TIP4P-Ew 52 water model. The ff99SB/ TIP4P-Ew combination does not overpredict the ROESY spectrum of the Aβ 21-30 fragment, 22 and this force-field combination is used widely with NMR simulations. 53 Simulations were performed in the NVT ensemble using the canonical sampling through velocity rescaling thermostat. 54 The Aβ 21-30 peptide was simulated in explicit water in a cubic box with a height of 48.6 Å for one-monomer, 61.9 Å for twomonomer, and 77.9 Å for four-monomer systems. One Na + was added to the system per monomer, giving approximately a 10 mM ionic strength, to neutralize any net charges in the simulation systems.

A time step of 2 fs was used to integrate the forces. Longrange electrostatic forces were calculated with the particle mesh Ewald method. Shifted van der Waals and short-range electrostatics were used with a cutoff distance of 10 Å. Hydrogen containing covalent bonds were constrained using the LINear Constraint Solver (LINCS) algorithm. 55 All simulations were performed in the MD simulator, Gromacs-5.1.4. 56 GromacsWrapper was used to automate Gromacs commands. 57 An initial extended structure of Aβ 21-30 was generated with PeptideBuilder. 58 Packmol was used to pack together multiple peptide chains in a single coordinates file (gro file). 59 The PLUMED2 package was used to calculate CVs (δ H , radius of gyration, torsional angles, and distances) and to bias the system via its implementation of the EDS method. 60 δ H is computed in PLUMED2 using a method called Camshift which was developed by Kohlhoff et al. 61 MDAnalysis was used when analyzing the contacts between peptides. 62 All systems used an energy minimization, 2 ns annealing and 2 ns equilibration at 300 K and 1 bar using the Parrinello-Rahman NPT thermostat. 63 Sixteen replicas of PT-WTE were then equilibrated without EDS for 400 ps to find a PT-WTE bias (note separate from EDS bias) that gave at least 20% replica efficiency. Then, 16 replicas of EDS simulation with a constant PT-WTE bias ran for 40 ns to find parameters of EDS that biased the system with experiments and sampled the equilibrium ensemble. All properties calculated during this simulation require reweighting to remove the effect of the PT-WTE bias and was done so for all analysis, unless otherwise noted. 36 The 16 replicas had the following temperatures: 278.0, 284.4, 290.9, 297.6, 304.5, 311.5, 318.6, 326.0, 333.5, 341.1, 349.0, 357.0, 365.2, 373.6, 382.2, and 391.0 K. Only the lowest temperature replica at 278 K was considered in the results and data analysis. The lowest temperature replica's temperature was 278 because the experiments were conducted at that temperature. PT-WTE parameters that had high replica exchange efficiency were found to be a Gaussian hill height of 2 kJ/mol, width of 140 kJ/mol, biasfactor of 10, and period of 250 time steps. 36 The EDS minimal bias was updated at each replica exchange attempt, every 400 time steps. Ten δ H were used as CVs in biasing with EDS to match experimental δ H . For each δ H that was biased, the EDS force constant had a range of 1 k b T. Each EDS CV was updated with a probability of 0.1 at each time step.

Control simulations were performed in the NVT ensemble at 278 K for 40 ns for structural property calculations. To test how MD computes the properties of Aβ 21-30 , we performed a plain MD (non-EDS and non-PT-WTE) simulation for a onemonomer system, labeled as MD (1) in Table 1. To test how the lack of EDS during PT-WTE simulations affects the conformations, we ran a control simulation of one monomer with the same PT-WTE setup and without EDS bias, labeled as PT (1) in Table 1. We also ran a four-monomer control system that was non-EDS but had the PT-WTE bias which is labeled as PT (4)

where k b is the Boltzmann's constant, T is the absolute temperature, η is the viscosity of the solvent, t is the time, r(t) is the position at time, and r H is the hydrodynamic radius of a molecule. Diffusion coefficients for control (non-EDS and non-PT-WTE) simulations are calculated in the NVE ensemble also for 40 ns, which is labeled as Diffusion MD. We also ran 500 ns production EDS simulations for one monomer, labeled Production EDS (1) in Table 1, and two monomers, labeled Production EDS (2) in Table 1. These production simulations were run using the EDS biases from equilibration steps in the NVE ensemble to thoroughly analyze the structure. The distance between centers of mass of chains over time was used to assess convergence of the simulations. Multiple binding and unbinding events during the simulation are indicators of good conformational sampling of the peptide interactions, as discussed in the Supporting Information. The total simulation of the replica exchange system was 640 ns including all replicas. Uncertainty for calculated NMR δ H from simulation was assessed via the block averaging method with a block size of 1 ns. 65 ■ RESULTS AND DISCUSSION Experimental Results. Figure 1a shows the TOCSY spectrum of Aβ 21-30 . It does not show cross-peaks between the lysine side chain 1 H ζN and its other side chain 1 H such as 1 H β , 1 H γ , and 1 H δ (denoted by K28HB, K28HG, and K28HD in NMR spectral images). The H N -Cα fingerprint region of the TOCSY spectrum of Aβ 21-30 is included as Figure S1 in the Supporting Information. This suggests that K28 is mobile with the side chain 1 H N being exchanged rapidly with the solvent. The result is consistent with the hypothesis that solvation of monomers hinders salt bridge formation and hence prevents formation of stable structures like the β turn in the Aβ 21-30 region. 66 The corresponding ROESY spectrum of Aβ 21-30 is given in Figure 1b. It exhibits sequential inter-residue d αN (i,i + 1) crosspeaks with relatively higher intensities than intraresidue d αN (i,i) cross-peaks, indicating a random coil structure with a stretched backbone. [67][68][69][70][71] Figure 1b does not show long-range interactions, such as the E22 1 H α -A30 1 H N interaction reported in a previous study 17 on Aβ 21-30 at 10 °C by Grant et al. 17 Later, a study by Fawsi et al. 22 using 900 MHz NMR showed a well-resolved spectrum where the previously reported E22 1 H α -A30 1 H N ROE cross-peak was missing. This discrepancy was attributed to misassignment of peaks which are too close to be distinguished by a lower NMR field. Grant et al. 17 also reported visible TOCSY cross-peaks between K28 H ζN and K28 1 H γ and K28 1 H δ and K28 1 H ϵ . These results were attributed to the presence of a β turn structure between V24 and K28 stabilized by electrostatic interactions between K28 and E22 or D23. However, Fawsi et al. 22 did not observe long-range ROE cross-peaks; only weak medium-range ROE cross-peaks were detected. A more recent study 67 of Aβ40 and Aβ42 also agreed with Fawsi et al. 22 and detected no long-range E22 1 H α -A30 1 H N ROE interaction. Moreover, the study did not find any significant evidence to suggest high β turn-forming propensity in the G22-K28 region. It should also be noted that buffers used in a study can largely affect the conformations and fibril formation. 72 A number of previously reported Aβ studies were conducted using ammonium acetate buffer. 15,17,22,29 Like the present study, Roche et al. 67 had used phosphate buffer.

The experimental value of the diffusion coefficient of Aβ 21-30 was calculated to be 1.8 × 10 -6 cm 2 /s by fitting PGSE data to the Stejskal-Tanner equation. To our knowledge, the diffusion coefficient of Aβ 21-30 at 278 K has not been reported. Danielsson et al. 73 measured the translational diffusion coefficient of Aβ40 and 14 shorter Aβ sequences, some with point mutations, at 298 K and at different concentrations (50 The Journal of Physical Chemistry B μM to 3.12 mM) and pH (6.33-7.47). Previously reported values of the diffusion coefficient of Aβ40 are ∼1.52 × 10 -6 cm 2 /s 73 and ∼1.2 × 10 -6 cm 2 /s 71 at 298 K and at 293 K, respectively. In a recently published work, the diffusion coefficient of Aβ42 at 298 K was reported to be ∼1.35 × 10 -6 cm 2 /s. 74 The diffusion coefficient of the monomeric fragment Aβ 12-28 has been found to be ∼1.09 × 10 -6 and ∼2.18 × 10 -6 cm 2 /s at 298 and 274 K, respectively. 75 Therefore, our measured diffusivity is consistent with previously measured Aβ fragment monomer diffusivities.

Simulation Results. EDS significantly improves agreement with experimental chemical shifts. Figure 2 shows the absolute difference between the cumulative average backbone δ H and experimental backbone δ H over the simulation period for the different simulation conditions. The experimentally measured δ H values are shown in the figure legend of Figure 2a. The horizontal line in Figure 2 is the standard deviation of the absolute difference of all δ H with respect to their experimental values. The EDS converges to agreement within this noise term within 40 ns. The average differences, δ H , are 0.07 and 0.03 ppm for amide and α hydrogens, respectively. The PT (1) simulation had average differences of 0.4 and 0.2 ppm for the amide and α hydrogens, respectively. The uncertainty for experimental δ H is estimated to be ±0.02 ppm based on the ambiguity of peak assignments from ROESY and TOCSY spectra. This is significantly smaller than the error from the predictive NMR model Camshift and thus is not considered in analysis. The uncertainties in Camshift are estimated to be 0.46 ppm for N H and 0.23 ppm for H. 76 Thus, the agreement is within the experimental uncertainty, sampling uncertainty, and NMR model uncertainty. One may argue that the EDS simulation is overfit; however, maximum entropy methods are designed to not overfit by matching averages and provably maximizing the entropy of the chemical shift distribution.

Two-[EDS (2)] and four-monomer [EDS (4)] EDS simulations with PT-WTE were performed to study interpeptide interactions. When there is more than one monomer, each backbone δ H is averaged over the number of monomers, as shown in eq 2. EDS significantly improves the agreement of calculated δ H to the experimental δ H for all monomer number a,c,d). The absolute difference in δ H in EDS simulations is smaller than the uncertainty computed from block averaging after the first 3.2 ns. Uncertainties of the predicted chemical shifts (Camshift) are estimated to be 0.46 ppm for N H and 0.23 ppm for H. 76 The outlier, H 4 , in (d) EDS ( 4) is discussed in the Supporting Information in Figure S2.

The Journal of Physical Chemistry B systems. The δ H of H 4 (red line) in four-monomer EDS simulation (Figure 2d) was an exception, which is discussed in Figure S2 in the Supporting Information.

Structure of Fragments. EDS allows a single monomer simulation to give insight into the aggregating structure because the experimental data are measured with aggregating peptides. This is shown in Figure 3 which shows the free energy surface for torsional angles of E22. E22 is important because it can form salt bridges, and E22Q mutations make the peptide more flexible and prone to aggregation. 77 E22 with EDS has the same local (L-α) and global minima (β) with one or four monomers in Figure 3. The δ H information, which is obtained with multiple monomers from experiments, biases the single monomer simulation with this information about the multiple monomer system. This concept of using EDS to incorporate information from a more complex system has been explored as "coarse-grained directed simulation" by Hocky et al. 78 Separately, the IDP character of this fragment can be observed by the similar free energy of both the α helix (D-α) and β regions.

The PT (1) results (Figure 3b) show that the change in the EDS secondary structure represents a change in free energy and is not a consequence of improved sampling. Even with additional sampling, PT never observes the two local minima in L-α. When we look at δ H of PT (4) in Figure S7 in the Supporting Information, the disagreement between experimental and computed δ H is greater compared to EDS (4). Thus, EDS improves the agreement between experimental and computed δ H not by biasing to explore new conformations but by improving the force field accuracy.

Convergence of our simulations was assessed in three ways. Note that the total simulation time of all replicas is 640 ns, and PT-WTE is an enhanced sampling method. Figure S3 shows that there are multiple binding and unbinding events between monomers. Figure S4 shows the free energy surfaces for E22 at 20, 30, and 40 ns and shows the number of local minima converges. Finally, we ran two 500 ns long single replica simulations without PT-WTE for one and two monomers. These longer simulations miss local minima because of the inability to cross free energy barriers because they do not have enhanced sampling and do not find any other local minima. The results for these systems are discussed in more detail in the Supporting Information.

To better understand what conformations and what interactions correspond to the global minima, we performed cluster analysis and calculated contact maps as shown in Figure 4. Clustering an MD trajectory partitions the observed conformations, and the most frequently observed clusters correspond to free energy minimum structures. Cluster analysis was performed with Gromacs using single-linkage hierarchical clustering of the peptide fragments. 79 Figure 4 shows the most probable structures for EDS (1), MD (1), and PT (1) from cluster analysis on the top panel. For the contact maps, we counted how frequently nonhydrogen atoms were closer than 3.2 Å. The most prominent cluster from EDS shows a β turn in the D23-S26 region in Figure 4a. The rest of the fragment is The Journal of Physical Chemistry B in a random coil conformation, which agrees with Chen et al. 28 and Baumketner et al. 20 The cluster center in Figure 4a shows that the β turn is formed by S26, D23, and K28 hydrogen bonding. EDS simulation shows both K28 and S26 to be concave inward which agrees with the results showing that phosphorylation of S26 stopped aggregation in the work by The Journal of Physical Chemistry B Rezaei-Ghaleh et al. 21 Our EDS cluster also agrees with the cluster center by Rezaei-Ghaleh et al., 21 where they show a hydrogen bonding between K28 and D23 at 3.3 Å with S26 concaved inward. Chen et al. 28 found the E22-K28 salt bridge to be an important interaction in forming this β turn; however, our results showed the above interactions to be more significant.

The MD (1) simulation cluster has a β turn, although at a different location in Figure 4b. This location is inconsistent with both past work 21,28,66 and the EDS (4) and EDS (1) simulations. The PT (1) simulation in Figure 4e leads to this β turn disappearance, demonstrating that it is a sampling artifact. Thus, the force field neither captures the existence of the β turn nor the correct location, except by undersampling in a local minimum.

Dynamics. It has been previously observed that by improving accuracy of structural properties in EDS, dynamic properties also improve. 34,35 Thus, the diffusion coefficient (D T ) was computed along with the closely related R g , as shown in Figure 5. The global minimum of R g for Aβ 21-30 during all EDS simulations (one, two, and four monomers) and PT simulations was around 0.73 nm in Figure 5a at 278 K in the NVT ensemble. Only the one-monomer MD simulation had a larger global minimum for R g at 0.90 nm. Cruz et al. 19 studied The Journal of Physical Chemistry B Aβ 21-30 at 283 K in the NPT ensemble using the OPLS-AA force field and TIP4P water and found R g to be around 0.55 nm. They also predicted a β hairpin conformation, which is much more confined than a β turn conformation. This explains smaller R g .S m i t he ta l . 25 noted that the different conformations arise in Aβ 21-30 in different force fields and found that OPLS-AA gives more β hairpin conformation and Amber ff99SB gives more α helix conformation and both of them in TIP4P water give about the same probability of a β turn conformation. Smith et al. 25 measured the R g for Aβ 21-30 to be 0.55 nm in the NPT ensemble at 300 K using Amber ff99SB and TIP4P.

The presence of K28-D23 contact allowed the fragment during one-monomer EDS simulation to have a compact shape, whereas the MD simulations had a more extended conformation in the absence of both EDS and PT-WTE biases. As shown by the intramolecular contact map of the one monomer EDS in Figure 4d, EDS (1) arrives at a more compact structure with a lower R g via the D23-K28 and D23-S26-K28 interactions. This compact conformation results in higher D T . D T computed during EDS simulation, 1.7 ± 0.7 × 10 -6 cm 2 /s, was close to our experimental value of 1.8 × 10 -6 cm 2 /s. As mentioned before, we did not find previously reported values of D T for Aβ 21-30 . Massi et al. 80 studied a different fragment Aβ 21-35 using both experiments and simulations and measured D T and R g . The average D T using four trajectories was 1.4 × 10 -6 cm 2 /s, and the average R g was 0.92 nm at 300 K. 80 Because EDS on average enforces experimentally relevant conformations even on the one-monomer system that does not have neighboring chains in the vicinity to interact with, D T was higher for diffusion EDS which had only one monomer. Higher D T means the molecule becomes more confined, compact, and dynamic; hence, it has smaller R g . During early stages of aggregation under aggregation-prone conditions, the diffusion of monomers becomes slower in Aβ42 and Aβ40 monomers. 81 Similarly, at low salt concentrations, which is more aggregation-prone, the diffusion of proteins is slower compared to higher salt concentrations. 82 Because both of our experimental and simulation investigations are conducted at 278 K and at low peptide concentration, which is significantly below the critical aggregation concentration of 90 nM, 74 the diffusion coefficient increased when the EDS bias was added. This may mean that for one-monomer diffusion EDS, the conformation actually is less aggregation-prone than the onemonomer diffusion MD simulation.

Aggregation of Fragments. Backbone interactions between two monomers in the EDS (2) simulation and among four monomers in EDS (4) led to aggregating structures. We analyzed the two monomer system with cluster analysis with a cutoff radius of 2 Å and found that two monomers self-assemble into parallel dimers as shown Figure 6a. We analyzed the two-monomer system with cluster analysis with a cutoff radius of 2 Å and found that two monomers selfassemble into parallel dimers as shown Figure 6a. The two monomers were held by D23i-A21j and K28i-K28j intermolecular interactions between chains i and j. K28i-K28j interactions are between backbone atoms of chains i and j displayed in Figure 6a. The cluster center highlights whose interactions are intramolecular and intermolecular in Figure 6a with a black dashed line. The past work with a two-monomer system of Aβ 21-30 with REMD showed the cluster center to be in an antiparallel conformation, held together by the E22 and K28 interaction. 83 Smith et al. 83 noticed that the dimer was in a random coil conformation most often, followed by an unaligned dimer, antiparallel and least often in parallel conformations. When the number of monomers increased, the propensity to form a stable structure decreased, and random conformations were more favored. 83 EDS (2), which agrees with the NMR experiments, identifies the most prominent structure to be a parallel dimer with a coil secondary structure with a slight β turn.

When the number of monomers increases in the presence of EDS, the chains self-assemble into different orientations and form more β turns. The cluster center shown in Figure 6b) for EDS (4) is less aligned though still aggregating. Cluster analysis of PT (4) simulation did not result in an aligned structure, as shown in Figure 6c. EDS (4) does not have a specific orientation between chains in both the intermolecular contact map and cluster center. Both EDS (2) and EDS (4) rely on D23-A21 and N27-N27 interactions between chains to aggregate. The EDS (4) simulation also visited some antiparallel conformations as indicated by the A21-A30 interactions in the contact map (Figure 6d). The β turn region is stabilized by intramolecular interactions in EDS (4). The intramolecular contact maps of EDS (2) and EDS (4) show that there are similar interactions in the 23-26 region and a less frequent D23-K28 salt bridge. The EDS (2) and EDS (4) intramolecular contacts shown in Figure 6e,f are similar except the strength of D23-S26 interaction, which explains the higher number of β turns in the EDS (4) cluster center shown in Figure 6b.

The PT (4) cluster center has both a β turn region and a random coil region in each of the four chains in Figure 6c. PT (4) had few off-diagonal contacts, showing no β turn structure and strong D23-S26 intermolecular interactions (Figure 6f). The intermolecular interactions were antiparallel. The main difference between both EDS (2) and EDS (4) intramolecular contact maps and that of PT (4) shown in Figure 6i is the strength of D23-K28 contact. These results show that a change to the force field, via EDS, is necessary to see structured aggregation and that undersampling is not the main limitation in accuracy.