Several fibrillary structures of the β-amyloid fibril are available in the Protein Data Bank. The most complete PDB structure, with code 2NAO (Wälti et al., 2016), corresponds to the sequence 1–42 of the amyloid for each strand, though it retrieves poor punctuation in all the PDB parameters (clashes score 22, Ramachandran 12.4% and Sidechain 26.7%). On the other hand, the PDB structure with code 2MXU (Xiao et al., 2015) presents far better punctuation (clashes score 0, Ramachandran 3.8% and sidechain 6.22%) at the expense of residues 1 to 11 being missing, an observation consistent with the high flexibility of the N_Ter region of the peptide. Therefore, for this study, we selected the latest and elongated each strand with the 11 missing residues using homology modeling through the UCSF Chimera Software Modeler tool (Pettersen et al., 2004) to obtain the N-terminal missing region. Chimera minimization was performed to reduce the clashes, with the Amber ff14SB force field (Maier et al., 2015) and 100 steps of steepest descent.

Once the complete fibrillary structure was obtained, several preliminary molecular dynamics calculations were performed. We considered the systems with different numbers of strands, namely 4, 6, 8, and 10 aiming at identifying the minimal number of strands necessary for the fiber to remain stable in the fibrillar form during the trajectory. This is also particularly important to avoid possible drawbacks inherent to force field propensity to stabilize certain secondary structures against others (Case et al., 2005). With the Amber force field, we observed that at least 6 strands are necessary to obtain a fibril, though at least 10 strands should be considered for a stable core. In this study, only the fibrillary structure of 10 strands is discussed.

Putative metal-binding areas were detected by applying the BioMetAll software (Sánchez-Aparicio et al., 2021) to the Aβ_11–42 crystallographic structure available in the PDB 2MXU, which was configured to include the backbone atoms and a minimum of three residues as possible coordinating groups. BioMetAll only considers the backbone and the residues that are within a suitable distance to coordinate the metal, but their orientation may not be optimized for coordination. Therefore, protein-ligand dockings calculations were performed with the program GOLD (Jones et al., 1997) using the GoldScore scoring (Jones et al., 1995) function, whose parameter file was modified to include atom types for metal ions and their possible coordinating amino acids (Sciortino et al., 2018). Genetic algorithm (GA) parameters were set to 50 GA runs and a minimum of 100,000 operations each. The remaining parameters, including pressure, number of islands, crossovers, or mutations, were set to default. Finally, docking solutions were analyzed through GaudiView (Rodríguez-Guerra Pedregal et al., 2017) an in-house developed GUI tool built as an extension of UCSF Chimera.

Metal parameters were obtained using the MCPB.py package (Li and Merz, 2016) from quantum mechanical calculations with DFT (B3LYP) (Lee et al., 1988; Becke, 1998) and adding Grimme’s D3 correction for dispersion (Grimme et al., 2010). For aluminum complexes, the 6-31 + G (d,p) basis set was used for all atoms (Rassolov et al., 2001). For Cu(II) complexes, we used the 6-31 + G (d,p) basis set for C, H, N, and O atoms; and the SDD pseudo-potential and the corresponding basis-set supplemented with f-polarization function for Cu (Ehlers et al., 1993). Solvent-polarizable dielectric continuum model (SMD) was considered to account for the solvent effects in water (Marenich et al., 2009). Force constants and equilibrium parameters for metal-coordinating atoms were obtained through the Seminarios’s method (Seminario, 1996) while point charges were derived using the RESP (restrained electrostatic potential) model (Bayly et al., 1993).

Molecular dynamics simulations were then performed with AMBER18 (Case et al., 2018) using the AMBER ff14SB (Maier et al., 2015) force field in the NPT ensemble, with a 1 fs integration time step. The initial models were embedded within a cubic box of pre-equilibrated TIP3P (Jorgensen et al., 1998) water molecules and Na⁺ ions were included to balance the total charge depending on the system: 30 for the metal-free system, 20 for the copper-bound system, and 15 for the aluminum-bound system. Constant temperature and pressure were set by coupling the system to a Monte Carlo barostat (Duane et al., 1987) at 1.01325 bar and a Langevin thermostat at 300 K (Loncharich et al., 1992). The SHAKE algorithm (Ryckaert et al., 1977) was used to constrain the bonds involving hydrogen atoms. Classical MD simulations were performed for 10 ns for all the systems, just to equilibrate the system before running Gaussian accelerated molecular dynamics (GaMDs). Moreover, MD simulations have also been performed on the lowest energy wells obtained in the GaMD simulations.

Enhanced conformational sampling was performed with the Gaussian accelerated molecular dynamics (GaMDs) (Miao et al., 2015). For these simulations, the AMBER ff14SB force field in the NVT ensemble was used, constraining the bonds involving H atoms with SHAKE, with an integration time step of 2 fs. A boost on both dihedral and total potential energy (igamd = 3) was applied. In total, three GaMD replicas of 1 μs were produced for each system.

The energy of the systems was calculated along the trajectory with the Generalized Born Implicit model (Onufriev et al., 2000, 2004), for 1,000 frames extracted from each 1 μs GaMD trajectory, stripping water molecules, and performing a short minimization (maxcyc = 500 steps) previously to the energy extraction. The radius of gyration (RoG) was computed over the same 1,000 frames, so a graphic with the different minima was obtained. For each system, the well with the minimum energy was selected for an MD, to test its stability, and for the analysis.

Principal component analysis (PCA) (Wold et al., 1987), RMSD all-to-all, and clustering were performed over both the GaMD and MD simulations to test whether the exploration was exhaustive enough. For the MD simulations performed on the lowest energy wells from the GaMDs, the RMSD with respect to the initial structure was also computed. The last analysis performed consisted of studying the PCA movements of the fibers along the trajectories with the VMD tool NMWiz and counterchecked by normal mode analysis performed with WEBnma software (Humphrey et al., 1996; Bakan et al., 2011) to study the global low-energy internal motions of the systems.

The initial structure of the metal-free fibril has been obtained from the solid-state NMR (ss-NMR) structure of Aβ_11–42 (pdb code 2MXU). Since the N_Ter end of the system is missing (residues 1 to 10), the structure was completed by homology modeling using modeller (Martí-Renom et al., 2000; Figure 2A). Then, long GaMD simulations were performed using a biased potential and the results were analyzed using several statistical tools.

The total 3 μs GaMD simulation spread over three replicas of 1 μs. The three replicas exhibit similar behaviors in terms of global breathing motions and exhaustive exploration of the conformational space. Accordingly, PCA and energetic profiles of the replica with the lowest energetic conformation identified are presented in the main text (Figures 2B, C) while the other ones are reported in ESI (Supplementary Figure 1). Principal component analysis shows a displacement of the system at the end of the first 100 ns, then, it explores a limited region of the space until approximately 600 ns, and finally, it explores back the previous region, so that the simulation is converged (Figure 2B). The potential energy profile clearly indicates that during the trajectory, a series of low-energy wells are explored with a radius of gyration slowly raising–meaning a geometric expansion of the system (Figure 2C). This is in line with the results reported by Miller et al. (2012), who suggest that fibers do not show a unique lower energy state but different conformers with barriers between them, which are expected to be more or less favorable upon metal binding.

The trajectories present a series of interesting features. The first one is that only the six central strands maintain the fibrillary structure while the two at each extreme tend to unfold. This is related to the fact that these extremes are more exposed to the solvent and naturally start to adopt a more molten globule geometry with some helical contents, a situation reminiscent of the simulations of monomeric peptides. The second is that the terminal regions of each strand represented by the amino acids from 1 to 10 are overall disorganized. This is consistent with the absence of atomic density in the pdb structure and suggests that the overall geometry of the amyloid is barely affected by the N_Ter region adopting quite random conformations. Based on this, the N_Ter region and the four strands at the extremes are discarded for the rest of the analyses, which focuses on the core formed by the six central strands and from residues 11 to 42.

The major motion observed is a fan-like movement of the Aβ_11–20 β-sheet, which appears to be the natural breathing movement of the system (Figure 2D) as counterchecked by normal mode analysis performed using WEBnma software (Tiwari et al., 2014). This fan-like motion consists of a hinge that can be defined as the angle between the principal axes of the sets of atoms of residues 11 to 15 and 39 to 41 (Figure 2E). The angle obtained between these two principal axes for a representative structure of the most populated cluster reaches 48.5° against the 7.8° of the ss-NMR structure. For the lowest energy structure, the same analysis shows a hinge of 46.7° (Figure 2F), which still represents the main distortion with respect to the ss-NMR structure. These simulations and analyses show that the metal-free fiber mostly retains the experimental geometry although it can display the dynamic motions not observed in the experimental structures.

The binding of metal ions to amyloids has been the source of an extensive number of studies (Faller, 2009; Mutter et al., 2018; Pham et al., 2018; Roldán-Martín et al., 2021; Miller, 2022). Today, a relatively clear picture of metal binding has emerged in the monomeric systems, either theoretically for Al(III) (Mujika et al., 2017) or experimentally for Cu(II) (Drew and Barnham, 2011), which has narrowed the location and type of amino acids that bind the metal ions with only few candidates, all part of the N_Ter region of the peptide. For aggregates, though, there is still a substantial level of uncertainty.

The previous part of our study sustains an intrinsically disordered protein (IDP) behavior of the N_Ter region of the aggregate. This flexibility indicates that the binding of the metal ions to this region would be similar to what occurs in monomeric species. It also shows that the core of the aggregate is quite rigid, and one would naturally wonder if this rigidity induces some pre-organizations for metal binding.

To identify putative metal-binding sites in the aggregates, we performed BioMetAll calculations, a software recently released by our group for the prediction of metal-binding sites in biomolecules. Based on the backbone pre-organization hypothesis, BioMetAll has shown excellent success rates in identifying metal-binding locations for a wide sample of proteins and ions (Sánchez-Aparicio et al., 2021). To search for site-specific binding areas for copper and aluminum, the search was first performed looking for a minimum of three coordinating amino acids and any amino acid known to participate in the coordination sphere of transition metals. With this screening, BioMetAll finds the geometries compatible with the simultaneous coordination of two His residues in the N_Ter region of the Aβ fibril as identified in the monomers (Alí-Torres et al., 2011, 2014) from the experiments and QM calculations. However, this only occurs in a few strands due to the intrinsic flexibility of that region. Since other metal-binding sites are more frequently found in the core region of the fibril, indicating that in this system, the situation is different from what happens in monomers, we discarded N_Ter monomer-like binding sites and focused on the core region, favoring Asp and Glu for Al(III) (Mujika et al., 2017) and His and Gln for Cu(II). For the latter, it is important to notice that Gln residues have been included in the search since, although rare, the coordination of Gln to Cu(II) is observed in a series of X-ray structures from the protein data bank, it has also been proposed as a possible coordinate residue (Parthasarathy et al., 2011), and Gln15 shows the proper location for metal coordination in the selected structures.

In total, two areas were found for metal binding in the core of the fibril: one consistent with copper and the other with aluminum. The former, named Cu_F, stands at the linker between the N_Ter loop and the β-strand regions (from residue 12–15) and is composed of histidine His13 of two adjacent strands and two glutamine Gln15 of the same strands (Figure 3A). The latter, named Al_F, is located at the beginning of the coil with residues Glu22 and Asp23 of adjacent strands (Figure 3B). This latter site agrees with the work reported by Yang et al. (2022) which shows that this area is likely to be bound to a metal ion. With these two sites identified, we pursued generating realistic three-dimensional models of the metal-aggregate systems by carrying out protein-ligand dockings. Those were performed with GOLD, assuming a square planar environment for Cu(II) and an octahedral one for Al(III) (Figures 3C, D).

The results indicate that Cu(II) can achieve a square planar coordination with two His and two Gln from adjacent strands. However, for Al(III), docking results only find an incomplete octahedral coordination in which two Glu and two Asp of adjacent strands are involved, thereby leaving vacant sites that would probably be occupied by the solvent molecules. Each model was then submitted to triplicated Gaussian accelerated molecular dynamic simulations to evaluate its stability and dynamics properties as well as to compare with the metal-free system.

Data of the triplicated simulations on Cu_F are given in Figure 4 and Supplementary Figure 2 of ESI. In all cases, the trajectories tend to converge around 400 ns (Figure 4B and Supplementary Figures 4D–F) and show the presence of a series of minima (Figure 4C). Such minima agree with the fact that amyloids are highly polymorphic, metal binding shifting the population toward a certain conformation (Miller et al., 2010). Overall, the 3 μs GaMD simulations depict the same fan-like motion observed in the metal-free fiber (Figure 4D). Strikingly though, the geometry variability of the N_Ter loop (residues 1–11) is significantly more pronounced in Cu_F; i.e., each metal-bounded pairs of strands have their motions uncoupled from adjacent pairs. Consequently, the β-sheet H-bonds between these adjacent pairs are significantly weakened, which favors the change in their secondary structure increasing the α-helix content in the N_Ter ends, particularly on the first strands. Such α-helix regions are reminiscent of the monomeric forms, a behavior more pronounced than in the metal-free system (Roldán-Martín et al., 2021). However, the hinge movement previously observed in the metal-free fiber appears to be reduced, the angle between the two principal axes (considering the six central strands as done previously for the metal-free fiber) in a representative structure of the most populated cluster being 6.1° (Figure 4E). Note that in this case, the motion is not globally spread all other the fibril since there is a disruption of the fibrillar form that arises from a higher displacement of the first strands with respect to the remaining ones, probably due to the repulsion between Cu(II) ions (Figure 4F). Such observations may indicate that Cu(II) binding tends to disrupt the fibrillar form to lower-aggregation structures or amorphous aggregates without fibrillary structure.

For Al(III), the preferred binding site in the fiber involves the residues Glu22 and Asp23 of adjacent strands (Figure 5A). In this last case, the GaMD explores a series of minima until finding a well it cannot escape during the last 400 ns (Figure 5B and Supplementary Figure 3). Despite the energy profile highlighting a remarkable minimum (Figure 5C), the whole trajectory displays a very stable fold with minor motions (Figure 5D). In particular, the amplitude of the fan-like movement is far more limited in Al_F than in the other two systems. For example, one of the most representative structures presents an inter-axial angle of only 2.7°, similar–even lower–than that in the ss-NMR structure (7.8°) (Figures 5E, F). Therefore, GaMD simulations show that the aluminum-bound aggregate displays a more compact structure than both free-metal and copper-bound ones and remains closer to the static view provided by experimental ss-NMR geometry (Figure 5F).

To further analyze how the metal ions impact the conformation landscape and stability of the aggregates, classical MDs of 100 ns were carried out starting from the structure of the lowest energy well obtained from the GaMD simulations. To test the exhaustiveness of the exploration and their stability, RMSD and PCA were performed on the trajectories (Supplementary Figure 4). Emphasis is given to how metal binding alters the fiber’s low collective modes and breathing motions. For that, five distances were evaluated using the carbon α of each residue and considering only the six central strands of the fibril. To ease the labeling and the reading, we use the X-ResNum identification format where X is the number of the strand the residue belongs to, Res is the type of the residue, and Num is its number in the sequence.

With these descriptors, it is possible to see that the three complexes behave similarly both in the horizontal and in the vertical axes (measures 1 and 2). For measure 1 (Figure 6E), the metal-free fiber shows a distance of 24.64 (±5.6) Å, while Al(III) and Cu(II) exhibit values around 25 Å, which are very similar to the 24.6 Å of the ss-NMR structure. Distances of measure 2 (Figure 6F) are also similar, with the values of 19.9 (±0.3) Å for the metal-free fiber, 18.7 (±0.3) Å and 18.8 (±0.2) Å for the Cu and Al bound complexes and 18.3 Å for the ss-NMR structure, respectively. The other distances, though, demonstrate the substantial differences between the metal-free fiber and the metal-bound complexes.

Measure 3 relates to the disruption of the fibrillary structure in the N_Ter region. It is observed that Cu(II) coordination produces an increase in this distance, from 24.3 Å in the ss-NMR to 26 (±1.1) Å. In contrast, both the metal-free fibril and Al(III)-bound structure show shorter distances, of 19.93 (±0.3) Å and 21.5 (±0.7) Å, respectively, in comparison with the ss-NMR structure (Figure 6G). Measure 4, on its side, is remarkably increased in the metal-free fibril [32.1 (± 0.4) Å] with respect to the ss-NMR (29.3 Å), due to the hinge fan-like movement. Cu(II) system is maintained in a similar value to the ss-NMR structure, with a distance of 29.6 (±0.5) Å. Note, however, that this distance is an average of six values and the displacement of the terminal strands is faded by the shorter distance retrieved by the other ones. Indeed, the two most external strands included in the average exhibit values of 35 and 42 Å, respectively, whereas the four inner ones show values that range between 20 and 25 Å. Finally, the Al(III) system is even more compact than the ss-NMR structure, with an average distance of 23.1 (±0.4) Å, supporting the fact that aluminum binding increases the stability of the fibrillary conformation (Figure 6H).

The distances between each pair of strands (measure 5) were also measured (Figure 7). Calculations confirm the trends previously observed. The metal-free fibril has approximately the same distance between each pair of strands. Such system has a higher flexibility on the vertical axis, as measures 2 and 4 also support, while the horizontal axis is more compact and shows a collective movement of all the strands, which corresponds to the hinge fan-like movement. However, upon metal coordination, the distance between the strands differs.

First, Cu(II) ion, which binds to the N_Ter region, increases the distance between each pair of strands, even between those bound to the same metal ion. The most remarkable is the increase in the 4-His14–5-His14 distance, which indicates a disruption in the collective movement of the strands (Figure 7). This only happens in the Cu-bound system. For Al(III), the distance between the pairs of strands is also increased, especially in strands 3 to 6, but the collective movement is preserved. This points to the fact that, although Al(III) binding increases the compactness and stability of the core region, as measures 1 to 4 support, the (N_Ter) region is still highly flexible.

The mentioned changes in metal coordination are also clearly observed if the systems are overlapped against the ss-NMR structure (Figure 8), which demonstrates that the metal-free fiber has a collective hinge fan-like movement, increasing the flexibility of the system but without its disruption. In contrast, Cu(II) leads to the partial dismantling of the fibrillary conformation in the (N_Ter) region by breaking the collective movement of the fiber. Finally, Al(III)-bound system is more compact in the core region, though the (N_Ter) region is still highly flexible.