Before the REST2 simulation, we performed short conventional molecular dynamic simulations for the Aβ40 peptide in solution to acquire ththe onset and progression of ADe pre-equilibrated initial conformations of the peptide for REST2. The initial structure of the full-length Aβ40 peptide is obtained from the Protein Data Bank (PDB ID: 2lfm; Figure 1A) with a helical peptide determined by nuclear magnetic resonance (Vivekanandan et al., 2011). The sequence of the Aβ40 peptide compared with the Aβ42 peptide is shown in Figure 1B. Compared to Aβ42, Aβ40 only lacks the last two hydrophobic residues I41 and A42 at the C-terminal (Tomaselli et al., 2006; Vivekanandan et al., 2011). This peptide was then under the calculation of the protonation state through the H++ web server (Gordon et al., 2005; Anandakrishnan et al., 2012) and put into a rectangle box filled with 0.15 M NaCl and TIP3P (Jorgensen et al., 1983) water molecules. The system underwent a 100 ns cMD at 343 K, and the final conformation was used as the initial peptide structure for the REST2 simulations. The peptide was placed 5.01 nm center of mass (COM) distance above the lipid bilayer and filled with 150 mM NaCl and TIP3P (Jorgensen et al., 1983) water molecules in a rectangular box. The lipid bilayer was constructed in CHARMM-GUI (Jo et al., 2009; Lee et al., 2016) with each leaflet containing 18 cholesterols (CHOL), 18 1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-L-serine (POPS) lipids, and 54 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) lipids, with a total of 36 CHOL, 36 POPS, and 108 POPC in the system. This POPC/POPS/CHOL bilayer in a ratio of POPC: POPS: CHOL = 3: 1: 1 is identical to our previous work (Wang et al., 2022).

All the simulations were carried out employing the GROMACS software, version 2020.6 (Van der Spoel et al., 2005), with the CHARMM36m force field (Huang et al., 2017) and a time step of 4 fs with hydrogen-mass repartitioning (Hopkins et al., 2015). The temperature coupling was dealt with using the Nosé–Hoover method while the Parrinello-Rahman barostat (Parrinello and Rahman, 1981; Nos’e and Klein, 1983) was employed for pressure coupling. The periodic boundary condition was used in all directions with the electrostatic interactions using the particle-mesh Ewald method (Essmann et al., 1995). The length of all bonds was constrained by the LINCS algorithm (Hess, 2008). The distance cutoff for the Lennard-Jones potentials and electrostatic interactions were all 1.2 nm. We also performed a 5000-step energy minimization with the steepest descent algorithm followed by 1 ns NVT and 1 ns NPT equilibration before final production runs.

The REST2 method coupled with PLUMED version 2.7.5 (Tribello et al., 2014) was used herein to explore the binding mechanism between full-length Aβ40 peptide and lipid bilayer, which provided an enhanced sampling method analogous to REMD, whereas fewer replicas were needed to achieve equal performance (Wang et al., 2011). The following expression was utilized to determine the temperature distribution in our REST2 (Jo and Jiang, 2015). T i = T min ⁡ exp   i ⁡ ln T max / T min N r e p −   1 (1a)

Peptides and ions were dealt with as the “hot” region while water molecules and membranes remained cold. Exchanges between the adjacent replicas were attempted every 2 ps. Eighteen replicas were used with effective temperatures ranging between 343 K and 500 K. Each replica runs for 500 ns, resulting in a total of 9 μs simulation times for this entire work.

Trajectories were analyzed using in-house codes and the GROMACS built-in programs. A free energy landscape (FEL) describing the binding process was constructed along the number of contacts and peptide-bilayer distance. Contacts between the Aβ40 peptide and lipid bilayer were determined between any heavy atoms of the peptide and bilayer within 0.5 nm. Peptide-bilayer distance was defined as the perpendicular COM distance between the peptide and the bilayer. For a closer look, the FEL describing each basin in the binding process was also constructed along the root mean square deviation (RMSD) of the peptide with respect to its initial conformation and β-sheet content. Distances between residues and membrane surfaces were generated by computing the perpendicular COM distance between each residue and bilayer and further subtracting half of the membrane thickness (2 nm). Hydrogen bond (H-bond) was defined using the criteria of 0.3 nm donor-acceptor distance cutoff and 20° angle cutoff. Free energies were defined using the following expression (1). G = − kT ⁡ ln P (1b) where G is free energy, k is the Boltzmann constant, T is the temperature, and P is the probability of the conformations appearing in one bin. Intra-peptide contacts were computed using the 0.5 nm cutoff within the heavy atoms of residues and the contact probabilities were the average of all peptides in the ensemble. For the REST2 simulation, the first 200 ns trajectory of each replica was discarded to avoid initial transients, that is, only the 200–500 ns trajectory of each replica was utilized for analysis. All the secondary structure content calculation, H-bond analysis, and snapshot generation were carried out using the VMD software (Humphrey et al., 1996).

The REST2 method was used herein to explore the binding mechanism between Aβ40 and the bilayer. The lipid bilayer is the POPC/POPS/CHOL membrane used in our previous work (Wang et al., 2022). Similarly, this membrane was used here based on its characteristic of being extensively studied and close to the authentic cell membrane, as well as its simplicity (van Meer et al., 2008; Lemkul and Bevan, 2011; Lai et al., 2017; Banerjee et al., 2021). As our previous study revealed, Aβ42 peptides showed little interaction with the bilayer, and the obtained binding models were rare at the physiological temperature of 310 K, while the high temperature of 343 K could make the membranes ready to be penetrated without modifying their basic architecture (Wang et al., 2022). The interactions raised a lot at 343 K and the analyses of Aβ42 peptide-bilayer interactions were all based on the results at 343 K. Hence, with the aim of promoting Aβ40 binding to the lipid bilayer and for the convenience of direct comparison, we performed the REST2 simulation at 343 K for the Aβ40-bilayer system herein. The temperature distribution is calculated using Eq. 1a. Exchange rates between adjacent replicas were all larger than 17% (Supplementary Figure S1). The average exchange rate was 20.80% averaged over 18 replicas. The temperature trajectory for each replica sufficiently visits all the temperatures from 343 K to 500 K (Supplementary Figure S2). All the above outcomes demonstrated the sufficient sampling of our REST2 simulations. Consistent with our expectations, the peptides showed a high tendency to form interactions with bilayers at 343 K (Supplementary Figure S3). For most residues, the contact fractions were larger than 10% and the average contact fraction was 25.29%.

First, the free energy landscape describing the position and interaction of Aβ40 with respect to the bilayer is mapped (Figure 2A). The two-dimensional (2D) FEL is constructed based on the peptide-bilayer distance and the number of peptide-bilayer contacts. Not surprisingly, as the peptide approaches the bilayer, the interactions between the peptide and bilayer increase, leading to the FEL in an L-form. Seven free energy basins, a to g, are identified. Basin c corresponds to the lowest free energy in the 2D FEL, where ΔG = 0. By projecting the 2D FEL onto the peptide-bilayer distance, we obtain the one-dimensional (1D) free energy profile (FEP, Figure 2B). From the FEP, it can be found that Aβ40 is most populated in the region of 2.65 nm (e, f, and g) and 3.33–3.56 nm (a, b, c, and d) distance. These results indicate that Aβ40 energetically favors binding with the bilayer.

To further identify the bound structures distributed at the seven minima a to g characterized in Figure 2, we extracted the structural ensemble of the peptides at each minimum to construct the 2D FEL of the ensemble for each minimum to describe their structural features (Figure 3A). The FEL is mapped along the positional RMSD and β-sheet content of the peptide. As can be seen from Figure 3A, as the peptides are getting closer to the bilayer (from a to g), the β-sheet content is progressively reduced. Conformations at each minimum can be explicitly classified into two groups, and their representative structures are shown in Figure 3B. One group, labeled as a1, b1, c1, d1, e1, f1, and g1 of Figure 3B, contains bound peptides with predominantly β-sheet structure. Another group, labeled as a2, b2, c2, d2, e2, f2, and g2 of Figure 3B, contains bound peptides with predominantly helix or coil structure.

Snapshots in Figure 3B can accurately represent the microstates in Figure 3A except for a2. a2 microstate contains peptides with few residues interacting with the membrane surface and their conformations are disordered and distinct from each other. The remaining 13 microstates can be classified into four binding models according to their structural features and membrane-interacting regions. Distances between each residue and bilayer surface of the four interacting models are shown in Figure 4A. Snapshots of peptides and bilayers are depicted in Figure 4B. Microstates a1, b1, c1, and d1 are classified into model 1, where barely several residues of CT are inserted inside the bilayer with other residues remaining in the solution. Conformations in this model adopt a β-rich structure of 30%–50% β-sheet content with three parallel β-strands as NT, CHC, and CT region each possessing one β-strand. Microstates b2, c2, and d2 are classified into model 2, where several residues of NT are lying on the membrane surface. Peptides in this model adopt two helical fragments at NT and residues at these two fragments are slightly touching the bilayer surface. Model 3 contains microstates e1, f1, and g1, where the peptides are mainly β-sheet structures with most residues of CT buried inside the bilayer, and several residues in NT and CHC are lying on the membrane surface. Compared to the β-rich structure in model 1, peptides of this model also adopt the structure with three parallel β-strands, whereas the β-sheet structures are of 10%–30% content lower than model 1 with NT, CHC, and CL region, each possessing one β-strand. Besides, residues at CT in model 3 were buried significantly deeper than in model 1. Microstates e2, f2, and g2 are classified into the most deeply buried model 4. Peptides in this model are unstructured and dominated by a random coil with the CHC and CT buried below the membrane surface while NT and CL stretch into the solution. Two to 15 residues of peptides in model 1 and model 2 contact with the bilayer in the 3.0 nm–4.0 nm peptide-bilayer distance region, whereas 16 to 25 residues in model 3 and model 4 contact and insert deeper in the 2.5 nm to the 3.0 nm distance region.

These binding models of Aβ40 were further compared with those of Aβ42 reported in our previous work (Wang et al., 2022). The results show two significant differences in their binding models. From their representative structures in Figure 4, it can be seen that Aβ40 can adopt more diverse structures of bound peptide including β-strand, helix, or unstructured structure, whereas bound peptides in Aβ42 are mostly unstructured. This difference can be attributed to the highly enhanced sampling efficiency of REST2 compared to cMD as the latter was used in the work of Aβ42. Another difference is that residues at CT frequently interact with the bilayer in Aβ40 (Models one to four in Figure 4) while showing little tendency to bind with the bilayer in Aβ42. This can be explained by the two additional residues Ile41 and Ala42 in CT of Aβ42. Due to the absence of these two residues in Aβ40, two Val residues with stronger hydrophobicity are exposed at the CT region, which can interact with the bilayer more frequently driven by the hydrophobic interactions.

To give a deeper view of the conformation for the binding peptides, the H-bonds formed between the peptide and membrane/water for minima a to g were obtained (Figure 5A). From minima a to g, with the contacts increased, the number of H-bonds formed by Aβ40 with bilayer also shows a roughly rising trend. However, there is no obvious rising or declining trend for the H-bond formed with water. We then computed the H-bonds formed between binding models and the membrane/water (Figure 5B). The peptide-bilayer H-bonds exhibit no correlation with the peptide-water H-bonds. Conformations in model 2 form the highest number of H-bonds with the membrane but also show a high tendency to form H-bonds with water. Conformations in model 1 form the least H-bonds with the membrane. The average number of H-bonds formed by each residue with membrane was then acquired to explore the role of residues in peptide-membrane interactions in four binding models (Figure 5C). Two residues forming the highest number of H-bond with bilayer were pointed in each model. Most of these residues were charged especially for the deepest binding model of model 4, and positively charged K16 and K28 form strong H-bonds with the bilayer. Positively charged R5 residue at NT also plays a key role in forming H-bonds in models 2 and 3. We have labeled these residues intuitively in the structure of binding models in Figure 5D. They are all distributed in the peptide-bilayer interacting regions. The free energy data were computed using Eq. 1b to compare the relative stability of each microstate (Supplementary Table S1). Microstates d1, g2, and c1 showed the lowest free energy, indicating the high stability of these microstates. Free energies of other microstates were relatively high, especially for b2, c2, and d2 of model 2, suggesting the instability of this binding model.

Aβ peptides, as numerous works previously characterized, are intrinsically disordered in bulk water (Viola and Klein, 2015; Ono, 2018). The initial Aβ40 conformation used here is a partially folded structure with a helix structure spanning residues H13-D23. To explore the conformational transition of Aβ40 induced by binding to the lipid bilayer, an analysis of the secondary structure was carried out, as delineated in Figure 6. Aβ40 peptides were considered bound when their peptide-bilayer contacts in Figure 2A are larger than zero; otherwise, they are considered unbound. As can be seen from Figures 6A, B, for both the bound and unbound peptides, β-sheet (E), turn (T), and coil (C) structures were dominant. These three structures took up at least 79.50% and 90.80% of the bound and unbound peptides, respectively, and alternatively exist in each residue. In Figure 6C, the seven structures obtained were classified into three groups: the sum of the fraction for coil, turn, and bridge structure (C + T + B) was taken as the unstructured structure fraction; the fraction for β-sheet structure (E) was taken as the β-sheet structure fraction; the sum of 3–10 helix, α-helix, and Pi-helix structure (G + H + I) was taken as the helix structure fraction. Fractions of helix structure were low in both bound and unbound peptides, whereas residues E3-R5 and H13-Q15 at NT could form helix structures with fractions larger than 10% in the bound peptides. Peptides possessed more helix structure in the bound state than in the unbound state. β-sheet content (E) of both bound and unbound peptides was much higher than that of helix, especially for unbound peptides. Residues D7, E11, Q15, K16, L17, and F20 in unbound peptides adopted the β-sheet structure with fractions larger than 80%, which were all distributed at the NT and CHC regions. In the bound state, the highest fraction of β-sheet was 49.84%, which was evidently lower than the unbound state, and the residues with β-sheet structure were widely distributed in NT, CHC, CL, and CT regions. In summary, the β-sheet structure is dominant in both bound and unbound Aβ40 peptides.

To explore the driving force for the peptide binding, the peptide-bilayer contacts and H-bonds were analyzed. It is clear from Figure 7A; Supplementary Figure S4 that binding events occurred mainly in the NT, CHC, and CT regions. The contact fractions of residues in hydrophilic CL were significantly lower, especially for residues E22-N27. Residues at CT showed the highest tendency to interact with the bilayer. The top 10 residues with the highest contact fractions were G38, A40, G37, V39, Q15, V36, K16, H14, M35, and L34, which were all hydrophobic or positively charged residues. We compared the average fractions of residue-bilayer contact between Aβ40 and Aβ42 in the four sequence regions in Supplementary Figure S4 (Wang et al., 2022). It is remarkable that residues at CT in Aβ40 frequently interacted with the bilayer, whereas they showed little tendency to bind with the bilayer in Aβ42. From Figure 7B, R5 formed the highest number of H-bonds much larger than other residues, followed by K16 and K28 residues. These three residues also exhibited the same behavior in Aβ42 (Wang et al., 2022). The highest 10 types of H-bonds formed between R5, K16, and K28 residues and bilayer are provided in Supplementary Figure S5. The most populated H-bonds are formed with the side chain of POPC or POPS lipids. Like in the Aβ42-bilayer system, cholesterols also do not favor forming H-bonds with these residues (Wang et al., 2022). Overall, residues in the CT region of the Aβ40 peptide showed the highest tendency to interact with the bilayer, and H-bonds formed between positively charged residues and bilayer may drive the binding. Both hydrophobic interactions and electrostatic interactions contributed to the binding of Aβ40 to the bilayer.

Besides peptide-bilayer interactions, intra-peptide and peptide-water interactions may also drive the peptide bind to the bilayer. Residue-residue contact maps for the bound and unbound Aβ40 peptides were plotted to investigate the intra-peptide interactions (Figure 8A). Unlike Aβ42, the conformations of Aβ40 in the unbound state were rather different from the bound state (Wang et al., 2022). The evident cross-diagonal formed by residues F4-V12 and residues H14-D23 in the contact map of the unbound state corresponded to the anti-parallel β-sheet structure formed at these regions. The contact probabilities in these regions were extremely high, corresponding to the high β-sheet structure of these regions in Figure 6C. A cross-diagonal represented that the β-sheet structure was also formed by residues H14-A21 and V24-V36 of the bound peptide, whereas the probabilities were not as high as the β-sheet structure in the unbound peptide. Besides, there was a band indicating the contacts formed between residues A2-H6 and residues G29-G33, which corresponded to the β-sheet interactions at these two fragments. The highest 10 contact fractions for the residue pairs (at least three residues apart) were used to characterize long-range interactions (Supplementary Figure S6).

To explore the peptide-water interactions, the average number of H-bonds formed between peptides and water molecules was counted (Figure 8B). Surprisingly, for most residues in the bound state, the average number of H-bonds formed with water was higher than residues in the unbound state except for a few residues at CT. This phenomenon occurred due to the fact that bound peptides were more exposed to water than unbound peptides and water molecules could drive the binding of peptide to the bilayer. Another reasonable explanation is that Aβ40 in extended solvent-accessible structure preferred binding with the bilayer than in the globular structure. These results suggested that water molecules played an important role in the binding of Aβ40 to membranes.