Molecular dynamics simulation on the effect of transition metal binding to the N-terminal fragment of amyloid-β

Abstract We report molecular dynamics simulations of three possible adducts of Fe(II) to the N-terminal 1–16 fragments of the amyloid-β peptide, along with analogous simulations of Cu(II) and Zn(II) adducts. We find that multiple simulations from different starting points reach pseudo-equilibration within 100–300 ns, leading to over 900 ns of equilibrated trajectory data for each system. The specifics of the coordination modes for Fe(II) have only a weak effect on peptide secondary and tertiary structures, and we therefore compare one of these with analogous models of Cu(II) and Zn(II) complexes. All share broadly similar structural features, with mixture of coil, turn and bend in the N-terminal region and helical structure for residues 11–16. Within this overall pattern, subtle effects due to changes in metal are evident: Fe(II) complexes are more compact and are more likely to occupy bridge and ribbon regions of Ramachandran maps, while Cu(II) coordination leads to greater occupancy of the poly-proline region. Analysis of representative clusters in terms of molecular mechanics energy and atoms-in-molecules properties indicates similarity of four-coordinate Cu and Zn complexes, compared to five-coordinate Fe complex that exhibits lower stability and weaker metal–ligand bonding. Communicated by Ramaswamy H. Sarma

Metal ion coordination has distinct effects on the structure and chemistry of Ab, including its propensity to aggregate. The anchoring effect of metal coordination provides structural restraint and potentially a seed for aggregation. Indeed, many studies report that metal ions act as strong inducers of Ab aggregation (Bush et al., 1994;Esler et al., 1996;Huang et al., 2004) though the results are diverse (Miller, Ma, & Nussinov, 2012), including a range of fibrillar and non-fibrillar, amorphous aggregates. Furthermore, these results are strongly dependent on reaction conditions. The redox activity of Cu and Fe may also lead to the formation of reactive oxygen species that place diseased nerve cells and synapses under oxidative stress (Greenough et al., 2013).

Computational methods
Ab1-16 was constructed in an extended conformation in MOE (Molecular Operating Environment (MOE) (2013.08), 2013) with the appropriate residue protonation states for physiological pH. Metals were coordinated to the peptide as shown in Figure 1 and subjected to brief LowMode conformational searches (Labute, 2010) to obtain reasonable starting structures. MD simulations were performed using the AMBER16 package (Case et al., 2016). The AMBER ff14SB forcefield parameter set (Maier et al., 2015) was used to model all standard amino acid residues, while metal parameters were obtained using the MCPB.py program (Li & Merz, 2016). Here, parameters are obtained from B3LYP/6-31G(d), and RESP charges for the metal-coordinating regions were obtained at the same level of theory using Gaussian09 (Frisch et al., 2009). The geometry of each system was optimised using 1000 steps of steepest descent and 1000 steps of conjugate gradient methods. MD simulations were carried out in the NVT ensemble, using a Langevin thermostat to control the temperature at 310 K. The generalised Born solvation model (Clark Still, Tempczyk, Hawley, & Hendrickson, 1990;Constancien & Contreras, 1984;Schaefer & Karplus, 1996) was used to solvate all systems, with electrostatic interactions neglected beyond a cut-off of 12 Å: this approach has been shown to enhance conformational sampling of flexible systems (Anandakrishnan, Drozdetski, Walker, & Onufriev, 2015). During all simulations, the SHAKE algorithm (Ryckaert, Ciccotti, & Berendsen, 1977) was used to constrain bonds to hydrogen. Simulations were performed using a 2-fs integration time step. Analysis of the trajectories was performed using CPPTRAJ v16.16 (Roe & Cheatham, 2013) and VMD 1.9.3 (Humphrey, Dalke, & Schulten, 1996).
Secondary structure assignment was made using the DSSP algorithm (Kabsch & Sander, 1983) via CPPTRAJ. Salt bridges were defined as any contact distance of less than 3.2 Å between oxygen and nitrogen atoms in charged residues. Ramachandran maps and hydrogen-bonding plots were made using MDplot (Margreitter & Oostenbrink, 2017). Atoms-in-molecules (AIM) analysis used the AIMAll package (Keith, 2017).

Effect of Fe(II) coordination mode
Three separate 500 ns MD simulations of each Fe-Ab complex (denoted Fe_1, Fe_2 and Fe_3) were carried out, starting from the same minimised structure but with different initial velocities, randomly sampled from the Maxwell-Boltzmann distribution at 310 K. Root-mean-square displacement (RMSD) relative to starting structure was used as our primary measure of equilibration (Huy et al., 2016;. Plots of backbone RMSD for each run are shown in Figure 2, which show that simulations reach stable values after 100-300 ns; equilibration points for each simulation are included in Table S1, supplementary material. All analyses reported are taken from data extracted from frames after each equilibration point, averaged over three separate runs, leading to 900, 950 and 950 ns of simulation data for Fe_1, Fe_2 and Fe_3, respectively. 'Spikes' in RMSD are seen in some plots, for instance around 350 ns in Fe_1C: closer inspection indicates that these correspond to relatively minor changes in backbone conformation for a few residues ( Figure  S1, supplementary material), without affecting metal coordination or the overall structure of the peptide.
Averages and standard deviations of RMSD, relative to a common starting point, over these trajectories (Table S2, (Table S2, supplementary material). Averages for the three coordination modes converge to similar values of between 7.2 and 7.6 Å, with the largest value observed for the six-coordinate Fe_3. Root-mean-square fluctuations (RMSF) for individual amino acids exhibit similar behaviour, as shown in Figure 3 (data in Table S3, supplementary material). Metal-binding residues exhibit the lowest RMSF valuesthe most mobile being Glu3 in Fe_3, with an RMSF of 5.69 Åwhile larger values of up to 10 Å are found for aromatic residues Phe4 and Tyr10, and also for the C-terminal Lys16 residue.
Metal coordination can exert profound effects on the structure of Ab, which has an important role in its neurobiology. Figure 4 shows secondary structure by residue for each Fe(II)-Ab system, as well as Ramachandran maps for the equilibrated trajectories; Table 1 shows total peptide secondary structure percentage. Plots are very similar, exhibiting large amounts of bend structure in residues 3-11 and a mix  of aand 3 10 -helix for residues 12-16. The six-coordinate structure favours turn over bend in the N-terminus, but is otherwise similar. Very small amounts (<0.25%) of b-strand structure are found to be present in all three systems.
Ramachandran maps reveal that all three complexes have the greatest concentration of conformations in the region expected for aand 3 10 -helix, along with smaller but significant populations of b-strand and left-handed helix, the former being more prevalent for the two five-coordinate complexes (Fe_1 and Fe_2). Residue contact maps are also broadly similar, exhibiting relatively short contacts between coordinated residues and longer distances between central and C-terminal ones. Clustering analysis of all equilibrated trajectories on the basis of peptide backbone dihedrals was also performed, leading to the data in Figure 5, including the minimised structures of the most representative clusters. For both Fe_1 and Fe_3, a single cluster accounts for over 90% of conformations, while for Fe_2, the dominant cluster accounts for 57% of simulation. All three show similar properties, with coil/bend/turn throughout the N-terminal region and a small amount of helix towards the C-terminus. These aspects of the secondary structure are summarised in Table 1, which reports that each simulation spends between 15 and 20% of the time in helix-like conformations, less than 0.25% in strand-like, and over 80% in turn/bend/coil conformations.   Overall, we surmise that the specifics of coordinating residues do not strongly affect the conformations visited by or the flexibility of the peptide chain.

Comparison of Cu, Fe and Zn binding
We therefore decided to compare data from the Fe_1 coordination mode to analogous structures bound to Cu(II) and Zn(II), since these share many of the same binding residues. All three metals are bound through Asp1, His6 and His14, with Glu11 as an additional binding residue for Zn(II) only. Figure 6 shows backbone RMSD values for each simulation, while Table S4, supplementary material, reports RMSD and R g values for postequilibration frames for the three systems studied. In each case, RMSD values are relatively small (2-3.5 Å) with small standard deviations (0.2-0.35 Å). Maximum values are very slightly larger for Zn than for Cu and Fe, but differences are rather small. A similar pattern is seen for R g data, in which Fe is notably more compact than the others and Zn slightly larger than Cu. These data suggest that the five-coordinate nature of the Fe complex restricts the peptide to a smaller, more compact set of structures. As expected, given the harmonic bonding model used for metal coordination, bond lengths involving metal ions oscillate around equilibrium values of around 2.0 Å with standard deviations of ca. 0.05 Å. Bond angles exhibit a similar pattern, exhibiting standard deviations of 2-3 , and reflect the approximately square-pyramidal, distorted square-planar and tetrahedral coordination geometries of Fe, Cu and Zn, respectively. Full data are reported in Table S5, supplementary material.
Further information on the impact of metal coordination on peptide conformational freedom can be drawn from RMSF data, as reported in Table S6, supplementary material, and Figure 7. The overall pattern for all three metals is similar: metal-bound residues are the least mobile: His6 and His14 in particular exhibit lower RMSF values. Reflecting the values above, Fe_1 has the lowest average RMSF and the lowest value for 9 of the 16 residues, while Zn has the highest average RMSF and is highest for nine residues.
Metal binding has important effects on the structure and chemistry of Ab: Figure 8 reports a summary of the secondary structure across equilibrated trajectories, and Table 2 shows total peptide secondary structure percentage. Secondary structures show a good deal of similarity between metals: amino acids Asp1-Tyr10 exhibit mainly coil/turn/ bend structures, while Glu11-Gln15 show a much more helical character, along with turn. Within this broad picture, however, subtle details differ: Zn induces more a-helical content than the other metals, which extends further towards the C-terminus of the peptide, whereas Cu and Fe show more 3 10 helix characters over a slightly shorter range. A small amount of parallel b-strand contact between Ala2 and Gln15 is also seen for Fe_1, and an even smaller extent of anti-parallel strand for Ser8 and Gly9 in the Cu complex. Ramachandran maps reflect this pattern, but also show more detail of conformations sampled by MD. All show the greatest concentration within the right-handed helix region of the maps, along with variable amounts of b-strand and lefthanded helical conformations. Notably, Cu shows a large number of conformations in the PPII region (ca. À75, 150) (Adzhubei & Sternberg, 1993;Hollingsworth & Karplus, 2010;Woody, 2009). The PPII structure is relatively open and contains little/no internal hydrogen bonding (vide infra) and as  such is likely characterised as turn/bend/coil. Fe_1 and Zn show progressively fewer conformations in this region. In addition, Fe_1 shows significant occupation of the bridge and ribbon regions of the plot, extending between the broad helical and b-sheet regions (Hollingsworth & Karplus, 2010;Karplus, 1996) Fe_1 also has a greater number of conformations in the b-sheet region of the Ramachandran map than either Cu or Zn, though this is not reflected in secondary structure analysis, likely due to their lack of requisite hydrogen bonds (vide infra). Overall though, we stress that the Nterminus of Ab remains rather unstructured despite the metal coordination imposing structural restraints on certain residues, as exemplified by the data in Figure 8. Salt bridges are crucial in determining peptide structure and stability: Figure 9 shows the percentage occupancy of all eight possible such contacts within Ab 1-16 (numerical data in Tables S7-S9, supplementary material). In all systems, the most commonly observed interaction is between Arg5 and Asp7, despite the fact that the metal ions are coordinated through His6 between these residues. This bridge is present for almost 100% of the sampled configurations of the Fe_1 complex, and between 75-85% for Cu and Zn. Glu3-Lys16 is also very prevalent (91%) for Fe_1, but is barely observed for Cu and Zn for which Glu3 tends to make contact with Arg5 (38 and 10%, respectively). Notably, the Asp7-Lys16 interaction is observed much more frequently (14%) in the Cu simulations than in either the iron or zinc. Asp1 and Glu11 rarely form salt bridge interactions in any of the metal systems.
The effect of these salt bridges, and of the metal coordination patterns, is evident in the contact maps shown in Figure 9. Fe_1 contains rather short distances for distant amino acids that share coordination, for instance, between Asp1 and His14 and also between Glu3 and Lys16 that participate in a high-occupancy salt bridge. In contrast, Cu and Zn complexes have rather fewer close contacts, reflecting the lower coordination number and fewer salt bridges in these complexes. Particularly long distances are observed for Glu3-Val12 in Cu, and Glu3-Tyr10, His13 and Lys16 in Zn.
As well as salt bridges, hydrogen bonds play an important role in determining stability and structure: Table 3 summarises the overall H-bond patterns for the three metals. All are rather transient, with frames containing zero hydrogen bond interactions observed as well as others with as many as 15. On average, Fe_1 contains the most: contacts between backbone C ¼ O of Glu11 and H-Ne of His14 (58% of trajectory) and between backbone C ¼ O of Phe4 and H-Ne of His6 (57%) are the most prevalent. In contrast, the most occupied H-bond in the Zn complex is between Tyr10 C ¼ O and H-Ne of His14 at just 33% of frames, while for Cu, the contact between Lys16 C ¼ O and H-Ne of His6 is observed for 42% of frames. Figure 10 adds more detail to this picture and highlights significant differences between metals: both Cu and Zn contain frequent interactions between Arg5 donor and Glu3 and Asp7, whereas Fe contains only the former interaction. Both Cu and Fe display frequent His6 (donor) to Glu3 and Phe4 (acceptors), which are not present in Zn. In contrast, Zn Glu11 accepts H-bonds with four different residues with greater than 5% incidence, the most of any residue. Interestingly, Glu11 is coordinated to Zn in this case, while the uncoordinated Glu11 in Fe and Cu cases forms Hbond interactions with no more than three residues, though at greater incidence than for Zn.
Clustering of all recorded frames on the basis of backbone dihedrals was also performed, leading to the data shown in Figure 11, including ribbon diagrams of the most representative cluster following minimisation. For Cu and Zn, a single cluster accounts for over 60% of the trajectory, with the remaining frames spread across clusters encompassing no more than 17%. The Fe trajectory, however, leads to essentially a single cluster accounting for the entire set of postequilibration frames recorded at 92%, with no other cluster covering more than 5%. This reinforces the picture that the extra metal-ligand bond in the five-coordinate iron complex  limits the flexibility of the peptide more than is the case for four-coordinate metals. We also report the total MM energy after minimisation, which shows that the Cu and Zn complexes have similar energy and are markedly more stable than the Fe_1 complex. Further inspection of the breakdown of this energy (Table S10, supplementary material) indicates that bond stretch, angle bend and van der Waals terms are within 2 kcal/mol for all three systems and that the overall energy difference stems from changes in electrostatic and generalised Born solvation terms. It seems that higher coordination number of the iron complex prevents the polar side chains from interacting electrostatically with either the rest of the peptide or the implicit solvent, compared to the four-coordinate Cu and Zn complexes. Optimised geometries of representative clusters were used to carry out AIM analysis, based on electron density from B3LYP/6-31G(d) calculations. This is focussed on the metal-ligand bonds: bond critical points (BCPs) are located for all coordination bonds: properties at these points are summarised in Table 4 Figure 11. Summary of clustering analysis: minimised geometry, percentage occupancy and molecular mechanics energy (kcal/mol) after minimisation. Figure 10. Plots of hydrogen bond occupancy for Cu (left), Fe_1 (centre) and Zn (right). Plots only display data for hydrogen bonds present for greater than 5% of simulation. Black circles indicate multiple interactions between relevant residues. .034 a q is the electron density at the bond critical point; ᭢2q is the Laplacian of the electron density here; e is the ellipticity of the bond; G is the kinetic energy density at the bond critical point; V is the potential energy density here; H is the total energy density, the sum of G þ V. All properties are reported in atomic units.
bonds being notably stronger than those of other atoms, exhibiting greater electron density (q), lesser positive Laplacian (r 2 q) and more negative energy density (H). By these measures, the strongest of these bonds are between Cu and Nd of His6 and His14, while those to Asp1 are around 10% weaker. This pattern is less marked for the other metals, but in general, these data indicate that coordination through oxygen is weaker than through N. Particularly weak, and hence more labile, bonds are observed for Fe-O(His6) and Zn-O(Glu11).

Conclusions
We have used atomistic MD to examine and compare the effect of coordination of copper, iron and zinc on the structure and dynamics of the N-terminal 1-16 fragments of Ab. These ions are believed to play a role in the aggregation of Ab and in the onset of AD, such that their effect on peptide conformational behaviour may be an important factor in disease. In order to compare different metals and coordination modes on an equal footing, we use Li and Merz's MCPB.py approach, which extracts force constants and non-bonded parameters for metal ions from DFT calculations. These parameters were calculated for three different proposed coordination modes for Fe(II) and for comparable coordination of Cu(II) and Zn(II).
This method is used to obtain 1.5 ls of MD trajectories for two five-coordinate Fe(II) complexes, which differ in coordination through either His13 or His14, and a six-coordinate complex involving also the side chain carboxylate of Glu3. Each system is found to reach pseudo-equilibration after 100-300 ns, allowing us to extract over 900 ns of equilibrated data. Over the course of these trajectories, the conformations adopted by the peptide are broadly independent of the coordination mode: all exhibit turn/coil/bend structures in the N-terminus and a significant amount of helical content towards the C-terminus. Similar patterns in the flexibility of individual residues, as measured by RMSF, and inter-residue contact distances are also seen.
Comparison of one of these Fe(II) coordination modes with analogous ones to Cu(II) and Zn(II) was then made, from which we find more notable differences between metals. Equilibration again occurs within similar time scales, allowing extraction of at least 900 ns of equilibrated trajectories. The overall picture of turn/bend/coil in N-terminus and helix in C-terminus is retained, but the amounts of each vary more between metals. In particular, Cu coordination leads to greater occupancy of the PPII region of Ramachandran maps more than other metals, while Fe leads to significant occupation of bridge and ribbon regions. Salt bridges are also affected differently by different metals: Arg5-Asp7 is ubiquitous, whereas Glu3-Lys16 is common for Fe complexes, but is hardly present for Cu/Zn. These changes are evident in inter-residue contacts: the five-coordinate Fe complex is more compact than the four-coordinate Cu/Zn ones, with Zn in particular allowing greater peptide flexibility. Potential biological consequences of these changes are difficult to estimate, especially since we have only considered the N-terminal portion of the peptide. Nevertheless, we speculate that the rigidity of the Fe complex might provide a better seed for aggregation than the more flexible Cu and Zn complexes.
Clustering allows us to extract a single structure that is representative of the equilibrated trajectories: after minimisation, the energies of the Cu and Zn complexes are similar, but that of the Fe complex is significantly higher, apparently due to the strain imposed by the extra coordination bond for this metal and the lower solvation energy that results from coordination of the carboxylate of Glu3. DFT and AIM analysis of the metal-peptide bonds within these minimised clusters indicates that metal-histidine bonds are strongest, especially within the Cu complex, and that bonds to N are slightly stronger than those to O.
This work focusses on the N-terminal fragment of Ab since this is known to be the site of metal binding, while the relatively compact nature of the complexes formed aid equilibration of MD trajectories. Having established suitable parameters for simulation of different transition metals, future work will address how the same methods can be applied to the full 1-40 or 1-42 sequences and their aggregation.