Benchmarking of copper(II) LFMM parameters for studying amyloid-β peptides

Ligand field molecular mechanics (LFMM) parameters have been benchmarked for copper (II) bound to the amyloid-β1–16 peptide fragment. Several density functional theory (DFT) optimised small test models, representative of different possible copper coordination modes, have been used to test the accuracy of the LFMM copper bond lengths and angles, resulting in errors typically less than 0.1 Å and 5°. Ligand field molecular dynamics (LFMD) simulations have been carried out on the copper bound amyloid-β1–16 peptide and snapshots extracted from the subsequent trajectory. Snapshots have been optimised using DFT and the semi-empirical PM7 method resulting in good agreement against the LFMM calculated geometry. Analysis of substructures within snapshots shows that the larger contribution of geometrical difference, as measured by RMSD, lies within the peptide backbone, arising from differences in DFT and AMBER, and the copper coordination sphere is reproduced well by LFMM. PM7 performs excellently against LFMM with an average RMSD of 0.2 Å over 21 tested snapshots. Further analysis of the LFMD trajectory shows that copper bond lengths and angles have only small deviations from average values, with the exception of a carbonyl moiety from the N-terminus, which can act as a weakly bound fifth ligand.


Introduction
Alzheimer's disease is one of the major health challenges facing modern society. Its etiology is complex but the importance of amyloid-β (Aβ) peptide in this has long been reported, (Masters et al., 1985) with aggregation of Aβ fibrils and plaques observed in the brains of Alzheimer's patients. The N-terminus of Aβ contains several residues that can act as ligands towards transition metals, (Bolognin et al., 2011;Kowalik-Jankowska, Ruta-Dolejsz, Wisniewska, & Lankiewicz, 2001;Zirah et al., 2006) such that the structural and/or chemical effects of metal coordination have been proposed to play a role in the onset of Alzheimer's (Bush, 2003;Spinello, Bonsignore, Barone, Keppler, & Terenzi, 2016). Naturally occurring metals in the brain include copper, iron and zinc, the redox activity of the former two indicating a possible mechanism for damage to brain cells (Huang et al., 1999;Jomova, Vondrakova, Lawson, & Valko, 2010). Moreover, non-natural metals such as platinum show promise as potential anti-Alzheimer's agents through disruption of native metal coordination Turner, Platts, & Deeth, 2016).
A large body of experimental work has shed light on copper-Aβ coordination, (Kepp, 2012;Warmlander et al., 2013) supplemented by extensive simulation work (Ali-Torres, Mirats, Marechal, Rodriguez-Santiago, & Sodupe, 2015). Classical molecular dynamics (MD) simulations have been a powerful computational tool in elucidating structural information of copper-Aβ complexes. Raffa and Rauk (2007) used DFT benchmarks to include copper parameters in the Gromacs force field for simulations of Aβ 1-42 monomers to study how binding perturbed the peptide secondary structure. Similar approaches were carried out by Dong et al. (2016) and Kozman and Tvaroska (2015) focusing on aggregation inhibition by small molecule drugs and curcumin, respectively. Huy, Vuong, La Penna, Faller, and Li (2016) recently used microsecond MD simulations to model the structure and dynamics of Cu-Aβ 1-42 dimer interactions. Despite this extensive work there is still contention over which residues are actively involved in binding copper at physiological conditions. Reports acknowledge the importance of histidine as nitrogen ligands but several potential oxygen ligands have been reported, including those within alanine, aspartic acid, glutamic acid and tyrosine residues (Dorlet, Gambarelli, Faller, & Hureau, 2009;Drew, Masters, & Barnham, 2009;Minicozzi et al., 2008;Streltsov et al., 2008).
Extensive QM/MM studies combined with homology modelling were carried out by Ali-Torres et al. (2011) to probe the Cu(II)-Aβ 1-16 coordination sphere. The 3N1O model was considered with His6, His13 and His14 as nitrogen ligands with both δ and ε binding modes, with a range of potential residues for the oxygen ligand. A comparison of implicitly solvated free energy calculations revealed the most stable complex to contain the backbone carbonyl oxygen of the ALA2 residue, known as component II (Drew, Noble, Masters, Hanson, & Barnham, 2008). This result was also in agreement with the experimental EPR results of Barnham and co-workers (Drew et al., 2009). These calculations confirmed the preferential histidine binding modes to be ε, δ and ε for HIS6, HIS13 and HIS14, respectively. In light of the above reports the [O c A2 ,N ε H6 ,N δ H13 ,N ε H14 ] model of Cu (II)-Aβ 1-16 was chosen as the focal point of this study, as shown in Figure 1.
An alternative to use of ad hoc molecular mechanics parameters for Cu is ligand field molecular mechanics (LFMM), which captures key d-orbital effects in a small number of transferable parameters. This approach has previously been used successfully to study copper-containing transition metal complexes in a variety of systems, including type 1 copper proteins, (Deeth, 2007) type 3 copper proteins, (Diedrich & Deeth, 2008) and bispidine complexes (Bentz et al., 2008). In this work, we test the suitability of the LFMM approach for describing Cu(II)-Aβ 1-16 coordination, then employ this approach to carry out molecular dynamics simulations of this model system.
The Cu(II)-Aβ 1-16 system was constructed in MOE. The N terminus at ASP1 was left uncapped in the zwitterionic form found at physiological pH, whereas the C terminus at LYS16 was capped with an amide group. The three histidine residues were built in the neutral form and all other residues were left as standard charges at physiological pH, leading to an overall charge of +1. AMBER parm94 parameters were assigned for the entire peptide using MOE's dictionary lookup, then charges modified to include the presence of copper using the scheme of Comba and Remenyi (2002). Altered charges are reported in Supporting Information (Table S1). LFMD simulations on Cu(II)-Aβ 1-16 were carried out with the LFMM/AMBER parm94 forcefield at 310 K, with reaction field electrostatics with a van der Waals cutoff of 10 Å and long range cutoff of 21 Å. A 10 ns simulation was run in the NVT ensemble, with a 1 fs integration timestep used throughout.
Benchmarking calculations were carried out with a variety of hybrid DFT methods, namely M06-2X, BHandH and B3LYP both with and without Grimme's empirical dispersion correction, B3LYP(-D2) (Grimme, 2006). The 6-31G* basis set, combined with LANL2DZ and associated pseudopotentials for copper was used throughout. This level of theory was deemed suitable by Sousa et al. (2013) for studying systems containing copper. Aqueous solvation was accounted for in all DFT calculations with the use of the CPCM polarisable conductor model. Semi-empirical calculations were carried out using the PM7 method, (Stewart, 2013) with the COSMO solvation model.

Results and discussion
The similarity of the reported distorted tetrahedral coordination of Cu(II) by Aβ to that seen in Type I copper proteins led us to use those reported previously for the latter as an initial set of LFMM parameters  (Deeth, 2007). Atomic charges for Cu-coordinated residues were calculated using Comba's scheme (Comba & Remenyi, 2002). Parameters were then refined against DFT structures of model complexes (vide infra), then tested for Cu(II)-Aβ 1-16 and smaller models representative of the metal binding within this system. LFMM parameters were initially tested against DFT results on six simple four-coordinate copper complexes. These models are representative of potential binding sites of copper to Aβ in the manner described above. Imidazole (imid) and formamide (form) are used to represent histidine and carbonyl backbone coordination, respectively. The models include all the potential combinations of these two ligands. Figure 2 compares these complexes optimised with LFMM, in yellow, and B3LYP-D2, in blue. Parameters (r 0 and α for Morse potential plus ligand-ligand repulsion terms) were adjusted in order to minimise the root mean square deviation (RMSD) of atomic positions from DFT structures. In general, parameters were very close to those previously reported for Type I copper proteins. New LFMM parameters, in the MOE format, are reported in supporting information.
All six models in Figure 2 show good agreement between LFMM and B3LYP-D2 optimised structures, although minor deviations relating to the orientation of imidazole ligands are present in B and C. LFMM bond lengths are typically within 0.1 Å of B3LYP-D2 values, similar to previous performance reported for LFMM (Deeth, 2007). Angles around the Cu centre are generally within 1°on average for pseudo-cis ligands, and within 5°for pseudo-trans ones, with the largest difference of 10°for model B. A very similar pattern of agreement between DFT and LFMM is found using different functionals: all data are reported in tables S2 and S3 in the supporting information.
A larger model system representative of the copper binding site in Aβ, with one alanine and three histidine ligands, was also used to test LFMM parameters. A comparison of LFMM (yellow) and B3LYP-D2 (blue) optimised geometries is shown in Figure 3, which shows good agreement between the LFMM and B3LYP-D2 optimised structures. Differences in bond lengths of 0.08 and 0.01 Å are found for average Cu-N and Cu-O, respectively, and 5°and 1°for O-Cu-N and N-Cu-N angles, respectively. These small differences show that the LFMM copper parameters model the copper coordination geometry to a high degree of accuracy against DFT results.
The generally excellent agreement for truncated model systems optimised by LFMM and DFT led us to use the former in MD simulations on the Cu(II)-Aβ 1-16 system. The good performance of PM7 might also offer a means to carry out dynamical simulation, but the extra speed of the LFMM approach makes it especially attractive for this purpose. The 16 residue peptide with copper was constructed in MOE and optimised with LFMM. To  find a suitable starting point for LFMD simulations, conformational searching was carried out with low mode molecular dynamics, (Turner et al., 2016) with the lowest energy structure from the low mode molecular dynamics used as the starting point. A plot of RMSD from the starting point for a 10 ns trajectory is shown in Figure 4(a). This shows large changes in the initial 2.5 ns of simulation, before reaching a plateau at approximately 1.9 Å for the bulk of the remaining 7.5 ns.
However, there is a sharp drop in RMSD to around 1.6 Å after 4 ns, before a return to oscillation around 1.9 Å over the rest of the trajectory. Even here though, occasional structures are observed to have small RMSD values: the possible structural origin of this is discussed in more detail below. Further evidence for equilibration is provided by the calculated radius of gyration, R g , which falls from its initial value of around 7 Å to reach a stable oscillation about 6.6 Å within 2.5 ns (Figure 4(b)). To examine the validity of structures generated through LFMD simulations, three snapshots were extracted from the trajectory at 5, 7.5 and 10 ns, and optimised with LFMM, PM7 and DFT methods. As with the smaller models above, Cu-N and Cu-O bond lengths from LFMM and B3LYP-D2 agree within 0.1 Å, and most angles around Cu to within 1°. An exception is the pseudo-trans O(Ala2)-Cu-N(His14), which is predicted to be 177°with LFMM but 164°with DFT.

(A) (B)
The RMSD between the three LFMM optimised snapshots and the structures re-optimised using the methods above are reported in Table 1. This shows that PM7 and all DFT methods tested are in good agreement with LFMM. Notably, DFT methods that explicitly account for non-covalent interactions generally show better agreement, with B3LYP giving the largest RMSD values, while B3LYP-D2 is the DFT method that gives structures most comparable to LFMM. Heavy atom RMSD values are slightly smaller than those for all atoms. Further analysis (Tables S4 and S5 in the supporting information) shows that non-coordinated residues make a larger contribution to the overall RMSD than do Cu and coordinated residues (0.36/0.22, 0.61/0.43 and 0.52/0.36 Å for non-coordinated/coordinated residues in snapshots 1, 2 and 3, respectively), such that the minor differences observed are apparently due to mismatches between AMBER and B3LYP-D2 for the peptide, rather than the effects of Cu-coordination. Despite some variance between the methods, RMSD values indicate good agreement against structures calculated using the copper LFMM parameters. This is particularly evident for the systems truncated to include only the copper and coordinating residues, showing very good agreement between DFT and LFMM optimised coordination spheres.
PM7 gives excellent agreement against LFMM optimised structures with very small RMSD values across all three snapshots and the corresponding atom subsets. To validate these results the PM7 optimised structures have been compared against the DFT optimised structures, with RMSD values reported in Table S7 in the supporting information. In general, the difference between the structures is very good with RMSD values ranging from 0.30 to 0.59 Å against B3LYP-D2. This performance of PM7 led us to use this method on further snapshots; a total of 21 snapshots, including snapshots 1, 2 and 3; were collected from the simulation at regular intervals of 0.25 ns starting at 5 ns. RMSD between LFMM and PM7 optimised structures averages 0.20 Å (maximum value = 1.37 Å, minimum = 0.08 Å). The difference between the two methods here is found to be small with very little rearrangement of the LFMM systems when using the semi-empirical method, showing that for this case the LFMM parameters perform very well compared to PM7.
As noted above, RMSD from the starting point indicates that equilibration is reached in approximately 5 ns, but than even beyond this point variations exist. To probe the significance of these, a range of geometrical parameters were calculated and their variation over the second half of the trajectory analysed. The resulting values are summarised in Table 2, indicating that most vary only slightly around average values, with standard deviations close to 0.1 Å. The obvious exception to this pattern is the 'floating' contact between Cu and O from the backbone carbonyl of ASP1, which exhibits a larger standard deviation and a much greater spread between minimum and maximum values. As shown in Figure 5, for the bulk of the 5 ns recorded this distance stays between 2.5 and 3.0 Å, but repeatedly increases to around 4.5 Å. This indicates that O1 can act as a possible fifth ligand, albeit weakly, but that an alternate form with no such weak coordination is also sampled during the LFMD trajectory. A histogram of the same data shows that O1 never comes closer than 2.30 Å to Cu, but spends over 95% of the trajectory within 3.0 Å. It should also be reiterated that the molecular dynamics simulation was carried out with implicit solvation; in the case of explicit solvation a neutral water molecule might act as a weakly coordinated fifth ligand instead of the carbonyl moiety observed here. Transient coordination modes, such as Cu...O1, were previously proposed by Raffa and Rauk (2007) although the apical ligands in their study were carboxylates from ASP and GLU residues. Coordination of oxygen from ASP1 is also observed in 'component I' species of copper Aβ, found at slightly acidic pH Drew et al., 2009). Figure 6 shows two snapshots from the LFMD trajectory with short and long Cu⋯O1 interactions, truncated to the coordination sphere. Within this implementation of LFMM there is no possibility for ligand exchange, such that associative ligand exchange reactions, for instance of ALA2 by ASP1, cannot be modelled. Table 2 reports R g over the final 5 ns of simulation. The average value of 6.6 Å, with a very small standard deviation of 0.04 Å, indicates a relatively compact structure. A low mode molecular dynamics conformational search of free Aβ 1-16 was carried out resulting in an average R g of 7.2 Å. This larger value shows, as would be expected, that the geometrical constraints arising from copper coordination lead to a slightly more compact peptide conformation on average.

Conclusions
A deeper understanding of copper's role within Aβ peptides and its significance in the progression of Alzheimer's disease is of great importance. Experimental studies have broached this, alongside extensive computational simulations but much is still unknown. A potential tool for this application is LFMM, which has already proven successful in calculations on copper containing proteins (Deeth, 2007;Diedrich & Deeth, 2008). It offers greater transferability than inclusion of ad hoc copper parameters within force fields and is much less computationally expensive than ab initio molecular dynamics, whilst still producing results of close to DFT accuracy (Cendic, Matovic, & Deeth, 2013).
We have tested copper LFMM parameters on a set of small molecule test systems with imidazole and formamide representing, N binding from histidine and backbone carbonyl coordination, respectively. Errors between angle and bond lengths of DFT and LFMM optimisations were found to be small, typically less than 0.1 Å or 5°. These parameters were then utilised in LFMD simulations of Cu-Aβ 1-16 with several snapshots extracted from the trajectory for further analysis. Optimisations of the snapshots using several DFT methods and PM7 resulted in good agreement against LFMM geometries. In particular, PM7 showed excellent results with very small RMSD values compared to the LFMM geometry.
Analysis of the LFMD trajectory showed that copper bond lengths and angles remained stable throughout the simulation, with small deviations away from average values. A possible fifth coordination mode was also observed, with the carbonyl group from the ASP1 residue acting as a transient ligand. Penta-coordinate Aβ copper systems have been previously observed in simulations, but usually attributed to carboxylate groups on ASP and GLU residues, Hureau et al., 2009;Raffa & Rauk, 2007) whilst this binding mode is more representative of the component I model of Cu-Aβ.
The overall strong performance of LFMM for copper bound to Aβ here shows that there is great potential in this technique. With good agreement between LFMM and DFT/PM7 geometries for not only test structures but also the Aβ 1-16 peptide, it shows that these parameters should be applicable to larger Aβ fragments as well the full 42 residue peptide and to dimers and beyond.

Disclosure statement
No potential conflict of interest was reported by the authors.