DFT Protocol for EPR Prediction of Paramagnetic Cu ( II ) Complexes and Application to Protein Binding Sites

With the aim to provide a general protocol to interpret electron paramagnetic resonance (EPR) spectra of paramagnetic copper(II) coordination compounds, density functional theory (DFT) calculations of spin Hamiltonian parameters g and A for fourteen Cu(II) complexes with different charges, donor sets, and geometry were carried out using ORCA software. The performance of eleven functionals was tested, and on the basis of the mean absolute percent deviation (MAPD) and standard deviation (SD), the ranking of the functionals for Az is: B3LYP > B3PW91 ~ B3P86 > PBE0 > CAM-B3LYP > TPSSh > BH and HLYP > B2PLYP > MPW1PW91 >ω-B97x-D >> M06; and for gz is: PBE0 > BH and HLYP > B2PLYP > ω-B97x-D > B3PW91~B3LYP~B3P86 > CAM-B3LYP > TPSSh~MPW1PW91 >> M06. With B3LYP the MAPD with respect to A tl z is 8.6% with a SD of 4.2%, while with PBE0 the MAPD with respect to g tl z is 2.9% with a SD of 1.1%. The results of the validation confirm the fundamental role of the second order spin-orbit contribution to Az. The computational procedure was applied to predict the values of gz and Az of the adducts formed by Cu(II) with albumin and two fragments of prion protein, 106–126 and 180–193.


Introduction
The coordination chemistry of copper has been extensively studied for its biological role [1][2][3][4], synthetic and catalytic applications [5][6][7], and recently, for its pharmacological action and the potential use of its compounds in medicine [8][9][10].The most stable oxidation states of copper are +I (3d 10 , diamagnetic) and +II (3d 9 , paramagnetic) [11], which often-and particularly in the living systems and in catalysis-give interconversion between each other.In normal conditions, the state +II is more stable than +I, which disproportionate in aqueous solution, and specific ligands are necessary to stabilize the lowest oxidation state [12].
Due to the non-spherical symmetry and Jahn-Teller effect expected for a d 9 ion, the stereochemistry of Cu(II) complexes is characterized by non-rigid structures and covers several distorted geometries (the phenomenon often named "plasticity") [13,14].Species with coordination numbers 4, 5, and 6, and geometries square planar, square pyramidal, bipyramidal trigonal, compressed and elongated octahedral with all the possible distortions, have been observed [15].Several instrumental techniques are available for the study of Cu(II) compounds, among which UV-Vis and circular dichroism (CD) spectroscopy, but the most used is electron paramagnetic resonance (EPR).
Electron paramagnetic resonance technique can be applied, in principle, to any paramagnetic systems with one or more unpaired electrons; the background of EPR is similar to nuclear magnetic resonance (NMR), but the magnetic field splits the energy of the α and β spin states of the electron instead of those of the nuclei [16].For a mononuclear paramagnetic compound with only one unpaired electron, two spin Hamiltonian parameters characterize its EPR spectrum: the g factor and the hyperfine coupling constant (HFCC), indicated with A, between the unpaired electron and nuclei (g iso and A iso in an isotropic spectrum, and g x , g y , g z and A x , A y , A z in an anisotropic spectrum) [16].For Cu(II) complexes the coupling of the electron is with the nuclei of the isotopes 63 Cu and 65 Cu, whose natural abundance is 69.1% and 30.9%, respectively.The g and A tensors are strongly perturbed by the chemical environment around Cu(II) ion and their variation can give information on the structure of a specific species.In particular, the values of g z and A z are more sensitive to the change in the coordination sphere than the x and y components.In general, with increasing the electron donor capability of the ligands, a decrease of the experimental value of g z and an increase of A z are observed.For example, Peisach and Blumberg [17] found a linear relationship between the experimental g z and A z values for CuO 4 , CuNO 3 , CuN 2 O 2 , CuN 3 O, and CuN 4 equatorial coordination with different electric charge of the species.The presence of ligands in the axial positions causes an increase of g z and a decrease of A z .Moreover, because of the plasticity of copper and the possibility to have several distortion degrees (for example, in an elongated octahedron the distance between copper and the axial ligands can be more or less long) a variation-not always easily predictable-in the values of g z and A z is expected.This results in more complicated situations compared to other metal ions such as vanadium(IV) (3d 1 ), for which the "additivity relationship" allows one to predict the value of A z within 3 × 10 −4 cm −1 [18].
Therefore, methods to help in the interpretation of the experimental spectra and to correlate the values of g z and A z with the structure of a specific Cu(II) compound are desirable.Among the possible tools, density functional theory (DFT) [19] has reached an enormous popularity and a number of reviews have been published over the last years [20][21][22][23][24][25][26].Density functional theory is a theory of the electronic structure based on the electron density ρ(r) rather than the many-electron wavefunction and allows to determine the properties of a system using a functional dependent only on ρ(r) with relatively short computational times.With an appropriate choice of the functional and basis set it is possible to perform calculations on large molecules with high efficiency.Up to now, unfortunately, for most of systems there are no rules to choose the correct DFT protocol to obtain the best agreement with the experimental data: an efficient functional is fundamental but it depends on the molecular/spectroscopic properties and on the metal under examination, the basis set must be large enough and include-if necessary-polarization and diffuse functions and relativistic effects, and in some cases the role of the solvent must be considered through continuum solvation models such as polarizable continuum model (PCM) [27] or solvation model based on density (SMD) [28].Therefore, on the basis of the actual knowledge in the literature, the best protocol for a DFT calculation must be obtained case by case, and a criterion to establish which level of theory will give the best agreement with the experiment under consideration is lacking.Density functional theory methods allows to calculate the structure, energy, and molecular properties with great accuracy.Concerning the spectroscopic behavior of the transition metal ions, IR, UV-vis, CD, and magnetic circular dichroism (MCD) spectra were simulated, but only for nuclear magnetic resonance (NMR) and EPR-and UV-Vis in some rare cases [29]-the agreement with the experimental data can be quantitative [30][31][32][33][34][35][36][37].
For transition metal species, the prediction of A is strongly influenced by the choice of the functional, and it has been shown that for each metal ion it is necessary to validate that which gives the best performance [38].At the moment, a unique functional which works efficiently for the entire series of transition metals is not available.For example, for some iron(III) complexes (e.g., with nitrosyl ligands) the OLYP functional performs better than meta-generalized gradient approximation (meta-GGA) and hybrid functionals such as B3LPY and TPSSh [39]; for Mo(V) species hybrid functionals with 30-40% of exchange admixture give a good agreement with the experimental data [40], and for V(IV)O compounds, the half-and-half hybrid BHandHLYP functional is able to predict A z for neutral species formed by (N,O) ligands with a mean absolute percent deviation (MAPD) from the experimental value lower than 5% [36].Generally, meta-gradient functionals underestimate the value of A and an elevated fraction of the Hartree-Fock (HF) exchange, which increases the spin polarization, is necessary to predict values in line with the experiments.Over the last decade, new meta-GGA, hybrid, and double hybrid functionals have been tested in the prediction of EPR spin Hamiltonian parameters; for instance, the hybrid version of meta-GGA TPSS, i.e.TPSSh, represents a good alternative to the "usual" B3LYP.In the double hybrid functional B2PLYP [41], which consists of B88 exchange combined with the LYP correlation functional through two empirical parameters controlling the amount of HF exchange and PT2 correlation, the HF exchange is quite large with respect to the most used hybrid functionals (53%, compared with 20% of B3LYP and 25% of PBE0, for example); it has been successfully tested on "bare" V(IV) complexes reaching a MAPD value of 3.3% from the experimental values of A [37].
For the calculations of HFCC A, three contributions must be taken into account: the Fermi contact (A FC ), the anisotropic or dipolar hyperfine interaction (A D ), and the second-order spin-orbit (SO) coupling (A SO ) [34].Generally, the most critical parameter is the Fermi contact, which is directly proportional to the spin density at the nucleus (ρ α−β N ), and, for this reason, a good model must be able to predict the spin polarization in the electron distribution of the core and valence shells due to the unpaired electron in the d orbitals [42,43].The SO contribution increases going toward the right of the first transition series and becomes important for late 3d elements, such as Co [44], Ni [45], and mainly Cu [34,44].A SO z is in the approximate range 4-8% with a mean of ca.6% compared to the sum (A FC + A D z ) for V(IV)O [36], while it increases for Cu(II) species; for example, a value of ca.25% has been calculated with BP86 functional for [Cu(acac) 2 ] [34].Therefore, as it was pointed out in the literature, the SO contribution should be included in the A calculation [34,44].Van Lenthe et al. [46] calculated the zero-order regular approximation (ZORA) to the Dirac equation to determine the SO coupling; this approach, that has the limitation that spin polarization and SO effects cannot be introduced contemporaneously, has been overcame by Neese in the coupled-perturbed Kohn-Sham implementation of SO contribution to HFCC [34].The validation of HFCC on a series of 3d transition metal complexes agreed well with the experimental data; the main error was related to the underestimation of the Fermi contact, whose prediction can be improved by varying the functional [34].
For the first row of transition metal complexes, the scalar relativistic (SR) contributions are of little importance or negligible [35], and they must be considered only for the elements of the second and third series, as demonstrated in a systematic study for 4s, 5s, and 6s radicals [47].
Over the last years, several papers have been published on the calculation of g and A values for Cu(II) complexes, but to the best of our knowledge the only systematic study is from Neese [34], who found a MAPD below 5% using B3LYP functional for 10 different species.However, he tested only B3LYP and BP86, and therefore, any comparison with other or more recent functionals is lacking.Larsen and co-workers [48][49][50] dedicated several papers to Cu(II) compounds but the number of cases considered was always below five structures.
In this paper a systematic study on the DFT prediction of g and A values for fourteen Cu(II) complexes with different total charge, coordination geometry, and type of donors is presented.The experimental values of A x , A y , A z , and g z for these complexes are listed in Table S1 of the Supplementary Materials.Eleven functionals were examined and their performance compared with a statistical analysis based mainly on the MAPD value.The most problematic cases were also evaluated in detail and some methods to improve the prediction are also proposed.Finally, the best computational conditions were applied to the calculations of g z and A z values to the species formed upon the binding of Cu(II) ion to the N-terminal region of human serum albumin (HSA) and to fragments 106-126 and 180-193 of prion protein (PrP).
The optimized geometries are shown in Figure 1, while the comparison between the calculated and experimental bond lengths and angles for three selected complexes are summarized in Table 1.Cartesian coordinates of the optimized minima are available in Tables S2-S12 of the Supplementary Materials.As emerged from an examination of the data in Table 1, generally, a good agreement between the simulated and the single-crystal X-ray diffraction structure was found and this ensures that the geometry optimization was reasonably accurate to predict the EPR spin Hamiltonian parameters.The superimposition of the optimized geometry of [Cu(bipy) 2 (NCS)] + , [Cu(ttcn)] 2+ and [Cu(S,S-mnpala)] with the respective X-ray structure is reported in Figures S1-S3 of the Supplementary Materials.
The most critical cases will be discussed in the Section 2.5. 1 Distances in Å and angles in degrees. 2 D stands for donor. 3Experimental data taken from Reference [54]. 4  Experimental data taken from Reference [55]. 5Experimental data taken from Reference [56].
The accuracy of the method was assessed using the mean percent deviation (MPD) and mean absolute percent deviation (MAPD) from the experimental value as in Equations ( 1) and ( 2), considered in the literature good quality parameters [71][72][73].The value of the pure percent deviation (PD) was also taken into account with the aim to identify possible overestimation or underestimation of the predicted spin Hamiltonian EPR parameters by the specific functional.Another parameter considered to rank the performance of the functionals was the standard deviation (SD) of the data.

Prediction of A z
The values of the MAPD obtained for each functional, tested on the fourteen Cu(II) complexes, are summarized in Table 2 and Figure 2. From an analysis of the data, it emerges that the best performance was achieved by the three-parameters hybrid functionals B3P86 (20% HF exchange), B3PW91 (25% HF exchange), B3LYP (20% HF exchange), PBE0 (25% HF exchange), and the long-range-corrected CAM-B3LYP (short range HF exchange 19%, long range HF exchange 65%).Worst performance was obtained with TPSSh (10% HF exchange), BHandHLYP (50% HF exchange), MPW1PW91 (25% HF exchange), the double hybrid B2PLYP (53% HF exchange), and the range separated ω-B97x-D (short range HF exchange 22.2%, long range HF exchange 100%).Finally, the M06 global hybrid functional (27% HF exchange), in contrast to its good performance in the prediction of thermochemistry, kinetic, chemical structures, and non-covalent interactions [74,75], cannot be suggested to simulate spectroscopic data, confirming previous published data [29].Considering the strong dependency of A z on the core-and valence shell electron polarization, the reason could be found in the low accuracy of the functionals with a large number of empirical parameters in the prediction of the electronic density around the nucleus with respect to the functionals less parameterized, as recently reported by Medvedev et al. [76].Taking into account also the error dispersion, represented by the standard deviation, it can be concluded that the highest accuracy is reached by the B3LYP functional with lowest MAPD from the experimental values of 8.6% and an associated SD of 4.2%.Summarizing the results, the ranking of the functionals is: 2).It can be noticed that from an analysis of the MPD values, B3P86 and B3PW91 resulted in an overestimation of A z , while B3LYP and PBE0 in an underestimation of its values; this is in agreement with the previous published data on the Cu(II)-peptides systems [77].These data confirm that a percent amount of HF exchange in the functionals is necessary even if a correlation between the percentage of HF and A calcd z cannot be found.It can be interesting to evaluate the contribution of A SO z to A calcd z for the Cu(II) species examined.While it is around 6% for V(IV)O complexes, it is significantly higher for copper(II); this is the reason why Gaussian software-which neglects the second-order term A SO z -simulates well A z of V(IV)O [36] and not of Cu(II) species.From the data of this study, the mean percent contribution of SO to A calcd z is 19.1% with a maximum value of 24.4% for [Cu(bipy) 2 (NCS)] + .This is in line with the data previously reported [34], suggesting that the calculation of A SO z is essential to accurately predict A z .For our Cu(II) species in Figure 1, the MAPD for the best functional B3LYP would increase from 8.6% to 32.9% without the inclusion of the SO coupling.the data previously reported [34], suggesting that the calculation of SO z A is essential to accurately predict Az.For our Cu(II) species in Figure 1, the MAPD for the best functional B3LYP would increase from 8.6% to 32.9% without the inclusion of the SO coupling.

Prediction of gz
The values of MAPD for the prediction of gz are summarized in Table 3 and graphically represented in Figure 2.They show a different behavior with respect to the prediction of Az.Overall, the results lie in a narrow range (4.6-2.7%),so slight differences in the prediction of gz are expected changing the functional.The values of SD are less scattered as well.The highest accuracy is achieved by the hybrid BHandHLYP and PBE0, the double hybrid B2PLYP, and the range separated ω-B97x-D (MAPD around 3%), while for Becke hybrid functionals B3P86, B3PW91, and B3LYP, and for the long-range-corrected CAM-B3LYP values around 3.5% are obtained.Among the best functionals, BHandHLYP slightly overestimates gz, and PBE0, B2PLYP, and ω-B97x-D underestimates it.The worst performance was achieved by TPSSh, MPW1PW91, and M06 functionals.Taking into account also the standard deviation, it can be concluded that the best functional is PBE0, 2.9% MAPD and 1.1% SD.The ranking of the functionals in the prediction of gz is: PBE0 > BHandHLYP > B2PLYP > ω-B97x-D > B3PW91~B3LYP~B3P86 > CAM-B3LYP > TPSSh~MPW1PW91 >> M06 (see Table 3).

Prediction of g z
The values of MAPD for the prediction of g z are summarized in Table 3 and graphically represented in Figure 2.They show a different behavior with respect to the prediction of A z .Overall, the results lie in a narrow range (4.6-2.7%),so slight differences in the prediction of g z are expected changing the functional.The values of SD are less scattered as well.The highest accuracy is achieved by the hybrid BHandHLYP and PBE0, the double hybrid B2PLYP, and the range separated ω-B97x-D (MAPD around 3%), while for Becke hybrid functionals B3P86, B3PW91, and B3LYP, and for the long-range-corrected CAM-B3LYP values around 3.5% are obtained.Among the best functionals, BHandHLYP slightly overestimates g z , and PBE0, B2PLYP, and ω-B97x-D underestimates it.The worst performance was achieved by TPSSh, MPW1PW91, and M06 functionals.Taking into account also the standard deviation, it can be concluded that the best functional is PBE0, 2.9% MAPD and 1.1% SD.The ranking of the functionals in the prediction of g z is: PBE0 > BHandHLYP > B2PLYP > ω-B97x-D > B3PW91~B3LYP~B3P86 > CAM-B3LYP > TPSSh~MPW1PW91 >> M06 (see Table 3).

Effect of the Geometry Optimization on the Prediction of A z
For sake of completeness, the effect of the geometry optimization on the prediction of A z was also investigated focusing on the cases for which percent deviations larger than 10.0% were obtained using the best functional B3LYP (see Table S1 of the Supplementary Materials) [78].In particular, the structure of [Cu(H -2 GGH)] − (PD = −13.8%),[Cu(salpn) 2 ] (PD = 11.5%), and [Cu(mnt) 2 ] 2− (PD = −11.8%)was re-optimized following two different strategies: i) increasing the number of basis functions using the def2-TZVP basis set for all the atoms, and ii) considering the solvation effects through the application of the SMD continuum model [28] during the optimization with def2-TZVP basis set.As reported in the Introduction, the Cu(II) coordination geometry exhibits a high degree of complexity due to the plasticity effect [13,14], and for this reason, the simulation of some species presents several difficulties.Since EPR HFCC's are highly sensitive on the metal environment, an incorrect optimization could result in the bad prediction of A z .Both the complexes [Cu(H -2 GGH)] − and [Cu(salpn) 2 ] should reasonably have a square planar geometry, but the flexibility of His side chain or the propanediamine backbone could give several plane distortions; instead [Cu(mnt) 2 ] 2− has a geometry intermediate between the square planar and the tetrahedral [79], which could vary depending on the solvent used to record the spectra.The deviation from square planar toward tetrahedral geometry can be expressed by the angle (θ) between the two planes defined by the metal (M) and the donor atoms D 1 and D 2 on one hand, and M, D 3 , and D 4 on the other: θ is 0 • for the square planar and 90 • for the tetrahedron (Scheme 1).In the case of the [Cu(H -2 GGH)] − , [Cu(salpn) 2 ], and [Cu(mnt) 2 ] 2− , θ is 22.7 • , 27.5 • , and 39.2 • , respectively.

Effect of the Geometry Optimization on the Prediction of Az
For sake of completeness, the effect of the geometry optimization on the prediction of Az was also investigated focusing on the cases for which percent deviations larger than 10.0% were obtained using the best functional B3LYP (see Table S1 of the Supplementary Materials) [78].In particular, the structure of [Cu(H-2GGH)] − (PD = −13.8%),[Cu(salpn)2] (PD = 11.5%), and [Cu(mnt)2] 2− (PD = −11.8%)was re-optimized following two different strategies: i) increasing the number of basis functions using the def2-TZVP basis set for all the atoms, and ii) considering the solvation effects through the application of the SMD continuum model [28] during the optimization with def2-TZVP basis set.As reported in the Introduction, the Cu(II) coordination geometry exhibits a high degree of complexity due to the plasticity effect [13,14], and for this reason, the simulation of some species presents several difficulties.Since EPR HFCC's are highly sensitive on the metal environment, an incorrect optimization could result in the bad prediction of Az.Both the complexes [Cu(H-2GGH)] − and [Cu(salpn)2] should reasonably have a square planar geometry, but the flexibility of His side chain or the propanediamine backbone could give several plane distortions; instead [Cu(mnt)2] 2− has a geometry intermediate between the square planar and the tetrahedral [79], which could vary depending on the solvent used to record the spectra.The deviation from square planar toward tetrahedral geometry can be expressed by the angle (θ) between the two planes defined by the metal (M) and the donor atoms D1 and D2 on one hand, and M, D3, and D4 on the other: θ is 0° for the square planar and 90° for the tetrahedron (Scheme 1).In the case of the [Cu(H-2GGH)] − , [Cu(salpn)2], and [Cu(mnt)2] 2− , θ is 22.7°, 27.5°, and 39.2°, respectively.The capability of the level of DFT theory to simulate the distortions could result in a good or poor prediction of Az.The re-optimization of the geometry using the extended def2-TZVP basis set Scheme 1. Definition of the square planar and tetrahedral geometry based on the angle θ.
The capability of the level of DFT theory to simulate the distortions could result in a good or poor prediction of A z .The re-optimization of the geometry using the extended def2-TZVP basis set gives a significant improvement of A calcd z for [Cu(salpn) 2 ] that reaches a deviation of 6.9% from the experimental value.Applying the def2-TZVP set and continuum solvation model (in CHCl 3 ), the PD of A z for [Cu(H -2 GGH)] − decreases −13.8% to −9.9%.For both the cases (Figures S3 and S4 of the Supplementary Materials), the increase of the optimization theory level reduces the distortion of the equatorial plane: the structures approach the square planar geometry (θ = 11.8 • and 2.8

8%).
As a final attempt to improve the prediction of A z for [Cu(mnt) 2 ] 2− , the A tensor was calculated from the experimental X-ray structure; in this case the deviation of A calcd z from A exp tl z lowers to −6.8%.These data indicate that the limitation of the method could be ascribed to the quality of the optimization.The structure optimization of transition metal compounds has been the object of excellent papers (see, for example, References [52,53,80]) and a more detailed discussion is out of the scope of the present paper.In general, for an increase of the optimization theory level an improvement in the calculation of A z is expected as a consequence of the better quality of the geometry prediction.

Application to the Interaction of Cu(II) with Human Serum Albumin and Prion Protein
Human serum albumin (HSA) is the most abundant protein in the blood (~0.63 mM) [81], with a number of physiological roles, such as the transport of metabolites (for example, the fatty acids and bilirubin) and the maintenance of the osmotic pressure.Furthermore, albumin plays a major role in the transport of metal ions, as Cu(II) [82], Zn(II) [83,84], and Ni(II) [82,85], Co(II) [85], and Cd(II) [85].For example, 5-10% of serum copper and 95% of nickel are bound to HSA [82].Albumin contains two metal binding sites: the N-terminal one (NTS or site I) [82] and the multi-metal binding site (MBS or site II) [86].While MBS is the primary binding site for Zn(II) that is coordinated by His67, Asn99, His247, and Asp249 [87], NTS site binds Cu(II) and Ni(II) with the donor set (NH 2 , N − , N − , His-N)-with NH 2 representing the terminal amino, N − the deprotonated amide groups of residues 2 and 3, and His-N the imidazole nitrogen of His3-in a square planar geometry with little degree of distortion.Even if the first two residues are different, bovine serum albumin (BSA) and rat serum albumin (RSA) binds Cu(II) with the same coordination mode since the histidine residue in third position is conserved (AspAlaHis for HSA, AspThrHis for BSA, and GluAlaHis for RSA).For other species, for example dog, porcine, and chicken serum albumin, His3 is replaced by Tyr and Cu(II) is not bound at the N-terminal region of the protein [82].Taking advantage of the results obtained in this study, the EPR spectrum of the system Cu(II)/HSA has been recorded at pH 7.4 and compared with that calculated for the model including the first five amino acids of HSA (AspAlaHisLysSerNHCH 3 ); B3LYP and PBE0 functionals were used to calculate A and g tensors with the procedure described above.
The optimized structure (with B3LYP-D3 functional and the same basis set used for the simple coordination compounds) of the adducts between Cu(II) and pentapeptide is represented in Figure 3.As it is possible to observe, the coordination consists of four nitrogen donors with His3 nitrogen which closes the last chelate ring; the structure approaches the square planar geometry and the distances Cu-NH 2 , Cu-N − , Cu-N − , Cu-N(His3) are 2.085, 1.944, 1.943, and 2.005 Å, respectively.Experimental and calculated spin Hamiltonian parameters are reported in Table 4; the value of g calcd z is comparable to g exp tl z with a deviation of −1.7%, while that of A calcd z is almost perfect (PD of −1.1%).The values of A are reported in 10 −4 cm −1 units. 2 Experimental parameters measured on the system  The values of A are reported in 10 −4 cm −1 units. 2 Experimental parameters measured on the system Cu(II)/HSA in this work. 3Calculated parameters using B3LYP functional for A and PBE0 for g tensor. 4Experimental g z and A z values measured in Reference [88] for the species [Cu(H -1 L)] + formed by the fragment PrP(106-126) with (Met109-S, N − , N − , His111-N) coordination. 5Values not reported. 6Experimental g z and A z values measured in Reference [88] for the species [Cu(H -2 L)] formed by the fragment PrP(106-126) with (N − , N − , N − , His111-N) coordination. 7Experimental g z and A z values measured in Reference [89] for the species [Cu(H -2 L)] formed by the fragment PrP(180-193) with (N − , N − , N − , His187-N) coordination.
The EPR spectrum generated through WinEPR SimFonia software [90] with the calculated parameters is compared with the experimental one in Figure 4; it is worth noting the good agreement except for the shift toward high fields of the center of the spectrum determined by the value of g calcd z (2.135) smaller than g exp tl z (2.173).If the spectrum is generated using the experimental g z (trace c of Figure 4), the agreement is excellent.Prion proteins (PrP c ) are normal cellular proteins involved in various diseases known as transmissible spongiform encephalopathies which include Creutzfeldt-Jakob disease in humans, bovine spongiform encephalopathy in cattle, and scrapie in sheep and goats [91,92].In these pathologies, prion protein undergoes a misfolding in the C-terminal region (128-231) to scrapie form (PrP Sc ) in a process in which the α-helical folded of PrP c becomes a β-sheet forming aggregates [93].One of the most interesting properties of PrP c is its capability to bind copper [94,95]; in particular, four Cu ions are bound in the "octarepeat" region (60-91) [96,97], and two additional ions in the region 92-111 upon the interaction with His96 and His111 [98,99].It has been proposed that the Cu binding in the non-octarepeat region is related to prion misfolding [100,101].In particular, the regions 106-126 and 180-193 could be involved in this transformation [88,102].
Over the last ten years, DFT methods were applied to calculate g and A tensors for the octarepeat .With the asterisk the superhyperfine coupling between the unpaired electron and 14 N nuclei (not generated with WinEPR SimFonia software [90]) is shown.
Prion proteins (PrP c ) are normal cellular proteins involved in various diseases known as transmissible spongiform encephalopathies which include Creutzfeldt-Jakob disease in humans, bovine spongiform encephalopathy in cattle, and scrapie in sheep and goats [91,92].In these pathologies, prion protein undergoes a misfolding in the C-terminal region (128-231) to scrapie form (PrP Sc ) in a process in which the α-helical folded of PrP c becomes a β-sheet forming aggregates [93].One of the most interesting properties of PrP c is its capability to bind copper [94,95]; in particular, four Cu ions are bound in the "octarepeat" region   [96,97], and two additional ions in the region 92-111 upon the interaction with His96 and His111 [98,99].It has been proposed that the Cu binding in the non-octarepeat region is related to prion misfolding [100,101].In particular, the regions 106-126 and 180-193 could be involved in this transformation [88,102].
Over the last ten years, DFT methods were applied to calculate g and A tensors for the octarepeat region with good results [103,104].Here we applied the computational procedure validated to confirm the coordination mode of Cu(II) with fragments 106-126 and 180-193.In such regions, the residues His111 and His187 play a pivotal role and favor the copper binding [105].In the region 106-126 Di Natale et al. [88] found two species at physiological pH: [Cu(H -1 L)] + with (Met109-S, N − , N − , His111-N) coordination, where N − are the two deprotonated peptide groups of Lys110 and His111, and [Cu(H -2 L)] with (N − , N − , N − , His111-N) coordination, where the binding of deprotonated peptide groups of Met109, Lys110, and His111 adds to His111 imidazole nitrogen.We optimized the structure of the Ac-AsnMetLysHisMet-NHCH 3 peptide (region 108-112), and subsequently, calculated g and A. The two structures optimized with the functional B3LYP-D3 are shown in Figure 5.The percent deviation of g calcd z and A calcd z from the experimental values is −2.9% and −8.7% for [Cu(H -1 L)] + and −2.5% and −1.9% for [Cu(H -2 L)] (Table 4).The data are in line with the validation on the fourteen Cu(II) structures and confirm the attribution of Reference [88].The large deviation of A z for [Cu(H -1 L)] + must be related to the Met-S binding (distance 2.541 Å), which causes a significant distortion of the equatorial plane of Cu(II) ion (θ = 31.2• ).
Concerning the fragment 180-193, the region with five amino acids around His187 was optimized: Ac-IleLysGlnHisThr-NHCH 3 (region 184-188, Figure 6).The residue His187 plays a pivotal role as an anchoring site for the metal binding; its coordination to Cu(II) favors the deprotonation of the peptide groups of Lys185, Gln186   The comparison between the four optimized models for HSA and PrP with the X-ray structures of similar Cu(II)-peptide complexes is presented in Table S13 of the Supplementary Materials.The comparison between the four optimized models for HSA and PrP with the X-ray structures of similar Cu(II)-peptide complexes is presented in Table S13 of the Supplementary Materials.The comparison between the four optimized models for HSA and PrP with the X-ray structures of similar Cu(II)-peptide complexes is presented in Table S13 of the Supplementary Materials.

Computational Details
For all the analyzed coordination compounds and protein models, geometry optimization and harmonic frequency were computed using Gaussian 09 software [67] at DFT theory level.The hybrid

Computational Details
For all the analyzed coordination compounds and protein models, geometry optimization and harmonic frequency were computed using Gaussian 09 software [67] at DFT theory level.The hybrid Becke three-parameter B3LYP functional [57,58] coupled with the Grimme's D3 dispersion was used combined with the split-valence Pople basis set 6-311g(d,p) for the main group elements, while Stuttgart-Dresden (SDD) implemented with f -functions and pseudo-potential was applied for Cu.These conditions have been successfully applied and discussed in the literature for the geometry prediction of first-row transition metal complexes [52,53].Cartesian coordinates of the fourteen Cu(II) complexes examined were generated by EsiGen software [106] (see Supplementary Materials).
The ORCA program [107,108] was used to calculate g and A tensors of 63 Cu center for each complex using the method developed and implemented into the package.The A tensor was obtained as a sum of the three contributions: the isotropic Fermi contact (A FC ), the anisotropic dipolar (A D x,y,z ), and the spin-orbit coupling term (A SO x,y,z ).To achieve an accurate estimation for A FC , the appropriate spin polarization needed to be included using the unrestricted Kohn-Sham (UKS) formalism.
It must be noticed that for the species examined A z value was usually negative, but in the literature its absolute value is usually reported; this procedure was also followed in this work.

Theory Background
The theoretical background and the formalism used by ORCA are well reported in References [34,43,108,109], and here we present only a brief summary.The HFCC in a Cu(II) EPR spectrum comes from the interaction between the spin angular momentum of the electron (S = 1/2) with the spin angular momentum of 63 Cu nucleus (I = 3/2, 69.15% natural abundance).The contributions to 63 Cu HFCC tensor A are: functionals BHandHLYP and B2PLYP, respectively, perform better than the others, and unlike Cu(II) species, the prediction of the second-order spin-orbit contribution is often not necessary.Therefore, the best combination of functional and basis set must be reached case by case and, up to today, there are no rules to know in advance which level of theory could give the best agreement with the experimental spectroscopic data.For Cu(II) complexes examined in this paper, the functionals B3LYP and PBE0 show the best performance in the calculation of A z and g z , respectively.They can be coupled with a triple-ζ basis set such as 6-311g(d,p) to obtain a satisfactory agreement with the experimental data.These computational conditions could be applied to clarify also the coordination environment of Cu proteins or in other biological systems.Finally, we would like to highlight that, when the quality of the optimization is not good, the prediction of A z worsens significantly, and a wider basis set must be used or the solvent effect must be considered.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2312-7481/4/4/55/s1,Table S1: Values of A x , A y , A z and g z for the fourteen Cu(II) complexes examined, Tables S2-S12: Values of A iso , A x , A y , A z and g z and percent deviation (PD) calculated at the level of theory Functional/6-311g(d.p), Table S13: Selected bond distances and angles for the models of HSA and PrP, Figures S1-S5 Funding: This research was funded by Spanish MINECO (project CTQ2017-87889-P), Universitat Autònoma de Barcelona (Ph.D. of Giuseppe Sciortino), Fondazione di Sardegna (project FdS15Garribba) and FFABR 2017 (Fondo per il finanziamento delle attività base di ricerca).

Figure 1 .
Figure 1.Geometry of the fourteen Cu(II) complexes simulated with B3LYP-D3 functional combined with triple-ζ Pople's 6-311g(d,p) basis set (main group elements) and Stuttgart-Dresden (SDD) implemented with f-functions and pseudo-potential (Cu).The numbers of the donors reported inTable 1 are also shown.

Figure 1 .
Figure 1.Geometry of the fourteen Cu(II) complexes simulated with B3LYP-D3 functional combined with triple-ζ Pople's 6-311g(d,p) basis set (main group elements) and Stuttgart-Dresden (SDD) implemented with f -functions and pseudo-potential (Cu).The numbers of the donors reported inTable 1 are also shown.

Magnetochemistry 2018, 4 , 20 Figure 2 .
Figure 2. Mean absolute percent deviation (MAPD) from the experimental values of Az (left) and gz (right) determined with the eleven functionals on the fourteen Cu(II) complexes included in the dataset.It can be interesting to evaluate the contribution of SO z A

Figure 2 .
Figure 2. Mean absolute percent deviation (MAPD) from the experimental values of A z (left) and g z (right) determined with the eleven functionals on the fourteen Cu(II) complexes included in the dataset.

Scheme 1 .
Scheme 1. Definition of the square planar and tetrahedral geometry based on the angle θ.
• for [Cu(H -2 GGH)] − and [Cu(salpn) 2 ], respectively) and an improvement in the calculation of A z is observed.In contrast, for [Cu(mnt) 2 ] 2− (Figure S5) very slight geometry variations are obtained with the change of the basis set or the application of SMD model (θ = 40.1 • and 39.5 • , respectively) and a not appreciable change in the A calcd z occurs (−11.7% vs. −11.

Figure 3 .
Figure 3. Optimized structure at the level of theory B3LYP-D3 of the adduct formed by Cu(II) with the pentapeptide AspAlaHisLysSer-NHCH3 (region 1-5 of HSA).The Cu-donor distances (in Å) are also indicated.Hydrogen atoms are omitted for clarity.

Figure 3 .
Figure 3. Optimized structure at the level of theory B3LYP-D3 of the adduct formed by Cu(II) with the pentapeptide AspAlaHisLysSer-NHCH 3 (region 1-5 of HSA).The Cu-donor distances (in Å) are also indicated.Hydrogen atoms are omitted for clarity.

Figure 4 .
Figure 4. EPR spectra of the system Cu(II)/HSA: (a) spectrum generated using the spin Hamiltonian parameters calculated by DFT methods and reported in Table 4; (b) experimental spectrum and (c) spectrum generated using A calculated by DFT methods and exp tl z g.With the asterisk the superhyperfine coupling between the unpaired electron and 14 N nuclei (not generated with WinEPR SimFonia software[90]) is shown.

Figure 4 .
Figure 4. EPR spectra of the system Cu(II)/HSA: (a) spectrum generated using the spin Hamiltonian parameters calculated by DFT methods and reported in Table 4; (b) experimental spectrum and (c) spectrum generated using A calculated by DFT methods and g exp tl z , and His187 to give a [Cu(H -2 L)] species with coordination (N − , N − , N − , His187-N) coordination.The values of of g calcd z and A calcd z are 2.148 and 191.5 × 10 −4 cm −1 , very close to those experimentally measured (PD of −1.7% and −2.3%) [89].
pivotal role as an anchoring site for the metal binding; its coordination to Cu(II) favors the deprotonation of the peptide groups of Lys185, Gln186, and His187 to give a [Cu(H-2L)] species with coordination (N − , N − , N − , His187-N) coordination.The values of of calcd z 10 −4 cm −1 , very close to those experimentally measured (PD of −1.7% and −2.3%)[89].

Table 1
are also shown.

Table 1
are also shown.

Table 2 .
Mean absolute percent deviation (MAPD) and mean percent deviation (MPD) from the experimental value A z of the fourteen Cu(II) complexes, and standard deviation (SD) as a function of the functional 1 .

Table 3 .
Mean absolute percent deviation (MAPD) and mean percent deviation (MPD) from the experimental value g z of the fourteen Cu(II) complexes, and standard deviation (SD) as a function of the functional 1 .

Table 4 .
Experimental spin Hamiltonian parameters for the interaction of Cu(II) with HSA and fragments 106-126 and 180-193 of prion protein, and comparison with the values calculated by DFT methods.

Table 4 .
Experimental spin Hamiltonian parameters for the interaction of Cu(II) with HSA and fragments 106-126 and 180-193 of prion protein, and comparison with the values calculated by DFT methods.