Converged Structural and Spectroscopic Properties for Refined QM/MM Models of Azurin.

Blue copper proteins continue to challenge experiment and theory with their electronic structure and spectroscopic properties that respond sensitively to the coordination environment of the copper ion. In this work, we report state-of-the art electronic structure studies for geometric and spectroscopic properties of the archetypal "Type I" copper protein azurin in its Cu(II) state. A hybrid quantum mechanics/molecular mechanics (QM/MM) approach is used, employing both density functional theory (DFT) and coupled cluster with singles, doubles, and perturbative triples (CCSD(T)) methods for the QM region, the latter method making use of the domain-based local pair natural orbital (DLPNO) approach. Models of increasing QM size are employed to investigate the convergence of critical geometric parameters. It is shown that convergence is slow and that a large QM region is critical for reproducing the short experimental Cu-SCys112 distance. The study of structural convergence is followed by investigation of spectroscopic parameters using both DFT and DLPNO-CC methods and comparing these to the experimental spectrum using simulations. The results allow us to examine for the first time the distribution of spin densities and hyperfine coupling constants at the coupled cluster level, leading us to revisit the experimental assignment of the 33S hyperfine splitting. The wavefunction-based approach to obtain spin-dependent properties of open-shell systems demonstrated here for the case of azurin is transferable and applicable to a large array of bioinorganic systems.

Introduction

Copper is the second most abundant transition metal in biological systems,[1] and Cu-containing enzymes are known to catalyze a variety of reactions, in addition to being involved in electron transfer processes. The copper centers can be categorized based on geometry and coordination.[2,3] Type I Cu centers, also called blue copper centers, feature an intense absorption at around 600 nm. Type II or “normal” low-molecular-weight copper coordination compounds lack this absorption.[4] Besides these, dimeric type III copper proteins and artificial classes (e.g., “type-0” copper[5]) exist.[6] Aside from their UV/vis spectra, these types can be distinguished by electron paramagnetic resonance (EPR) spectroscopy.[7] Because of the unusual active-site geometry compared to synthetic, tetrahedral, or square planar Cu(II) complexes, the blue copper proteins gave rise to the question if the protein is following the active site, or if the active site geometry is dictated by the protein matrix (entatic state principle).[8−11]

Azurin is a representative type I copper protein with a distorted trigonal bipyramidal active-site geometry.[12] It mediates one-electron transfer in bacteria by switching the oxidation state of copper between Cu(I) and Cu(II). The copper ion is ligated by a cysteinate (Cys112) and two histidine residues (His46 and His117) in the ligand plane, as well as a weakly bound methionine (Met121) and the backbone carbonyl of Gly45 as axial ligands (Figure ). This active site geometry, specifically the thiolate sulfur coordination, gives rise to unique spectral features,[13] which are known to be modulated by the covalency of the Cu–SCys bond. The origin of the bright absorption in the Cu(II) state has been clarified in pioneering in-depth studies by Solomon and co-workers and is attributed to ligand-to-metal charge transfer (LMCT) involving the cysteine pπ and the half-occupied copper d-based molecular orbitals.[14,15] This is contrary to molecular, normal, copper complexes, often represented by the square planar CuCl4.[16] Their absorption spectrum is dominated by transitions from the ligand pσ orbitals to copper d, while the π transition is of smaller intensity—an inverted intensity pattern compared to blue copper proteins.[17]

Figure 1

Structure of azurin (a), schematic depiction of the copper site including the SOMO (b), and molecular orbital picture (c) showing the name-giving LMCT absorption.

The strong LMCT bands in azurin indicate a strongly covalent bond between copper and cysteine sulfur. The quantification of the covalency at the copper center has been the subject of various studies utilizing a wide range of spectroscopic techniques. One method is sulfur K-edge X-ray absorption spectroscopy (XAS), where the intensity of the pre-edge region is determined by the 3p/3d mixing.[18,19] Another approach that gives a broader picture of the coordination environment of copper is EPR spectroscopy and associated techniques. The hyperfine coupling constants (HFCs) obtained from these experiments report on the interaction of the nuclear spin of a specific isotope with the electron spin of the system. In addition to copper,[20] the histidine nitrogens Nδ and Nε, and hydrogen HFCs have been measured.[21,22] Several studies also focused on the HFCs of the cysteine β-hydrogens (HA, HB) to elucidate the distribution of spin density on the cysteine ligand.[22−25] Overall, diverse interpretations of experimental data are encountered in different studies.

Using XAS Cu K-edge spectroscopy,[18,19] the pre-edge intensity indicates how much s character is mixed into the Cu 3d orbital. To extract the covalency of the ligand–metal bond, in the sense of the amount of ligand character, the pre-edge area is fitted and compared to a reference.[26] In wild-type azurin, the singly occupied molecular orbital (SOMO) was determined to have 40 ± 3% sulfur character.[18] Since there are two sulfur-containing ligands to the copper (Cys112 and Met121), the experiment was redone using selenomethionine (M121SeM) mutants, showing a sulfur character of 37.5 ± 3%, which is arising only from the Cys112 sulfur. Hence, the covalency of the Cu–SCys112 bond was characterized as 37.5% covalency at the sulfur. Using point mutations, also the second coordination sphere was investigated.[19] In that study, in wild-type azurin, a sulfur character of 45 ± 3% was determined. Removing the hydrogen bond to cysteine sulfur in an F114P mutant leads to an increase of the sulfur character to 54 ± 3%. Further changes in the second solvation sphere, which have been attributed to changes in the electrostatic environment of the copper center, lead to a decrease in the sulfur character, to 31 ± 3% for N47S and 43 ± 3% for F114N azurin, respectively.

In contrast to the isolated fitting process involved in the determination of covalency using XAS, EPR spectroscopy requires fitting of spin Hamiltonian parameters for the whole copper ligand system. This leads to a steep increase in complexity. Therefore, in the latest electron–electron double resonance detected nuclear magnetic resonance (EDNMR) study,[27] the results from quantum mechanics/molecular mechanics (QM/MM) calculations,[28] one-dimensional (1D) EDNMR, and two-dimensional (2D) EDNMR are compared. Unfortunately, this leads to different sets of spin Hamiltonian parameters, with the sulfur HFC differing by 10 MHz. Also, to improve the fit, the principal values of the 14N HFCs have been varied, leading to an increase of 2 MHz with respect to previous experiments.[21] To extract the spin population from the 33S HFCs empirical comparisons were used,[29] which yield a total sulfur spin population of 29.1–30.4%, a value that suggests significantly lower covalency of the Cu–S bond compared to 38% inferred from XAS.[18,19] It is also much lower than the spin population determined by previous spectroscopy-connected computational studies (36–62%).[28,30−32]

To rationalize these differences in interpretation, the concepts involved in the transfer from the experiment to a descriptor of the electronic structure, such as the spin population are reviewed. There are three main conceptual points here: (1) The experimental procedure and post-treatment, especially in the case of the HFCs for the fitting procedure of the spin Hamiltonian (SH) parameters; (2) extraction of the spin population from the experiment, in the case of the HFCs from the SH parameters; and (3) the concept of spin population and its generality. As pointed out previously,[33,34] spin populations are not physical observables, and hence there is not necessarily a uniquely best definition. The spin density distribution is a three-dimensional (3D) function of space and is an observable property from which other observables like hyperfine couplings can be unambiguously deduced. The spin population is a more or less arbitrary assignment of spin density to individual atoms in the system according to some prescription. It also takes no notice of the actual radial shape of the spin density distribution and rather represents an integrated quantity. Consequently, the relationship between spin population and physical observables involves approximations that are different for each spectroscopic method. Hence, it is not surprising when the values for spin populations deduced from different experimental methods do not match. Nevertheless, we use spin populations for qualitative interpretative purposes. In particular, in the case of the type 1 copper site, values of the sulfur spin population deduced from density functional theory (DFT) calculations tend to be much higher (50%) than values deduced either from XAS (38%)[19] or EPR/ENDOR (30%).[27] This is a significant discrepancy, on which we hope to shed more light in the present study that uses state-of-the-art quantum chemical modeling of the type 1 copper site in Azurin.

The first step in obtaining an accurate picture of the electronic structure at the copper center is to describe the underlying geometric structure as accurately as possible. The Cu–SCys distance is an important structural parameter that has been ill-defined in crystallographic studies, which historically present a wide distribution of values from 2.20 to 2.30 Å. Extended X-ray absorption fine structure (EXAFS) studies[35] suggested a shorter Cu–S distance and the most recent EXAFS study indicates a length of 2.12 Å for the Cu–S bond,[36] shorter than most crystallographic models. Computational studies have addressed many of the structural and spectroscopic properties of copper proteins using a variety of theoretical approaches, ranging from semiempirical[14,37,38] to multireference methods[10,39,40] and excited-state dynamics.[41−43] The protein has been sometimes treated using molecular mechanics[44,45] or QM/MM.[28,46−48] Spectroscopic properties have been calculated using density functional theory (DFT). One example is the azurin g-tensor and hyperfine couplings.[32] In a QM/MM framework, optical and X-ray absorption spectra as well as hyperfine coupling constants have been determined.[28] However, methods more accurate and reliable than DFT were not generally feasible for realistic models of azurin that include the protein environment.

Here we revisit the geometric and electronic structure and properties of azurin in the Cu(II) oxidation state using a QM/MM approach to include an explicit description of the protein environment. The complete solvated protein is included in the MM region. To ensure the best possible geometries are obtained, the performance of different DFT functionals is first benchmarked against coupled cluster theory for structural parameters. This is possible thanks to the domain-based local pair natural orbital implementation of coupled cluster singles doubles and perturbative triples, DLPNO-CCSD(T),[49] that is available for open-shell systems.[50] The convergence of spin densities is carefully studied and ligand HFCs are computed by both DFT and DLPNO-CCSD methods.[51] The results derived from a sequence of QM/MM models describe how the spin density distribution responds to specific components of the protein, while the spectroscopic parameters are analyzed in detail and compared with experimental data to resolve the conflicting interpretations of experimental studies regarding the Cu–S covalency and provide insight into the structure–property correlations for blue copper sites.

Models and Methods

Construction of the QM/MM Model

The MM models are based on an average geometry of the four tetramers in the 4AZU[52] structure.[53] Starting from the crystal structure, protons were added for pH 7 using the CHARMM36 force-field parameters.[54] Standard protonation states were assumed for all residues. The histidines were assumed to be protonated in the Nε position, unless indicated otherwise by nearby residues. TIP3P water molecules were added to the protein, forming a sphere with at least 5 Å of buffering between the protein surface and the surface of the sphere. This adds 1459 water molecules in total. Hydrogen positions were minimized with a 10 000-step conjugate gradient algorithm. Water positions were relaxed using a 10 000-step NVT simulation at 100 K, while the oxygens of the outer water layer were kept fixed at their initial positions to prevent “breathing” motions of the system. Finally, both hydrogen and water positions were optimized using a 10 000-step conjugate gradient algorithm, with the same constraints employed.

QM/MM calculations are performed using electrostatic embedding, with linking hydrogen atoms to saturate the QM area. A set of models between 57 and 245 atoms were used, to study the changes of properties with increasing system size. The number of QM atoms was systematically increased by adding functional groups to the model, as shown in Figure . The smallest model (57 atoms, model A) only includes copper, the cysteine Cys112 ligand, both histidine ligands His46 and His117, and the glycine Gly45 backbone. Model B (85 atoms) additionally includes the methionine ligand Met121 and an extension to the Gly45 backbone. Adding this backbone extension separately has been tested, but does not affect geometry or electronic structure. For model C, residues which form hydrogen bonds to Cys112 were added. The inclusion of Thr113 and Asn47 leads to a model with 124 atoms. In model D (154 atoms), hydrogen bonds from both histidine ligands were added, which leads to inclusion of the backbone of Gly9 to Gln12 and Met44 with a connecting H2O. In model E, the backbone connecting Cys112 and Met121 is added, leading to 206 atoms. To complete the second sphere around the Cu center, an additional methionine Met13 was added (model F, 232 atoms). In the final model (model G, 245 atoms), the so far omitted side chain of Phe114 is explicitly included in the QM part. Further increases in the QM size were not attempted, since these would result in inclusion of a large number of explicit water molecules beyond the protein surface. This would require conformational sampling to a point where size/functionality–property relations can no longer be unambiguously determined.

Figure 2

QM regions used in this study with labels indicating the residues included with respect to the previous model. The asterisk indicates that only the backbone of the residue is included.

Computational Details

All calculations were performed with the ORCA program package. QM/MM calculations used an interface of ORCA with NAMD.[55] DFT calculations used D3BJ dispersion corrections.[56] Relativistic effects were included with the use of the DKH2 Hamiltonian.[57−60] It is noted that the ZORA Hamiltonian[61−63] was also tested and no significant differences were observed from DKH2 in structural parameters. Appropriately recontracted DKH-def2-TZVP(-f) basis sets[64,65] were used for the DFT calculations, along with decontracted def2/J basis sets[66] for the RI approximation to the Coulomb integrals. The chain of spheres approximation (COSX)[67] was used for the exchange. Grids were increased to “Grid5” and “GridX5” in ORCA nomenclature. QM gas-phase optimizations were done using the same settings, with Cartesian constraints on the Cα atoms.

DLPNO-CCSD(T) energies were calculated using TightPNO thresholds on top of a TPSSh[68] reference. A DKH2 Hamiltonian was employed, with DKH-def2-TZVP basis sets on everything except copper, for which DKH-def2-TZVPP is used. The auxiliary basis sets were generated using the Autoaux feature.[69] DKH2 was employed with a finite nucleus using a Gaussian core model.[70] Picture change effects were accounted for.[71]

HFC parameters from DFT were calculated using the TPSSh functional with the DKH2 Hamiltonian with inclusion of picture change and finite nucleus effects.[70,71] The DKH-def2-TZVP(-f) basis set was used for all atoms except copper and sulfur. For Cu and S, the basis set employed was the DKH-def2-TZVP basis with fully decontracted s-functions and three additional s-functions created by scaling the tightest exponent of the original basis set by 15.625, 6.25, and 2.5.[72] No RI approximation was used in the calculation of spectroscopic properties. To calculate HFCs from DLPNO-CCSD densities, a previously established protocol was used[51] that utilizes very tight thresholds for the PNO generation and decontraction of the basis sets. In addition, the multifragment approach within the DLPNO framework is used,[73] with different thresholds for the PNO generation in the outer fragments. Here, thresholds vary from LoosePNO to TightPNO. Unrelaxed densities were used for the spectroscopic properties derived from the DLPNO-CCSD calculations.

Simulations of EPR and ENDOR spectra were performed using the Easyspin program.[74] Simulations of the EDNMR spectra were performed with a home-developed program by the Goldfarb lab that was already used for the interpretation of the experimental spectra.[27]

Results

Geometric Parameters

Since the crucial parameter to the Cu–SCys112 interaction is the distance between them, we require our computational models to yield a Cu–S distance close to experiment. In a first step, the experimental values for the copper ligand distances are analyzed (Table and Figure S1). The X-ray diffraction (XRD) structures show a large variety in the active-site distances. Possible reasons are the empirical restraints used in the refinement of the structure, differences in their resolution, or photoreduction during data collection. In the earliest crystal structure, the Cu–SCys112 distance was deduced to be 1.79 Å, which is unrealistically short. In subsequent XRD models, the distance varies between 2.12 and 2.27 Å. XRD has problems in showing the position of sulfur atoms close to electron-rich atoms, like copper. Here, EXAFS measurements helped to more reliably refine the distance to 2.12 Å.[35,36] The Cu–NHis distances vary similarly, between 1.94 and 2.11 Å for His117, and 1.93 and 2.15 Å for His46. However, there is an observed difference between the two Cu–NHis distances in each protein. The Cu–SMet121 distance is much longer than the Cu–SCys112 distance, but much better resolved. The Cu–O distance on the other hand varies between different structures.

Table 1

Distances (Å) Between Cu and Surrounding Atoms in the Reduced State of Azurin as Obtained from Structural Studies Reported in the Literature

	Cu–SCys112	Cu–NδHis117	Cu–NδHis46	Cu–SMet121	Cu–OGly45
1AZU[6]	1.79	2.42	2.15	3.21	2.47
4AZU[52]	2.27	2.11	1.99	3.18	2.84
	2.27	1.98	2.06	3.16	2.95
	2.24	1.96	2.13	3.21	3.05
	2.17	2.00	2.12	3.05	3.03
1DZ0[75]	2.16	2.02	2.03	3.26	2.75
1NWO[76]	2.13	1.94	2.01	3.01	3.16
	2.14	1.96	1.93	3.14	2.95
2CCW[77]	2.21	2.00	2.02	3.26	2.94
2AZA[78]	2.12	2.01	2.08	3.12	3.16
	2.17	1.99	2.09	3.10	3.09
XRD AVG	2.15	2.03	2.05	3.16	2.93
XRD STD	0.12	0.13	0.06	0.08	0.20
EXAFS[36]	2.12	1.86/1.94	1.86/1.94	3.39	2.82

Calculated distances depend on model used and the type of calculation: gas-phase DFT optimizations show the shortest Cu–SCys distances.[77] For QM/MM optimization, it depends on the type of coupling and embedding used.[31,77] The error of DFT was quantified in large-scale QM/MM calculations.[28] A summary of previous geometries obtained by computational modeling is given in Table S1. We initially tested a few commonly used density functionals for model B (for details, see Supporting Information Table S2). Based on the Cu–SCys distance as a main criterion, BP86, TPSS,[68] and TPSSh[79] were considered. The results are collected in Table . Additionally, we tested B3LYP,[80] because it was used in previous studies.[28] These functionals were tested against DLPNO-CCSD(T) energies (with DLPNO-CCSD(T) geometry optimizations are not feasible). The results show that TPSSh provides the structure with the lowest DLPNO-CCSD(T) energy, and hence this functional was used in further studies.

Table 2

Evaluation of Optimized Geometries Obtained with Different Functionals Showing the Cu–SCys112 Distance and the Relative DLPNO-CCSD(T) Single Point Energiesa

	Cu–S (Å)	ΔE (kcal/mol)
BP86	2.16	2.55
TPSS	2.16	0.81
B3LYP	2.18	0.61
TPSSh	2.16	0

The TPSSh geometry provides the lowest DLPNO-CCSD(T) energy.

In Table , the active-site distances computed with TPSSh for the model series are presented. The Cu–SCys distance decreases with increasing model size. Still, the EXAFS distance of 2.12 Å is not reproduced. Inclusion of dispersion effects in the form of D3 corrections[56] leads to slightly shorter distances. The latest D4 dispersion correction[81] was also tested, but did not change these results in any significant way. The influence of dispersion was additionally assessed by decomposition of the DFT energies. This indicated that the dispersion energy is about 1 order of magnitude smaller than the QM–MM interaction energy. Therefore, in geometry optimizations, the treatment of the environment by the QM/MM approach is much more significant than the inclusion of empirical dispersion corrections.

Table 3

Cu–Ligand Distances (in Å) from QM/MM Optimized Models Compared to Averaged EXAFS Distances

	Cu–SCys112	Cu–NδHis117	Cu–NδHis46	Cu–SMet121	Cu–OGly45
A	2.18	1.94	1.94		2.75
B	2.17	1.94	1.94	3.05	2.91
C	2.16	1.93	1.93	3.09	2.88
D	2.16	1.90	1.91	3.09	2.88
E	2.15	1.89	1.91	3.00	2.95
F	2.14	1.89	1.90	3.07	2.95
G	2.13	1.87	1.90	3.10	2.90
EXAFS[36]	2.12	1.90	1.90	3.39	2.82

Looking at the evolution of all Cu–ligand distances (Table ), it is found that also the Cu–NHis distances decrease with increasing QM size. The weakly bound ligands Met121 and the backbone of Gly45 however vary, depending on the steric hindrance that is included in the respective model.

While the distances show a limited picture, the analysis of angles (Table S3 and Figure S2) allows us to obtain a complementary view of the changes at the copper center upon increasing model sizes. While the Cu–SCys distance decreases, the Cα–Cβ–SCys angle, which defines the position of the side chain with respect to the backbone, remains constant. However, the Cβ–SCys–Cu angle, which was found crucial for description of the Cu–SCys interaction,[16,37] decreases with increasing model size. With increasing model sizes, also the Cu–NHis distances decrease. The corresponding angle between the His117 Nε–Nδ–Cu decreases, which describes a movement of the imidazole group of His117 toward the copper center. Such a movement is not observed for His46. Instead, the dihedral angle between the imidazole plane and the Cu–SCys vector decreases. This means that the imidazole group of His46 moves out of the plane defined by Cu, SCys, and Nε of His46. To compare this behavior to experimental geometries, the distances from the QM/MM optimized models are compared to the crystal structure averages in Figure S3. The Cu–SCys distances are within the range of distance reported in various crystal structures. However, for these first coordination shell distances, EXAFS is presumably providing more accurate numbers. Here, the EXAFS data give an average Cu–NHis distance of 1.90 Å.[36] The QM/MM results are in good to excellent agreement with this result. The QM–MM interaction was found to affect the Cu–SCys112 distance. Scans along the Cu–SCys coordinate indicate that neglect of the QM–MM interaction energy while maintaining the QM/MM structural constraints shifts the minimum by 0.05 Å toward longer distances. Interestingly, QM-only cluster optimizations that completely neglect the protein environment also lead to short Cu–SCys distances (Table S4) in agreement with previous gas-phase calculations,[77] but this is accompanied by other structural changes within the copper coordination sphere that are not consistent with the QM/MM geometries. Therefore, sufficient treatment of the environment in the QM/MM approach is essential for the description of the geometry of the copper site. These QM/MM structures were used as the basis for the calculation of spectroscopic parameters discussed in the following.

Hyperfine Coupling Constants

In the following, we present calculations of all hyperfine parameters that are relevant for understanding the electronic structure of the copper site and for establishing connections to experimental observations. We employ the highest level of theory available to us, DLPNO-CCSD. To evaluate the DLPNO-CCSD results, we present simulated spectra obtained with the calculated parameters and compare to experiment. This allows us to evaluate whether the calculated spin distribution over the active site, and in particular over copper and sulfur centers is accurate. Comparing the calculations with all experimental hyperfine information will allow us to develop a more global picture of the spin distribution and pinpoint where the calculations may fall short. The final outcome will be calibrated spin populations (in a given population analysis scheme) that may serve as reference for other theoretical methods. In the Supporting Information, we also provide a comparison with hyperfine couplings and spin populations from DFT.

The calculation of HFCs using DLPNO-CCSD was tested in a benchmark study for various molecules including transition-metal complexes by Saitow et al.[51] Due to the size of the azurin models, a multifragment approach was utilized.[73] This ensures to keep the accuracy of the established protocol (for a test on model A, see Table S5) but exceed its limits regarding the model size. In a two-fragment scheme, the atoms for which HFCs will be computed are included in fragment 1. For these atoms, the required tight thresholds are used.[51] The rest of each model is included in fragment 2. Using LoosePNO settings for fragment 2, the HFCs of all model sizes up to model G can be calculated. In Table , the results for the largest HFCs are shown.

Table 4

Model A to G DLPNO-CCSD Hyperfine Coupling Constants (MHz)a

	SCys112	NδHis117	NδHis46	HA	HB
A	19.6	23.8	21.3	21.6	9.0
B	19.1	24.6	20.8	17.8	13.7
C	19.4	24.9	21.8	17.3	14.8
D	15.6	31.6	21.3	10.6	12.5
E	14.2	32.7	22.1	9.5	10.5
F	14.8	32.4	22.2	10.3	11.0
G	15.2	32.8	22.1	10.4	11.9

Tight thresholds were applied in fragment 1 as defined in the text; LoosePNO thresholds were applied in fragment 2. A graphical representation is given in Figure S4. A similar table with the principal values is provided as Table S6.

This approach allows access to all sizes for the QM subsystem, but the drawback of the two-fragment scheme is the discontinuity between HFC atoms and the amino acid side chain that can lead to an incorrect description of spin delocalization. This can be observed, for example, in the case of HFCs of the cysteine β-hydrogens, but is also expected for the histidine ring. Therefore, we adopted instead a three-fragment scheme, where the amino acid side chain of each ligand is treated at an intermediate level of accuracy (fragment 2), and the rest of the model is treated as a lower-threshold fragment 3. This improves the quality of the ligand description (see Table ) albeit leading to increased computational cost, which means that not all model sizes can be treated this way. The best compromise between model size and accuracy was found for QM/MM model C. The corresponding Mulliken and Löwdin spin populations for model C are given in Table S7.

Table 5

Model B/C DLPNO-CCSD Hyperfine Coupling Constants (MHz)a

	SCys112	NδHis117	NδHis46	HA	HB
B (tight/loose)	19.9	24.8	20.7	22.3	17.4
C (tight/loose)	20.1	25.2	21.8	21.2	18.5

Tighter thresholds in fragments 1, 2, and 3 as indicated.

g-Tensor and Copper Hyperfine Coupling

Although the focus of our work is the ligand HFCs and their comparison to experiment, we briefly discuss the question of the g-tensor of the system and of the copper HFC. The g tensor is needed to define the orientation of the HFC tensors since the eigensystem of T  serves as the reference frame for the EPR simulations. As observed previously[32] and also shown in Figure S5, the calculated g tensor and copper HFC deviate from experiment. This is a known limitation of currently available theoretical approaches. Particularly, the copper HFC is one of the most challenging quantities to compute accurately, regardless of the level of quantum chemical approximation.[32,46,82] Therefore, to avoid any ambiguity, the experimental g values are used in the following, whereas the orientation is taken from the DFT calculations. It has generally been found to be in good agreement for the case of plastocyanin, where the g-tensor orientation had been deduced from single-crystal EPR experiments.[4,46] The orientation of the g tensor in the molecular frame is shown in Figure S6. The agreement of the calculated HFCs with experiment is shown by the comparison between simulated spectra using the calculated ligand HFCs and the experimental ENDOR[83] and EDNMR[27] experiments. By directly comparing simulated and measured spectra, we avoid any misinterpretation that may arise from HFC parameters deduced from simulating the experimental spectra.

Nitrogen Hyperfine Couplings

We consider four nitrogen atoms, the coordinating δ nitrogens of the histidines His117 and His46 and the remote ε nitrogens of the imidazole ligands. The remote nitrogen HFCs are of the order of 1 MHz,[21] which is perhaps too small to be analyzed with confidence given the intrinsic uncertainties of the computational method.[51] Nevertheless, the calculated values, 1.2 MHz for NεHis117 and 1.0 MHz for NεHis46, accurately reproduce the experimental observations, 1.3 MHz for NεHis117 and 0.9 MHz for NεHis46, as deduced from electron spin echo envelope modulation (ESEEM) experiments.[21]

The HFCs of the coordinating histidine δ-nitrogens are shown in Table . The NδHis117 HFCs agree very well with the ESEEM Aiso, but are smaller than the EDNMR results. Given the limited accuracy of the EDNMR results, the calculated NδHis117 HFCs show reasonable agreements. The NδHis46 HFC is about 1 MHz larger than the ESEEM value, but within the range of the EDNMR values. Again, for the EDNMR HFCs, a larger error is observed than for the ESEEM HFCs. This difference in the error might stem from the fact that the ESEEM was measured on a single crystal,[21] while the EDNMR data were obtained from frozen solution samples. A close-up on this issue is presented in the simulation of the experimental spectra below.

Table 6

Components of the NδHis HFCs Calculated with DLPNO-CCSD from Model C Compared to the Experimenta

NδHis117	Axx	Ayy	Azz	Aiso	NδHis46	Axx	Ayy	Azz	Aiso
DLPNO-CCSD	29.4	22.7	23.4	25.2	DLPNO-CCSD	25.7	19.6	20.2	21.8
ESEEM[21]	27.8 (±0.4)	24.0 (±0.3)	23.6 (±0.3)	25.1 (±0.3)	ESEEM[21]	19.1 (±0.3)	18.0 (±0.4)	17.2 (±0.4)	18.1 (±0.4)
EDNMR[27]	32.8 (±1.5)	25.0 (±1.5)	24.5 (±1.5)	27.4 (±1.5)	EDNMR[27]	24.0 (±0.8)	21.0 (±0.8)	17.8 (±0.8)	20.9 (±0.8)

Note that notation was adjusted to align with experiment, as A is the largest component of the HFC.

Along with the δ-nitrogen HFCs, the underlying interactions can be analyzed by decomposition of the HFC onto the multicenter components (Table ). Here, it is shown that the one-center terms, that is, the local contributions, are dominant among the other interactions. This indicates that the spin density can be approximately assigned to individual atoms. Hence, one can attempt to correlate a given HFC with an atomic spin population value. This is done here using the Mulliken population analysis scheme, which leads to DLPNO-CCSD spin populations of 0.0532 for NδHis117 and 0.0501 for NδHis46.

Table 7

Decomposition Analysis for the Nitrogen Hyperfine Coupling Constants

NHis117	Amin	Amid	Amax	NHis46	Amin	Amid	Amax
1-center	1.68	0.63	–2.31	1-center	0.98	0.82	–1.80
2-center pc	0.49	0.28	–0.77	2-center pc	0.29	0.38	–0.67
2-center bond	–0.25	–0.18	0.43	2-center bond	–0.14	–0.22	0.37
3-center	–0.03	0.01	0.02	3-center	–0.01	–0.01	0.02
total	1.89	0.73	–2.62	total	1.12	0.97	–2.08

The comparison with the ENDOR simulation is shown in Figure . In the experimental Rapid Passage ENDOR, the majority of the intensity lies in the signals corresponding to , while the transition is almost not detectable. The signals of His117 show good agreement with experiment, while the signals of His46 are slightly too high in energy. The calculated HFC for NδHis46 is 21.8 MHz, while experiment shows an HFC of 18.1 MHz. At the same time, the width of the signal deviates from experiment, which could arise from inhomogeneous broadening, or differences in the orientation of the calculated quadrupole interaction.

Figure 3

Comparison of the 14N ENDOR simulation with hyperfine and quadrupole coupling constants obtained from DLPNO-CCSD as described in the text (red line), compared to the simulation using parameters determined from experiment[83] (black dotted line) and the experimental ENDOR spectrum[84] (black line). Experimental parameters: ν = 35.2 GHz, B = 1113 mT.

To evaluate the agreement between the DLPNO-CCSD 14N HFCs and experiment in a more quantitative way, the DLPNO-CCSD HFCs are scaled in 5% increments. The results are shown in Figure . For the NδHis117 HFCs (top), the best agreement is observed if 5% of the value is added to the DLPNO-CCSD values. This corresponds to an HFC of 26.4 MHz, which is roughly midway between the experimental values obtained by Coremans et al.[21] and by the Goldfarb group.[27] Reducing the NδHis117 HFC on the other hand leads to a clear disagreement with experiment. For the NδHis46 HFC (Figure bottom), a more complicated picture is observed. The best agreement with experiment is observed if the DLPNO-CCSD Aiso is reduced by 20%. This corresponds to an HFC of 17.5 MHz, which is closer to the Coremans et al. value of 18.1 MHz.[21]

Figure 4

Scaling of the calculated DLPNO-CCSD HFCs for the Nδ of His117 (top) and His46 (bottom).

Given this difference between the DLPNO-CCSD HFC and the scaled HFC, one might wonder how this translates back to the spin population on the δ-nitrogens. The DLPNO-CCSD Mulliken spin population on the δ-nitrogen of His117 is 0.0532, and after scaling, it is 0.0558. On His46, the DLPNO-CCSD calculation yields a spin population of 0.0501, but after scaling the HFC down by 20%, only 0.0400 of the spin is found on the δ-nitrogen of His46. These numerical differences are not large, yet the simulations demonstrate that the spectra are sensitive to subtle differences of this magnitude. These differences can arise from slight variations in the geometries, such as a small rotation of the imidazole ring. Hence, high-accuracy calculations are needed to describe the system properly.

The second spectrum that can be used for comparison is the EDNMR of unsubstituted azurin. The simulation for the EDNMR spectrum at two different magnetic fields is given in Figure . Again, the subspectra of the individual histidines are shown. In both spectra, the general features below 30 MHz are reproduced very well. At 3048 mT, the calculated spectrum is narrower around the 20 MHz area than the experimental spectrum and a shoulder at 25 MHz is missing, which could be attributed to the differences between computed HFCs and experiment. The calculated NδHis46 HFCs of 21.8 MHz are slightly larger than the experimental results of 18.1 MHz[21] or the refitted value of 20.1 MHz.[27] While the DLPNO NδHis117 HFCs of 25.2 MHz agree very well with the ESEEM results,[21] they were refitted to 27.2 MHz in the EDNMR.[27] However, in the 40 MHz region, the calculated signals are visibly higher than the experimentally observed signals, i.e., both HFCs disagree. Since this signal is due to double quantum transitions, a twice large discrepancy is to be expected. Additionally, contributions from 63,65Cu lead to a broad signal around 40 MHz, which leads to difficulties in the interpretation. A similar picture is given by the simulated spectrum at 3302 mT, where the signals observed at around 45 MHz in the calculated spectrum cannot be found in the experiment. Again, a broad copper band situated between 30 and 50 MHz is masking the signals. To test the effect of the scaling, the results of the ENDOR comparison are used in EDNMR simulations in Figure . In the 3048 mT spectrum, the scaled nitrogen HFC of NδHis46 provides a too early onset in the 16 MHz region. In the −18 MHz region, however, the scaling improved the agreement with experiment. In the region beyond ±30 MHz, it is difficult to judge whether the scaling leads to an improvement. In the 3302 mT spectrum, the onset of the ±16 MHz signals has been improved by scaling the HFC values. Here again, the NδHis46 HFC is the main origin for the improvement. Similar to the 3048 mT spectrum, the region beyond 30 MHz cannot be judged due to the broad experimental lines.

Figure 5

Comparison of the simulated EDNMR spectra obtained from QM/MM calculations (red) with the subspectra of NδHis46 (blue) and NδHis117 (green) for 3048 mT (left) and 3302 mT (right) to the experiment[27] (black lines). Microwave frequency ν = 94.9 GHz.

Figure 6

Comparison of the simulated EDNMR spectra obtained from QM/MM calculations (red) with the subspectra of NδHis46 (blue) and NδHis117 (green) for 3048 mT (left) and 3302 mT (right) to the experiment[27] (black lines). Dotted lines indicate scaled Nitrogen HFCs according to the ENDOR discussion before. Microwave frequency ν = 94.9 GHz.

Proton Hyperfine Couplings

The proton HFCs can provide a complementary perspective on the electronic structure of the copper site. We first look at the HFCs of imidazole ring protons of histidines His46 and His117, which are of the order of 1 MHz.[22] The DLPNO-CCSD results are given in Table . Although these values are small, there is convincing correspondence of the computed values with those deduced from NMR experiments.[22]

Table 8

DLPNO-CCSD Hyperfine Coupling Constants (in MHz) of the Histidine Ring Protons, Compared to Experiment

	NεH(His117)	CεH(His117)	NεH(His46)	CεH(His46)
DLPNO-CCSD	1.0	0.8	1.0	0.9
NMR[22]	1.1/1.5	0.56	1.1/1.5	0.56

The situation with the two cysteine β-protons, experimentally termed H1 and H2, is more complicated. In contrast to the histidine ring protons, selective isotope substitution is not possible for the protons on the cysteine β-carbon. Hence, while it is possible to determine individual HFC values for each β-proton, no clear assignment to the individual atoms was possible.[22−24] However, orientation-dependent 1H ENDOR measurements provide information on the orientation of the A tensor of each cysteine β-proton.[25] Using this information, a direct assignment with the calculated proton HFCs is possible, by comparing the calculated orientations to the HFC orientation determined experimentally (cf. Figure ).[25] From this analysis, it can be deduced that DLPNO-CCSD HA corresponds to 1H-ENDOR H2,[25] and vice versa (HB to H1). Looking at the absolute values of the HFCs for the cysteine β-protons (Table ), this assignment seems opposite to what is suggested numerically, but given the possible error of the DLPNO-CCSD method and experimental uncertainties, the directions of the principal axes of the hyperfine tensor are probably more reliable in this case, and therefore we suggest that the correspondence of nuclei established from the tensor orientations is more reliable as well.

Figure 7

Orientation of the Cys112 β-proton HFC tensors at the azurin T1 center. Directions: z in blue, y in red, and x in green.

Table 9

DLPNO-CCSD Hyperfine Coupling Constants (in MHz) of the Cysteine β-Protons, Compared to Experiment

	Axx	Ayy	Azz	Aiso
DLPNO-CCSD HA	18.7	19.8	25.0	21.2
DLPNO-CCSD HB	15.9	17.1	22.6	18.5
ENDOR[25] H1	20.4	21.3	26.2	22.6
ENDOR[25] H2	14.4	19.1	23.0	18.8

The proton ENDOR spectrum can similarly be simulated and compared to experiment (Figure ). In general, a very good agreement is observed between experiment and the spectrum obtained on the basis of DLPNO-CCSD computed values. In the g direction, perfect agreement is observed for the HA subspectrum, although it is not as intense as experiment. HB agrees very well with experiment, especially considering the simulation of the experimental parameters. In the g direction, again good agreement is observed. Here, the signals of the DLPNO-CCSD simulation are clearly separated, while the experiment is convoluted in a single peak. In the g direction, the opposite is observed. Here, the DLPNO-CCSD parameters of the individual protons overlap, while the experiment shows a wider signal with two individual peaks. These slight disagreements might originate from differences in the orientation of the tensor, in combination with the position of the β-protons in the model.

Figure 8

Comparison of the simulated 1H ENDOR using the DLPNO-CCSD parameters (red line), consisting of the HA (orange) and HB (gray) subspectra with experiment[25] (black full lines) and simulation of the experimental parameters (black dashed lines), for the individual directions of the g tensor. Field of 3328.69 mT (g, νZeeman = 141.73), 3300.37 mT (g, νZeeman = 140.52), and 3016.98 mT (g, νZeeman = 128.46). Microwave frequency: ν = 94.9 GHz.

At this point, the origin of the HFC can be discussed. As shown from the decomposition of the HFC (Table ) in the multicenter parts, the majority originates from the two-center nonbonding interaction, which translates to the dipolar terms. Hence, a translation back to the spin population on the β-hydrogens as done for the δ-nitrogens does not make sense. These small differences clearly show that accurate methods are needed to determine the electronic structure and also that small differences in the geometry could dramatically change the results in this case.

Table 10

Decomposition Analysis for the Cys112 β-Proton Hyperfine Coupling Constants

HA	Amin	Amid	Amax	HB	Amin	Amid	Amax
1-center	0.03	–0.06	0.04	1-center	0.04	–0.06	0.02
2-center pc	1.64	1.98	–3.62	2-center pc	2.04	2.27	–4.32
2-center bond	0.37	–0.50	0.13	2-center bond	0.41	–0.52	0.10
3-center	–0.63	–0.51	1.15	3-center	–0.64	–0.59	1.22
total	1.40	0.91	–2.31	total	1.86	1.11	–2.96

33S Hyperfine Couplings

A crucial spectroscopic quantity is the sulfur HFC, which has been limited by the applicability of 33S labeling to the Cys112 ligand. Recently, a 1D-EDNMR spectrum was recorded, where the sulfur HFC could be determined.[27] Unfortunately, the resolution of these 1D-EDNMR spectra is rather limited; hence, 2D-EDNMR spectra were also recorded. In the following, first, we compare the computed HFC values, and then we proceed with simulations of the EDNMR spectra. Unfortunately, the two experimental approaches yield slightly different fitted parameters for the sulfur HFC (Table ), both for the isotopic value and for the principal values. While the general trends are the same, a larger A value but smaller A value is observed for the 1D-EDNMR. The calculated 33S HFCs show good agreement with the 1D-EDNMR experiment, but yield an even larger Aiso. Hence, the agreement with the fits from the 2D-EDNMR is rather limited.

Table 11

33S Principal HFC Values (MHz) from DLPNO-CCSD, Compared to Experiment

sulfur	Axx	Ayy	Azz	Aiso
DLPNO-CCSD	–11.5	–13.8	85.5	20.1
1D-EDNMR[27]	–15.4 ± 0.6	–17.0 ± 0.6	89.0 ± 0.6	18.8 ± 0.6
2D-EDNMR[27]	–15.4 ± 1	–27.0 ± 1	67.5 ± 2.5	8.3 ± 1.5

The sulfur HFCs are compared to the EDNMR experiment with 33S-labeled azurin.[27] An overview is given in Figure . The 14N and 33S signals overlap strongly. At 3048 mT, the sulfur transitions are observed as resolved features in 20–30 MHz, next to the nitrogen transitions at 20–25 mT other clear features appear at 40 and 60 MHz. At higher fields, shown in the case of 3302 MHz, the sulfur signals are broader, less resolved, and span the area between 10 and 50 MHz. If the subspectra are summed, a good agreement in the shape of the features is observed. However, the detailed signature of the 33S is lacking from experiment, especially at higher energies. At 3302 mT, the calculated spectrum rises below 20 MHz, while the experimental onset is found at slightly higher energies. The interpretation of the experimental spectra, and consequently the comparison to the calculation, is limited due to the complexity of the experiment and the large span of the 33S EDNMR signals of up to 100 MHz. Because of this, a unique assignment of the calculated parameters to the experiment is not straightforward.

Figure 9

Comparison of the simulated 1D-EDNMR spectra obtained from DLPNO-CCSD calculations (red, broadened to better approximate the experiment) with the subspectra of SCys112 (orange), NδHis46 (blue), and NδHis117 (green) for 3048 mT (left) and 3302 mT (right) to the experiment[27] (black dashes). Microwave frequency: ν = 94.9 GHz.

Similarly to the nitrogen HFCs, the sulfur HFCs were scaled and compared to the experiment. As shown in Figure , the comparison is not straightforward due to the limited resolution of the experiment. However, looking more closely at the signals around 20 MHz, it becomes obvious that values higher than the calculated sulfur HFCs are not in agreement with the experiment. This is particularly interesting since the original fits show lower HFC values for the sulfur HFC.

Figure 10

Scaling of the calculated DLPNO-CCSD HFCs for 33S at a field of 3048 mT (top) and 3302 mT (bottom) compared to the experiment[27] (black lines). The sum of the nitrogen spectra (dashed lines) is given as a guide. Microwave frequency: ν = 94.9 GHz.

To overcome the resolution issues in the 1D-EDNMR, 2D-EDNMR spectra were recorded.[27] The simulation of the 2D-EDNMR spectra using the DLPNO-CCSD computed parameters was attempted for completeness. The results are discussed in the Supporting Information. As shown in Figure S7, the computed values are in overall agreement with experiment; however, the complexity of the spectra and the uncertainties involved in defining their information content do not permit a detailed analysis in this case.

It is interesting to examine how the above translates to the spin population on the sulfur. The Mulliken spin population obtained from the DLPNO-CCSD calculation is 0.37. The unscaled 33S HFCs show a reasonable agreement with experiment, and all scaling tests show that a higher 33S HFC leads to better agreement with experiment, which would correlate to a higher spin population. This would then lead to an increased covalency between the copper and sulfur.

Additional insight is offered by the decomposition of the DLPNO-CCSD HFCs into multicenter terms. The one-center term describes the dipole interaction between the nucleus and the spin at atom A. The two-center terms are subdivided into bonded and nonbonded interactions. They describe the dipole interaction between the nucleus of atom A and spin at atoms B and C, bonded and nonbonded respectively. Similarly, the three-center term describes the interaction of the nucleus of atom A and the spin at atom C. The results for the HFC decomposition analysis, presented in Table , indicate that the one-center terms are indeed dominating the results for sulfur. However, for a quantitative analysis, the two-center bonded contributions should definitely be included, as they make up about 10% of the final HFC value.

Table 12

DLPNO-CCSD HFC Decomposition for 33S HFCs (in MHz)

	Aiso	Amin	Amid	Amax
1-center	21.94	28.97	13.49	–42.46
2-center pc	–0.17	–0.70	0.28	0.42
2-center bond	–1.82	0.01	–0.09	0.08
3-center	0.12	0.04	–0.01	–0.03
total	20.07	28.32	13.66	–41.98

Overall, the above analysis indicates that while the one-center terms are of general importance for the sulfur nucleus, it is not possible to neglect the multicenter terms in the analysis of 33S HFCs. Thus, the spin population should be discussed by comparing the experimental HFCs to sufficiently accurate quantum chemical calculations.

Discussion and Conclusions

In this work, the geometry and electronic structure of the azurin T1 copper center were assessed by a series of computational models. The copper environment was studied with a QM/MM model series of increasing size, to obtain a reasonable copper–sulfur bond distance. The hyperfine coupling constants of 13N, 1H, and 33S nuclei of the copper ligands were computed using the highest level of theory currently available to us, namely, the DLPNO implementation of CCSD. The results and the accompanying analysis demonstrate the need for high numerical accuracy. By simulating the respective spectra using the DLPNO-CCSD parameters, a direct comparison between calculation and experiment is possible. In contrast to the performance of density functional theory (DFT, results presented in the Supporting Information), the calculated DLPNO-CCSD parameters provide meaningful insights into the electronic structure of the azurin Cu(II) center and yield results in good agreement with experiment, or at least within the experimental uncertainty, even though in certain cases, the width and complexity of experimental spectra do not allow highly accurate comparisons and conclusive deductions. Here, spectroscopic studies on well-defined copper model complexes might be of aid to improve the interplay between calculation and spectroscopy. Importantly, we showed that the DLPNO-CCSD method is able to describe the highly covalent Cu–SCys interaction in the azurin copper center very well by comparing the calculated HFCs directly with experiment instead of relying on empirical relationships.

An important question that has been a central point of contention in the literature on the azurin Cu site is the concept of spin population and the covalency of the Cu–SCys bond. We note that often the term “spin density” is inappropriately used in experimental works to refer to what is correctly termed “spin population”, i.e., the more or less arbitrary partitioning of the continuous spin density (a global physical quantity) and its assignment to individual nuclei or atoms. From a quantum chemical perspective, there is no unique way of performing this partitioning. Hence, the numerical values depend on the protocol used, for example, the Mulliken or Löwdin population analysis scheme. In fact, from the numbers collected in the Supporting Information (Table S7), the Mulliken and Löwdin spin populations on copper or sulfur deduced from the same many-particle wavefunction differ by up to 6%. This is a strong reminder that these spin populations do not represent physical observables, and if one wants to judge the quality of a given calculation, one should resort to actual physical observables rather than derived quantities like spin populations. Indeed, a large variety of different values for the spin populations in type 1 copper sites have been reported in the literature on the basis of different flavors of DFT (Table S9). If we had to discuss the Cu–SCys covalency in these terms, our best estimates deduced from fitting the DLPNO-CCSD values to experimental spectra suggest that the average of the two spin population schemes (Table S7) on copper and sulfur amount to 55.3 ± 2.9 and 34.3 ± 2.9%, respectively. For sulfur, this is larger than the 30% suggested by EDNMR,[27] but possibly slightly smaller than the 37% obtained from XAS.[18,19] Note that either technique relies on a number of strongly implying assumptions, e.g., one-center approximation in the case of magnetic resonance techniques or reference compound calibrations in the case of XAS. The range of validity of these assumptions has been discussed in the present work and shown to be rather limited. More in-depth discussions can be found in the literature for XAS[85−87] and magnetic resonance spectroscopies.[33,34,88]

Despite the fact that the present study perhaps represents the theoretically most rigorous attempt to compute metalloprotein active site spectra, there is still ample room for improvement on the computational side, since neither the structural models that are used to represent the active site nor the theoretical method itself are perfect. For example, the DLPNO-CCSD framework for the calculation of hyperfine coupling constants currently does not include triple excitations, nor can it directly address at the moment the second-order spin–orbit contribution to the hyperfine interaction.[51] Yet, we have shown that it is possible to obtain a very accurate electronic structure picture of all ligand atoms surrounding the azurin copper site. We have also demonstrated what is in our view the best possible use of the raw computed values, that is, the simulation of the spectra and the direct comparison with experimental spectra. Therefore, the DLPNO-CCSD method can be of great assistance in the interpretation of highly complex experiments, and we hope that the present study served to demonstrate the principles of such an approach. Further refinements in this line of research are underway, and we expect that they will have wider implications in raising the standards of what is considered as “mainstream” quantum chemical approaches for the EPR spectroscopy of metalloenzymes in the near future.