Decoding the Ambiguous Electron Paramagnetic Resonance Signals in the Lytic Polysaccharide Monooxygenase from Photorhabdus luminescens.

Understanding the structure and function of lytic polysaccharide monooxygenases (LPMOs), copper enzymes that degrade recalcitrant polysaccharides, requires the reliable atomistic interpretation of electron paramagnetic resonance (EPR) data on the Cu(II) active site. Among various LPMO families, the chitin-active PlAA10 shows an intriguing phenomenology with distinct EPR signals, a major rhombic and a minor axial signal. Here, we combine experimental and computational investigations to uncover the structural identity of these signals. X-band EPR spectra recorded at different pH values demonstrate pH-dependent population inversion: the major rhombic signal at pH 6.5 becomes minor at pH 8.5, where the axial signal dominates. This suggests that a protonation change is involved in the interconversion. Precise structural interpretations are pursued with quantum chemical calculations. Given that accurate calculations of Cu g-tensors remain challenging for quantum chemistry, we first address this problem via a thorough calibration study. This enables us to define a density functional that achieves accurate and reliable prediction of g-tensors, giving confidence in our evaluation of PlAA10 LPMO models. Large models were considered that include all parts of the protein matrix surrounding the Cu site, along with the characteristic second-sphere features of PlAA10. The results uniquely identify the rhombic signal with a five-coordinate Cu ion bearing two water molecules in addition to three N-donor ligands. The axial signal is attributed to a four-coordinate Cu ion where only one of the waters remains bound, as hydroxy. Alternatives that involve decoordination of the histidine brace amino group are unlikely based on energetics and spectroscopy. These results provide a reliable spectroscopy-consistent view on the plasticity of the resting state in PlAA10 LPMO as a foundation for further elucidating structure-property relationships and the formation of catalytically competent species. Our strategy is generally applicable to the study of EPR parameters of mononuclear copper-containing metalloenzymes.

Introduction

Lytic polysaccharide monooxygenase (LPMO) enzymes are intensely studied for their ability to oxidatively cleave the glycosidic bond of recalcitrant polysaccharides such as cellulose[1−3] and chitin.[4,5] Given the societal impulse to find new sustainable alternatives to meet global energy demands, there is great interest to apply the LPMO reactivity in the production of second-generation biofuels from biomass.[6,7] Understanding the structure–function properties of the LPMO active site will give great insights into the development of future biomimetic models.[8] The active site of LPMOs is composed of a mononuclear copper center coordinated by the N-terminal histidine amino acid in a bidentate fashion and a second histidine side chain to form a motif known as histidine brace. Two water molecules complete the first coordination sphere of the Cu(II) metal center (Figure ).[1,4] This coordination motif is very rare and is only observed in a few copper-containing metalloproteins such as the PmoB subunit of another enigmatic enzyme, the particulate methane monooxygenase,[9] CopC[10,11] and PmoF,[12] two bacterial copper-transport proteins, and LPMO-like fungal copper-transport proteins.[13,14]

Figure 1

Structure of the PlAA10 LPMO enzyme (PDB code: 6T5Z)[15] and of its active site, with special focus on the copper center and its first coordination sphere. Bond distances shown in Å.

In the enzyme resting state, the copper center is at the +II oxidation state (3d9 configuration), leading to an S = 1/2 spin state.[1] To investigate copper-containing systems, electron paramagnetic resonance (EPR) spectroscopy is the technique of choice as it enables probing the electronic and geometric configuration of the paramagnetic metal center with the information encoded in the g and hyperfine coupling, A, tensors.[16,17] Despite the presence of a well conserved Cu site in the various LPMO enzymes, their EPR parameters are quite diverse,[6,18] often complicating the reliable structural interpretation of EPR signals. In the chitin-active LPMO from Photorhabdus luminescens, named PlAA10, and produced in our group, two co-existing species were detected by X-band EPR spectroscopy at 120 K using a frozen solution of PlAA10 in MES [2-(N-morpholino)ethanesulfonic acid] buffer at pH 6.5.[15] The first major species displayed the expected rhombic EPR signature (g ≠ g ≠ g, ACu ≠ ACu ≠ ACu), in agreement with the EPR signature of other chitin-active LPMOs.[5,19−21] This rhombicity is indicative of a distorted geometry falling between square base pyramidal and trigonal bipyramidal. On the contrary, the second, minor, species had distinct axial features (g > g ≈ g, ACu > ACu ≈ ACu). By EPR simulations, the major/minor species distribution was fitted to 80% rhombic and 20% axial.

Another chitin-active LPMO, SliLPMO10E from Streptomyces lividans, showed similar EPR characteristics.[21] The copper bound protein in a Tris [tris(hydroxymethyl)aminomethane] buffer at pH 7.0 containing NaCl exhibited two distinct EPR signals at 10 K: one of rhombic character and another of axial character. Chaplin et al. suggested that the origin of the minor species with axial character in the SliLPMO10E mixture may be due to the decoordination of the primary amine of the N-terminal histidine.[21] In addition, the mixed signal may originate not only from this decoordination alone but also from a possible coordination of chloride anions present in excess (150 mM NaCl) in the Tris buffer (Scheme ). It was also suggested that because the heterogeneity was observed in the apo-protein, the two configurations of the N-terminal coordination arm reflect a natural flexibility of LPMO enzymes for copper loading. The same rationale, that is, decoordination of the primary amine of the histidine brace, was suggested for the formation of the minor species in the PlAA10 enzyme mixture. This heterogeneous solution may form during freezing of the sample in the absence of a glassing agent, for example glycerol, before the EPR measurement.[15,18]

Scheme 1

Proposed Copper Coordination in SliLPMO10E (left) and in BlAA10[21−23]

Recently, Lindley et al. reported the study of a chitin-active LPMO, BlAA10 from Bacillus licheniformis, that was investigated in different pH conditions using EPR spectroscopic techniques and computations. Their work focused on the protonation states of the coordination sphere of the protein active site. They proposed the species observed at moderate basic pH to result from the deprotonation of the ligating water molecules and the partial decoordination of the −NH2 group from the histidine brace motif.[22]

Finally, Serra et al. investigated a new enzyme that belongs to the AA10 family, PpAA10 from Pseudomonas putida, by the means of Fourier transform infrared spectroscopy and EPR spectroscopies. Probing various buffer conditions, they investigated the influence of the pH on the EPR signature of the active and observed the presence of two different contributions to the spectra, having a rhombic and axial character, respectively. For the latter, they suggest that its formation would likely result from either alteration of the ligand water being replaced by a chloride ligand from the buffer or deprotonated, or coordination of an extra residue to the copper center (Scheme ).[23]

In light of these recent works, we seek to establish a robust methodology to elucidate the origin of the intriguing EPR signal distribution in the AA10 family using the PlAA10 LPMO enzyme by building theoretical model systems of the active site and calculating their EPR parameters using benchmarked protocols based on density functional theory (DFT) methods. These models are constructed using the PlAA10 crystal structure and considering numerous subsequent modifications around the copper center, such as decoordination of water molecules and of the primary amine, along with changes in protonation states of all relevant groups.

The pre-requisite of this approach is to be able to predict accurate EPR parameters from a theoretical model system. Therefore, the quantum chemical methods employed must be of proven reliability for both the g and A tensors.[24] Our group has recently presented a highly accurate protocol for the calculation of Cu(II) hyperfine coupling constants based on a specific density functional and a customized flexible-core basis set.[25] No method of similarly high accuracy has been established for the calculation of Cu(II) g-tensors;[26] therefore, in the present study, we first undertake a thorough methodological evaluation using a large set of reference synthetic complexes. On this basis, we define a new modified functional that yields the most accurate g-tensor values compared to any other method reported in the literature, either DFT- or wave-function-based, and hence will serve as a reference method for any future study of Cu(II) systems. Applying our benchmarked protocols on various atomistic models of the PlAA10 enzyme, we elucidate the chemical nature of the species resulting in the axial and rhombic EPR signals by comparing our computational results to the EPR measurements performed under different pH conditions. This work showcases the success of combining EPR spectroscopy with computational chemistry and contributes to a precise description of biologically relevant species in the family of AA10 LPMO enzymes.

Methodology

Experimental Details

Cu(II)-loaded PlAA10 was produced as previously described.[15] Protein samples (concentration around 200 μM) for continuous-wave EPR were prepared in either 50 mM MES buffer (pH 6.5) or 50 mM Tris·H2SO4 buffer (pH 8.5). EPR spectra were recorded on a Bruker Elexsys E500 spectrometer (Bruker, Karlsruhe, Germany) operating at X-band at 120 K (BVT 3000 digital temperature controller) with the following acquisition parameters: modulation frequency, 100 kHz; modulation amplitude, 5 G; conversion time, 90 ms; sweep time, 92.1 s; and microwave power, 20 mW. EPR spectra were simulated using the EasySpin toolbox developed for MATLAB.[27] The optimum Hamiltonian parameters have been obtained using the second order perturbation, and then, an exact diagonalization has been used for the final simulations. Pseudo-modulation treatment of the spectra was performed to graphically extract the exact number of nitrogen centers contributing to the superhyperfine pattern.[28,29]AN constants were considered isotropic and used in the final simulation of the first-derivative spectra. The Hamiltonian used for the simulations is the following equationwith i = the number of nitrogen centers.

Computational Details

Structural Models

For the methodological benchmarking part of this work, the coordinates of the 20 studied Cu complexes (Figure S1) were obtained from the Cambridge Structural Database (CSD) and were edited for completeness and chemical accuracy, for example, by removing solvent molecules or non-coordinating counterions and adding hydrogen atoms. The curated structures were subsequently employed in geometry optimizations. The set of compounds has been recently used in benchmarking quantum chemical methods for hyperfine coupling constants.[25]

Structural models of the PlAA10 active site were constructed from the 6T5Z X-ray crystal structure with 1.6 Å resolution.[15] The models include the amino acids that directly coordinate the Cu ion, His26 and His115, and second-sphere residues Ser59, Gly27, Ile113, Phe188, Trp179, and Ile181. Distinct interactions between the Cu coordination sphere and the protein envelope were conserved by also including fragments of the amino acids Tyr28, Gln58, Leu60, Met112, Gln114, Lys116, Thr117, Thr180, and Ala187. Two copper-coordinated crystallographic water molecules were included. Outer-layer crystallographically resolved solvent molecules were not included directly in the DFT models; long-range effects were accounted for through implicit solvation. The resulting model, denoted A1 in the present work (Figure ) has a five-coordinated Cu ion with two water ligands and consists of 197 atoms. All other structural modifications that were considered and evaluated in the present study were derived from this initial model.

Figure 2

DFT model A1 with the labels of the amino acids as defined in the 6T5Z X-ray structure of PlAA10 LPMO.[15] Hydrogen atoms bonded to carbon atoms are omitted for clarity. Truncated peripheral residues are denoted in square brackets.

Geometry Optimizations

All calculations were performed with Orca quantum chemistry software.[30,31] Geometry optimizations were carried out using the BP86 functional[32,33] with the all-electron def2-TZVP basis sets.[34] For the Coulomb fitting, the def2/J auxiliary basis sets were used.[35] To investigate the dependence of the optimized geometry on the type of functional, geometry optimization of model A1 was also performed with the hybrid functional B3LYP.[36,37] Key structural parameters obtained with the BP86 and B3LYP functionals are compared in Table S1. Differences in Cu–ligand bond lengths of less than 0.03 Å and in N2–Cu–O2 and N1–Cu–N3 angles of less than 0.4° are observed, indicating negligible sensitivity of the calculated structural parameters on the choice of functional. The conductor-like polarizable continuum model with ε = 20, typically used to study enzymatic systems with active sites at the surface in buffer solutions,[38] and the default refractive index of 1.33 was used in geometry optimizations of the enzyme models. Increased angular and radial integration grids (Grid4 and IntAcc 6.0, respectively, in Orca convention) and tight SCF convergence criteria were employed. Slow SCF convergence settings were applied in combination with recalculation of the full Fock matrix for each SCF step (DirectResetFreq 1) because this approach was found to ensure convergence even for uncommonly—electronically and computationally—difficult cases. In geometry optimizations of the PlAA10 LPMO models, constraints were applied on selected backbone carbons and terminal hydrogen atoms (details provided in the Supporting Information) in order to maintain the structural effect of the protein matrix on the active site geometry.

Calculations of EPR Parameters

Calculations of the copper hyperfine coupling constants and g-tensors employed increased general and radial integration grids (Grid6 and IntAcc 6.0), and specially enhanced grids for the copper center (SpecialGridIntAcc 11). The spin–orbit coupling operator was treated by an accurate mean-field (SOMF) approximation to the Breit–Pauli operator (SOCType 3).[39,40] The potential was constructed to include one-electron terms, compute the Coulomb term in a semi-numeric way, incorporate exchange via one-center exact integrals including the spin–other orbit interaction, and include local DFT correlation (SOCFlags 1,2,3,1).

All EPR calculations used the previously defined aug-cc-pVTZ-Jmod basis set[25] for Cu, while the def2-TZVP basis sets[34] were used for all other atoms. The Cu basis set, a modified version of the property-optimized basis set originally proposed by Hedegård et al.,[41,42] was shown to yield converged results for DFT calculations of hyperfine coupling constants.[25] In contrast to the hyperfine coupling tensor that require a highly flexible and accurate description of the core region, the g-tensor is less sensitive to the choice of basis set and shows faster convergence, provided that the valence region is sufficiently well described. Even though a smaller and more standard basis set on Cu might be sufficient for g-tensor calculations, the large aug-cc-pVTZ-Jmod, which is flexible and accurate also in the valence space, is employed here for g-tensors as well because the computational overhead is acceptable, the results are converged (Table S2), and for maintaining consistency over all EPR property calculations.

A-tensors were computed with the B3PW91 functional,[36,43] which was shown to be the best functional for copper hyperfine coupling constants in our recent evaluation study.[25] An essential element of most modern DFT approximations is the admixture to the density functional of an adjustable amount of the electron exchange contribution from the Hartree–Fock theory, also known as “exact exchange” (the term exchange referring here to the quantum mechanical energy component arising from the antisymmetry requirement of the electronic wave function). For the calculation of g-tensors, we defined a modified version of B3PW91 with 40% exact (Hartree–Fock) exchange, following an extensive benchmarking of various DFT methods that is described in the first part of the present work. Other DFT functionals evaluated for the calculation of g-tensors include the standard generalized gradient approximation (GGA) functional PBE,[44] the meta-GGA TPSS,[45] the hybrids B3LYP,[36,37] PBE0,[46] and BHandHLYP, the hybrid-meta-GGA TPSSh,[47] the long-range corrected or range separated functionals LC-BLYP[48] and CAM-B3LYP,[49] and the double-hybrid functionals B2PLYP,[50] DSD-PBEP86,[51] and PBE-QIDH.[52] The RI approximation was used for the MP2 part in combination with a very large automatically generated correlation fitting auxiliary basis sets.[53] The NoFrozenCore option and relaxed MP2 densities were used in the double-hybrid calculations. The exact exchange admixture was adjusted using the keywords ScalHFX and ScalDFX in Orca. Whenever we increased the contribution of exact exchange (ScalHFX), we decreased the default DFT exchange percentage by exactly the same amount. The use of scalar relativistic Hamiltonians, specifically the second-order Douglas–Kroll–Hess (DKH2)[54−60] and the zeroth-order regular approximation,[61−63] does not improve on average the quality of DFT calculations of A-tensors for copper systems.[25] In the present work, we tested both Hamiltonians also for the calculation of g-tensors for a subset of the complexes (Table S3), applying picture change effects and recontracted versions[64] of the def2 basis sets for the ligands. Similarly to the hyperfine coupling constants,[25] an insignificant and not clearly beneficial effect was observed with the use of scalar relativistic Hamiltonians in the calculation of g-tensors; therefore, this approach was not adopted in the evaluation of the enzyme models.

Results and Discussion

DFT Methodology for Accurate Cu(II) g-Tensor Calculations

DFT has proven successful in the calculation of accurate hyperfine coupling constants for mononuclear copper systems. Previous benchmarking studies by Sciortino et al.[26] and by us[25] have identified density functionals that perform consistently well and with predictive accuracy. Moreover, the critical role of the basis set used for copper was investigated in detail and an optimal methodology has been proposed that combines the B3PW91 functional with the aug-cc-pVTZ-Jmod basis set.[25] This approach was shown to consistently outperform applicable wave function methods for A-tensors.[25] However, calculations of g-tensors for Cu systems present a more complex challenge.[65−67] The evaluation study by Sciortino et al.[26] documented the performance of representative functionals for g-tensors of several copper complexes, but the size and spread of errors suggest that no functional is a clearly and systematically superior choice. Remarkably, wave-function-based methods perform in general more poorly than DFT for Cu g-tensors,[68−71] which underlines the difficulty of defining a robust and generally applicable theoretical protocol. In view of this problem, and given the need to have a reliable method for evaluating the many structurally different LPMO models, we begin this study by an extensive methodological benchmarking on a large reference set of spectroscopically characterized mononuclear Cu complexes. As described in the following, this enabled us to define a superior theoretical approach that gives strong confidence in the EPR parameters predicted for the LPMO models.

We employed the set of 20 mononuclear Cu(II) complexes depicted in Figure S1, as used also for the benchmarking of hyperfine coupling constants.[25] These complexes have experimentally resolved g-tensors obtained from EPR measurements (Table ). The set exhibits great variability in the properties of ligands: the coordinating atoms are of three different types (N, O, and S) and the coordination to the Cu center is in a 4N, 4S, 4O, 5N, or 6S fashion, or with a combination of atom types: 2N2O, 2N2S, or 3N1O. Some ligands are chelating, ranging from bidentate to tetradentate, and of varying sizes. The coordination geometry is also changing and includes square planar, tetrahedral, distorted square base pyramidal, and octahedral. The complexes also display a wide range of gmax values, between 2.085 and 2.285. We note that in this work, we adopt the convention of equating gmax with g in the case of computed values, but we maintain both nomenclatures for greater clarity in specific settings. The variety in coordination and properties provides us with a broad range to evaluate the performance of different computational methods.

Table 1

Experimental g-Values of the Complexes Considered in the Calibration Studya

	complexb	gx	gy	gz	refs.
1	[Cu(NH3)4]2+	2.047	2.047	2.241	(72)
2	[Cu(dtc)2]	2.025	2.025	2.085	(73)
3	[Cu(acac)2]	2.060	2.060	2.285	(74)
4	[Cu(en)2]2+	2.040	2.046	2.202	(75)
5	[Cu(mnt)2]2–	2.023	2.023	2.093	(76)
6	[Cu(gly)2]	2.052	2.052	2.267	(72)
7	[Cu(kts)]	2.030	2.030	2.140	(77)
8	[Cu(sac)2]	2.050	2.050	2.240	(78)
9	[Cu(im)4]2+	2.047	2.047	2.262	(72)
10	[Cu(py)4]2+	2.053	2.053	2.263	(72)
11	[Cu(eta)]2+	2.030	2.030	2.160	(79)
12	[Cu(epa)]2+	2.053	2.053	2.213	(80)
13	[Cu(atpt)]2+	2.100	2.100	2.235	(81)
14	[Cu(GGH)]−			2.173	(26)
15	[Cu(GGG)]−			2.202	(26)
16	[Cu(salpn)]	2.060	2.060	2.261	(82)
17	[Cu(S,S-mnpala)]	2.060	2.060	2.240	(83)
18	[Cu(salen)]	2.052	2.052	2.192	(84)
19	[Cu(bipy)2(NCS)]+	2.088	2.088	2.259	(85)
20	[Cu(ttcn)2]2+	2.045	2.045	2.117	(86)

Complexes 14 and 15 were only considered for the prediction of the g component. Molecular structures of the complexes are depicted in Figure S1.

Ligand abbreviations: dtc = dimethyldithiocarbamate, acac = acetylacetone; en = ethylenediamine; mnt = maleonitriledithiolate; gly = glycine; kts = 2-keto-3-ethoxybutyraldehyde-bis(thiosemicarbazone); sac = salicylaldehyde imine; im = imidazole; py = pyridine; epa = N,N′-ethylenebis(pyridine-2-aldimine); eta = N,N′-ethylenebis(thiophene-2-aldimine); atpt = 3,4-bis(3-amino-1-thiopropyl)toluene; GGH = glycine–glycine–histidine; GGG = glycine–glycine–glycine; salpn = N,N′-bis(salicylidene)-1,2-propanediamine; (S,S)-mnpala = 2,5,8-trimethyl-5-nitro-3,7-diazanonanedioate; salen = bis(salicylidene)ethylenediamine; bipy = 2,2′-bipyridine; and ttcn = 1,4,7-trithiacyclononane.

To evaluate the performance of the different functionals, the difference (D) between the calculated and experimental g-tensor components was used, defined as follows: D(g) = gcalc – gexp and D(g) = gzcalc – gexp, where g = g⊥ = (g + g)/2 and g = g∥ is the maximum tensor component. In addition, the parameter Δg, defined as the difference between the parallel and the perpendicular g-tensor components: Δg = g – g, is also used as a criterion of the performance of the different functionals according to the following: D(Δg) = Δgcalc – Δgexp. For each examined functional, the mean difference (MD) between the calculated and experimental parameters defined above, is estimated as the average of the D values of the N = 20 Cu(II) complexes. For the case of g, the equation iswhere i runs over the Cu(II) complexes 1–20. The mean per mille difference (‰ MPD) is the average of the per mille differences of complexes 1–20 for each functional and is given by

The calculation of the mean absolute per mille difference ‰ MAPD, by comparing to the ‰ MPD, enables us to identify possible systematic over- or underestimation of the g-tensor components by each functional. It is given by

Finally, the mean percent difference (MPD) and mean absolute percent difference (MAPD) of the g shift of the z component (gs), defined as the deviation from the free-electron value (ge = 2.002319): gs = g – ge are also presented.

The performances of 19 functionals that belong to all different rungs of DFT in terms of the most important MD, ‰ MPD, and ‰ MAPD values for the set of Cu(II) complexes are compared in Table . The complete sets of computed values for all complexes and functionals are given in Tables S4–S22. We first observe that MD(g) values are significantly smaller than MD(g) values for all functionals, which means that the g parameter is much more sensitive than g to the chosen method. The MD(g) follows the same trends as MD(g); therefore, we focus on the average D(g) and D(Δg) values to compare the different methods. Figure summarizes the major conclusions graphically.

Table 2

Comparison of the Performance of the 19 Studied Functionals, Including the B3PW91 Series with Variable Hartree–Fock Exchange Admixture, for the g-Tensor Calculations of the 20 Cu(II) Complexes in Terms of the Evaluation Parameters Defined in Eqs –4 and in the Text

functional	MD(gx,y)	MD(gz)	‰ MPD(gx,y)	‰ MPD(gz)	‰ MAPD(gz)	MAPD(gzs)	MPD(Δg)	MAPD(Δg)
PBE	–0.023	–0.115	–11.3	–51.9	51.9	55.7	–57.8	57.8
TPSS	–0.023	–0.118	–11.2	–53.2	53.2	56.8	–59.4	59.4
TPSSh	–0.018	–0.100	–8.6	–44.9	44.9	47.2	–50.1	50.1
B3LYP	–0.011	–0.072	–5.5	–32.3	32.3	33.5	–35.8	35.8
PBE0	–0.007	–0.055	–3.2	–24.5	24.6	24.6	–26.8	27.8
BHandHLYP	0.015	0.032	7.6	14.8	16.7	24.7	20.7	27.4
B3PW91 (20%)	–0.011	–0.071	–5.3	–31.7	31.7	32.8	–35.1	35.1
B3PW91-25%	–0.007	–0.056	–3.4	–25.0	25.0	25.0	–27.3	28.2
B3PW91-30%	–0.003	–0.039	–1.3	–17.5	19.5	20.1	–18.7	22.8
B3PW91-35%	0.002	–0.021	1.0	–9.3	14.8	16.6	–9.1	18.8
B3PW91-40%	0.007	–0.002	3.4	–0.5	12.4	15.7	1.5	18.7
B3PW91-45%	0.012	0.019	5.9	9.0	13.7	19.6	13.1	22.3
B3PW91-50%	0.017	0.041	8.4	18.9	19.5	28.2	25.7	29.8
B3PW91-55%	0.022	0.064	10.9	29.2	29.2	40.5	39.0	40.0
LC-BLYP	–0.041	–0.175	–20.0	–78.6	78.6	85.1	–85.5	85.5
CAM-B3LYP	–0.033	–0.144	–16.0	–64.7	64.7	69.3	–69.9	69.9
B2PLYP	–0.018	–0.070	–8.8	–31.7	31.7	39.2	–36.8	36.8
DSD-PBEP86	–0.001	0.005	–0.5	2.4	13.8	18.6	5.3	19.9
PBE-QIDH	0.003	0.019	1.7	8.6	14.6	19.3	13.4	23.4

Figure 3

Comparison of the MDs of the calculated g and Δg from the experimental values between different DFT functionals for the set of 20 Cu(II) complexes. The effect of variable exact exchange is shown in the series of results (from 20 to 55% exact exchange) with the B3PW91 functional. The red bars indicate standard deviations.

The magnitude and anisotropy of computed g-factors, the spin population on the Cu ion, and the admixture of the exact (Hartree–Fock) exchange in a density functional are all closely interrelated. This connection provides a convenient perspective for the analysis and discussion of results. PBE and TPSS show significant underestimation of both g and Δg. The ‰ MPD(g) and ‰ MAPD(g) for those functionals have equal absolute values, which shows that the underestimation is systematic. Progressive admixture of the exact exchange is the principal factor that affects the results qualitatively. Thus, increasing percentages of exact exchange in TPSSh (10%), B3LYP and B3PW91 (20%), and PBE0 (25%) progressively limit the underestimation, albeit without eventually correcting it to a satisfactory extent. The BHandHLYP functional, with 50% exact exchange admixture, performs much better, but in this case, systematic overestimation of g and Δg is observed.

Adjusting the percentage of exact exchange in global hybrid functionals provides a way to improve the predictions. To explore this aspect, we chose the B3PW91 as a platform because this functional in its default definition (20% exact exchange) is the method of choice for computing Cu(II) hyperfine coupling constants. We increased the percentage of exact exchange from 20 to 55% in steps of 5% and recomputed the g-tensors for all complexes. The results show a definitive, almost quantitative correlation between the percentage of exact exchange and the errors in g values (Figure ). The average underestimation obtained with the default definition of the functional is essentially eliminated at 40% exact exchange, while further increase in exact exchange results monotonically in overestimation of both the g-shifts and the g-anisotropy. These trends are mirrored on the correlation of the deviation, Δg, of the calculated g from the experimental value with the Cu spin populations computed with the respective functional (Figure S2) that follow closely the change in exact exchange admixture, not in a global sense but specifically for any individual compound. Therefore, the adequate description of covalency in Cu(II) complexes, a critical requirement for the reliable prediction of g-tensors, is where most common functionals falter because of the systematic underestimation of spin localization on copper or underestimation of the paramagnetic contribution to the target quantity. The increase to 40% exact exchange in the B3PW91 functional provides the best average errors in relevant metrics. This is consistent with the prior literature on the subject, particularly with studies by Kaupp and co-workers. B3PW91 has been associated with good performance for transition metal systems[87] and increasing the admixture of exact exchange has been reported to be favorable for g-tensor calculations.[87−89] Other benchmarking studies on 3d, 4d, and 5d transition metal complexes also converged to an optimal value of 40% exact exchange, while the choice of the pure DFT components was deemed less important.[90,91] A comparable approach was used recently for the calculation of EPR parameters in other LPMO enzymes.[22,92]

Figure 4

Correlation of the percentage of exact exchange in the B3PW91 functional with the mean per mille difference (‰ MPD) from experiment of the computed g-tensor component across the reference set of Cu complexes.

Neither the range-separated CAM-B3LYP nor the fully long-range corrected LC-BLYP (which goes to 100% exact exchange at long range) show improved performance over standard global hybrid functionals. On the contrary, they perform on average even more poorly than PBE and TPSS. We notice that these functionals also deviate considerably from the correlation between computed g values and the Cu spin population observed for the other functionals (see Figure S2). The poor performance may imply that the optimal exact exchange percentages and range-separation parameters deviate significantly from the standard definition of these functionals for the target property and system. For example, a higher percentage of exact exchange might be more appropriate as the minimum (short-range limit) value, or a significantly higher percentage should already be reached at shorter interelectronic distances than the default ones. In either case, the two functionals of this family tested here are not applicable to the problem.

Double-hybrid functionals deserve particular mention. These functionals mix perturbation theory (MP2) contribution on top of a Kohn–Sham (KS) determinant with a high percentage of exact exchange.[50] Therefore, virtual KS orbitals are also utilized in an attempt to achieve a balanced description of static and dynamic correlation effects. Because dynamic correlation is known to play an important role in the calculation of g-tensors,[93] double-hybrid functionals should provide a good start. According to our results, the original B2PLYP functional does not perform competitively in agreement with the observations of Sciortino et al.[26] On the other hand, the good performance of PBE-QIDH and, even more so, of DSD-PBEP86 shows that double-hybrid functionals hold promise in terms of achieving systematically improved computational predictions of g-tensors. On average, DSD-PBEP86 is second-best compared to the modified B3PW91. Overall, it does not manage to achieve better error control compared to the tuned global hybrid. Therefore, its use in the present context is not justified in view of the considerable additional computational cost. Nevertheless, the encouraging performance of double-hybrid DFT deserves a more thorough and extensive investigation.

Taking into account the present results as well as previous benchmark studies on calculations of Cu hyperfine coupling constants and g-tensors, the question remains whether one can conceive of a single method that performs equally well for both EPR properties. Such a method remains elusive, but we can surmise that a DFT-based approach would require a very particular handling of exact exchange admixture. A specifically optimized range-separated functional or double-hybrid might be realizable, but one would have to carefully balance such an approach with compromises in the performance of the functional for other properties. It is likely that local, as opposed to global, hybrid functionals, as described by Kaupp and co-workers,[94] represent a more promising development in this direction. For the moment, the results presented above lead us to conclude that doubling the default 20% exact exchange of B3PW91 provides the best solution for the calculation of g-tensors in the reference set of Cu complexes; therefore, this method is adopted for computing g-tensors of the LPMO models, while the default definition of the functional is used for calculations of hyperfine coupling constants.

pH Dependence Study of PlAA10 by EPR Spectroscopy

Following previous work in the group on PlAA10 which evidenced the contribution of two species to the EPR signature of the enzyme recorded at pH 6.5,[15] new spectroscopic measurements were undertaken. X-band EPR spectra were recorded at 120 K using frozen samples of PlAA10 under two different pH conditions at either pH 6.5 in MES buffer or at pH 8.5 in Tris-H2SO4 (2-amino-2-hydroxymethyl-propane-1,3-diol) buffer. EPR parameters simulated by using the EasySpin program package are reported in Table , and the corresponding simulated spectra are shown in Figure together with the experimental ones.

Table 3

X-Band EPR Parameters of PlAA10 at pH 6.5 and 8.5 Putting in Evidence the Existence of Two Signals in a Major/Minor Distribution

		g-tensor			A-tensor (MHz)
pH	species distribution	gx	gy	gz	AxCu	AyCu	AzCu
6.5	major (95%)	2.025	2.103	2.262	220	90	355
	minor (5%)	2.042	2.061	2.230	10	80	560
8.5	major (90%)	2.042	2.061	2.230	10	80	560
	minor (10%)	2.020	2.103	2.262	220	90	355

Figure 5

Experimental (red lines) and simulated (blue dashed lines) cw X-band EPR spectra of a 200 μM solution of PlAA10 at T = 120 K in MES buffer (pH = 6.5, top) and in TRIS buffer (pH = 8.5, bottom).

The EPR signal previously published was reproduced when working at pH 6.5. The major species is present at around 95% with rhombic EPR parameters. The minor species is observed in the remaining 5% of signal with axial EPR parameters. When the pH is increased to 8.5, the observed minor species is exactly reproducing the major species obtained at pH 6.5, that is, the one with a rhombic EPR character. In addition, the major species (around 90%) at pH 8.5 is reproducing very similar EPR values than that of the minor species observed at pH 6.5. The interconversion behavior of PlAA10 is thus solely dependent on the pH. The axial species is major in basic conditions, while the rhombic species predominates in acidic conditions.

The rhombic/axial shifts observed are evidenced in two parameters for both the g- and A-tensors: the maximum absolute value (gmax and |AmaxCu|) and the anisotropic value (Δg, ΔACu), obtained by calculating the difference between the largest and the smallest values (Δg = gmax – gmin, ΔACu = |AmaxCu| – |AminCu|, Table ). In terms of the g-tensor, the gmax is very slightly decreased from rhombic to axial. However, the effect is more important when taking the Δg (0.237 vs 0.188). In terms of the A-tensor, the AmaxCu values are greatly increased in the axial species when compared to the rhombic species (355 vs 560 MHz). A similar effect is observed in ΔACu (265 vs 550 MHz). For our comparative analysis with respect to the constructed active site models, we will compare our computational results to those of the major species, that is, the rhombic EPR parameters for pH = 6.5 and the axial parameters at pH = 8.5.

Table 4

Experimentally Derived EPR Parameters (gmax, Δg, AmaxCu, ΔA) of PlAA10 as a Function of the pH

		g-tensor		A-tensor (MHz)
pH	species distribution	gmax	Δg	AmaxCu	ΔACu
6.5	major (95%)	2.262	0.237	355	265
	minor (5%)	2.230	0.188	560	550
8.5	major (90%)	2.230	0.188	560	550
	minor (10%)	2.262	0.237	355	265

Interestingly, we observe a well-defined superhyperfine pattern in the EPR spectrum recorded at pH = 8.5. To determine to origin of this interaction, we included nitrogen centers in our simulations. To increase the spectral resolution of the narrow hyperfine lines of the nitrogen atoms, we performed second-derivative analysis using the pseudo-modulation method. The result of our second-derivative simulation is presented in Figure .

Figure 6

Experimental (red line) and simulated (green and blue lines with two and three nitrogen atoms, respectively) of the second-derivative cw X-band EPR spectrum of a 200 μM solution of PlAA10 at T = 120 K in TRIS buffer (pH = 8.5).

From Figure , we observe that the simulation obtained when considering three nitrogen atoms with isotropic coupling constants of 45 MHz is in close agreement with the experimental data. This is highlighted by the satellite lines that are found at 322 and 335 mT (see vertical dashed black lines in Figure ). On the contrary, such features are absent when using only two nitrogen atoms in the simulation. This point will be further discussed in light of the results from the computational investigation.

Evaluation of Models for the PlAA10 Active Site

Structures and Energetics

Starting from the parent model A1 (Figure ), we constructed seven types of model of the PlAA10 active site by varying the copper coordination sphere and the protonation states of the ligands with ionizable protons, which are the two oxygen ligands, O1 and O2, and the N-terminal amine of His26. Several geometry optimizations were attempted for each type of model starting from slightly different geometries. These investigations resulted in the seven core structures that are presented in Figure . Key structural parameters of all models are shown in Table . The DFT model A1, whose structural parameters are most consistent with the crystallographic structure (Table ), features a five-coordinated copper in a N3O2 coordination sphere with distorted square-pyramidal geometry. It has O1 and O2 in the aquo form, while the ligating N-terminal histidine bears two protons, that is, [H2O, H2O, −NH2].

Figure 7

Core structures of the optimized models derived from A1.

Table 5

Selected Structural Parameters (Distances in Å, Angles α and β in Degrees, and Structural Indices τ) of the Optimized Models of the PlAA10 LPMO Active Site Considered in the Present Study

model set	structure	Cu–N1	Cu–N2	Cu–N3	Cu–O1	Cu–O2	α	β	τ
A: CuN3O2	A1	1.958	2.143	1.979	2.332	2.363	134.0	174.2	0.67
B: CuN3O	B1	1.961	2.089	1.977	2.177		158.7	164.6	0.26
C: CuN2O2	C1	1.973	3.601	1.961	2.179	2.285	151.1	171.0	0.27
	C2-i	1.996	3.519	1.994	2.100	1.937	165.7	173.9	0.14
	C2-ii	2.028	3.581	2.005	1.902	2.035	159.5	167.7	0.23
	C3-i	2.098	3.509	2.082	1.934	1.893	175.0	169.0	0.11
	C3-ii	2.088	4.000	2.058	1.930	1.909	156.5	160.9	0.30
exp.	6T5Z	2.038	2.149	1.997	2.404	2.322	131.7	166.1	0.57

In the presentation of the derivative models, we begin with model set B that includes structures with a N3O coordination sphere. Structure B1 is obtained by decoordination of a single water molecule from model A1, which results in [H2O, −NH2] protonation states. B2 is obtained from B1 by deprotonation of the single copper-coordinating ligand, which results in the [OH–, −NH2] protonation state. It is important to note here that geometry optimizations of structures derived from A1 with deprotonation of a single water ligand, that is, [H2O, OH–, −NH2], leads to a distorted square pyramidal structure with a water axial ligand that decoordinates from the metal and results in structure B2.

Structures that belong to model set C possess a N2O2 coordination sphere. Parent structure C1 is obtained by decoordination of the N-terminal histidine and protonation of the amine group resulting in protonation states [H2O, H2O, decoordinated −NH3+]. In C1, the NH3+ group forms a hydrogen bond with a ligated water. In models C2 (labeled C2-i and C2-ii), the copper bears a water and a hydroxo ligand [H2O, OH–, decoordinated −NH+]. Model C2-i is obtained from parent structure C1 by removing a proton from the water ligand that is hydrogen-bonded to the NH3+ group. A1 and C2-i are structural isomers. Model C2-ii results from deprotonation of the decoordinated NH3+ group and has protonation states [H2O, OH–, decoordinated −NH2]. Finally, in models C3, the Cu ion bears two hydroxo ligands [OH–, OH–, decoordinated −NH(+)]. C3-i [OH–, OH–, decoordinated −NH3+] is obtained by the deprotonation of one of the water ligands of C2-i and is therefore a structural isomer of C2-ii. Finally, model C3-ii [OH–, OH–, decoordinated −NH2] is obtained from C3-i by deprotonation of the −NH3+ group.

In terms of computed relative energies, A1 and C2-i are isomers and C2-i is calculated to be 6.0 kcal mol–1 higher than A1. Among the structures that result from C2-i deprotonation, C2-ii is favored by 9.8 kcal mol–1 relative to C3-i, which shows that in model C2-i, the presence of the hydrogen bond favors the deprotonation of the primary amine. Model C2-ii is estimated at 11.3 kcal mol–1 higher than B2 when directly compared. Therefore, the present results show that model B2 is more stable than C3-i by 21.1 kcal mol–1. This is a crucial result because in a recent study by Lindley et al.[22] the axial EPR signal of the BlAA10 LPMO formed at pH 8.5 was assigned to a decoordinated histidine brace structure similar to the present model C3-i. This scenario is clearly disfavored in the present case. After evaluation of the models with respect to energetics, we proceed with evaluation of their EPR parameters against the available experimental EPR data.

EPR Parameters and Electronic Structure

The calculated g- and A-tensor components for all DFT models, using the B3PW91 functional with 40 and 20% exact exchange admixture, respectively, are shown in Table . EPR parameters of structures A1, B1, and C1, which have protonated water ligands, maintain a rhombic character. Structures B2 and C2, where copper is coordinated to a single hydroxo ligand, produce axial EPR parameters. Structures C3, where both water ligands are in their hydroxyl form, also exhibit rhombic character. To better assess the agreement with the experiment, we examined the rhombic/axial shift parameters (gmax, AmaxCu, Δg, and ΔACu) of all considered structures, as shown in Table . The rhombic character of the A1, B1, C1, and C3 EPR parameters is evident from the AmaxCu and ΔACu values which range between 309 and 446 MHz and between 269 and 375 MHz, respectively. On the contrary, structures B2 and C2 have AmaxCu and ΔACu values which range between 569 and 604 MHz and between 546 and 599 MHz, respectively. At the same time, B2 and C2 have slightly smaller gmax and Δg values than the rest of the models.

Table 6

Calculated g-Factors and Cu Hyperfine Coupling Constants (MHz) for all Models

model set	model	gx	gy	gz	AxCu	AyCu	AzCu
A: CuN3O2	A1	2.032	2.138	2.257	101	–299	370
B: CuN3O	B1	2.053	2.092	2.254	71	165	–446
	B2	2.059	2.077	2.243	23	–37	–569
C: CuN2O2	C1	2.041	2.160	2.329	15	181	309
	C2-i	2.064	2.066	2.243	–15	–16	–597
	C2-ii	2.057	2.078	2.246	5	–53	–604
	C3-i	2.047	2.088	2.248	85	200	–436
	C3-ii	2.045	2.100	2.258	91	227	–393
PlAA10
pH 6.5	rhombic	2.025	2.103	2.262	220	90	355
pH 8.5	axial	2.042	2.061	2.230	10	80	560

Table 7

Axial/Rhombic Shifts [gmax (=g), AmaxCu (=ACu) and Δg, ΔACu] for all Models

	model	gmax	Δg	|AmaxCu|	ΔACu
A: CuN3O2	A1	2.257	0.225	370	269
B: CuN3O	B1	2.254	0.201	446	375
	B2	2.243	0.184	569	546
C: CuN2O2	C1	2.329	0.288	309	293
	C2-i	2.243	0.179	597	582
	C2-ii	2.246	0.189	604	599
	C3-i	2.248	0.201	436	351
	C3-ii	2.258	0.213	393	303
PlAA10
pH 6.5	rhombic	2.262	0.237	355	265
pH 8.5	axial	2.230	0.188	560	550

The results confirm that structure A1 best reproduces the major species in the EPR spectrum of PlAA10 at pH = 6.5. Among the models B2 and C2 that produce axial EPR parameters, structure B2 has AmaxCu and ΔA values closer to the experiment than both protonated (C2-i) and deprotonated (C2-ii) C2 forms. Overall, the computed EPR parameters of structures A1 and B2 are in best agreement with experimental values obtained for PlAA10 at pH = 6.5 and 8.5, respectively.

The observed rhombic and axial spin-Hamiltonian parameters can be qualitatively rationalized based on the geometry of the copper coordination sphere. The singly occupied molecular orbital (SOMO) of square planar and square-pyramidal copper complexes has mostly d character. Distortion from the perfect square planar and square-pyramidal geometry results in increasing extent of mixing of the d orbital into the d SOMO. The rhombic spin-Hamiltonian parameters arise from this d orbital mixing. The extent of geometrical distortion can be expressed using the geometry index τ (defined in Scheme S2),[95] given in Table for each structure. Structure A1 has a geometry index value of 0.67 and therefore shows a more pronounced deviation from the square-pyramidal geometry compared to B2 which has a τ value of 0.34. The d orbital mixing into the d SOMO of A1 is clearly seen on the molecular orbital diagram of the SOMO shown in Figure . On the contrary, the SOMO of B2 has a dominant d character, which explains the computed—and experimentally observed—axial EPR parameters.

Figure 8

SOMOs of models A1 and B2 (for clarity, only the immediate environment of the Cu center is depicted).

A final point concerns the 14N isotropic hyperfine constants. The analysis of spectroscopic data presented above established the presence of three nitrogen centers in the immediate surrounding of the unpaired electron. However, in line with structure-based expectations, DFT calculations on models C2-i and C3-i confirm that in these coordination geometries, only the coordinating imidazole nitrogen atoms (N1 and N3) have large isotropic hyperfine coupling values (ca. 35 MHz for C2-i and 45 MHz for C3-i). In both cases, the computed hyperfine coupling constant for the decoordinated N2 nitrogen nucleus is negligible (ca. 0.5 MHz). Therefore, these models cannot be reconciled with the spectroscopic observations, in contrast to model B2 that indeed has three large 14N hyperfine couplings as required by the experiment (Table S23). These results further support that the major species observed in high pH conditions is consistent with a Cu(II) ion bound to three nitrogen centers and a hydroxo ligand featuring a N3O coordination sphere.

Discussion

The theoretical methodology we developed recently for Cu(II) hyperfine coupling constants[25] and, in this work, for accurate g-tensors, can be reliably applied in the present case to provide a solid identification of the coordination sphere of copper species when experimental EPR data are available. This work is thus filling a methodological gap that was recently pointed out and prevented a clear identification of unknown copper coordination spheres in proteins.[22,24,92] Armed with these methodologies, we focused on resolving the ambiguous, pH-dependent EPR signals from the LPMO enzyme PlAA10. In the construction of the presented set of active site models, many plausible chemical modifications that could occur when changing the experimental pH conditions were considered.

Our results point to a N3O2 coordination sphere for the species observed at low pH, that is, model A1 with two H2O ligands on Cu. This model best reproduces the experimental rhombic EPR signal that can be attributed to the PlAA10 active site from the 6T5Z X-ray crystal structure used as a starting point for our computational investigation. The rhombic-to-axial shift of the observed EPR signal upon increasing pH indicates that the axial signal is induced either by deprotonation of the active site or by a pH-induced structural change of the protein. Our results show that the axial signal can be assigned uniquely to a structure with N3O coordination where the copper ion is coordinated to a single hydroxyl ligand, that is, model B2 of the present work. Neither the mere decoordination of a single water ligand (B1) nor the decoordination of the protonated amine(C1) creates the shift to axial EPR parameters. These models exhibit well defined rhombic signatures. Alternative structures of N2O2 coordination and two hydroxyl ligands on copper, C3-i and C3-ii, also produce rhombic EPR parameters. Only models B2, C2-i, and C2-ii produce axial EPR signals. Of these candidate structures, models C2-i and C2-ii are both energetically unfavorable and inconsistent with the experimental requirement for three copper-ligated nitrogen nuclei. Therefore, model B2 provides the only convincing and fully experimentally consistent interpretation of the axial EPR signal produced at the higher pH value.

The proposed hypothesis of an A1 and B2 equilibrium is consistent with a recent EPR study of the PpAA10 enzyme, where the rhombic signal was assigned to the distorted square pyramidal structure described by crystallography and the axial signal at higher pH was assigned to the formation of a nearly planar species that resulted from the replacement of a copper-coordinating water molecule with a hydroxide or Cl– ion along with decoordination of a water ligand.[23] On the other hand, the A1 and B2 equilibrium is in apparent contradiction with a recent work on the BlAA10 enzyme which investigated the protonation states of the active site as function of pH by means of EPR spectroscopy and DFT calculations.[22] The preferred description of the low pH species in that work was of a N3O2 model comprising one water and one hydroxy ligand on Cu, even though comparison of the computed EPR parameters with the experimental data in acidic conditions shows non-negligible deviation. For the major species observed at moderate-high pH, the structure proposed contained a four-coordinated copper center in an N2O2 coordination sphere with two hydroxy ligands. The model featured decoordination of the amine from the N-terminal histidine along with its protonation. However, the reported simulation of the observed superhyperfine pattern for the species associated with the axial signal at moderate-high pH was carried out including only two nitrogen nuclei and does not achieve satisfactory agreement with the experiment.[22] This type of interpretation receives no support in the present case of PlAA10. While we acknowledge that the enzymes are distinct and that different phenomena might be at play, the excellent agreement between experimental and computed EPR parameters for the models presented in this work allow us to propose an alternative structure for the axial species observed in basic conditions. Our suggested N3(OH) structure, model B2, displays a fully coordinated histidine brace and a single hydroxy ligand, an assignment strongly supported by our refined simulations of the experimental superhyperfine pattern: using the pseudo-modulation method, we show that the agreement between the experimental and the simulated second-derivative EPR spectrum can only be reached upon inclusion of three nitrogen centers. This conclusion is further supported by DFT computations on both energetic and spectroscopic grounds. Therefore, the combination of EPR measurements with the properly calibrated DFT methodology employed in this work provide a definitive picture of the experimentally observed pH-dependent species and allows us to identify unambiguously the coordination spheres of the copper active sites of the PlAA10 family, with implications for re-evaluating relevant assignments in other AA10 LPMOs.

Conclusions

In the present work, we described an accurate quantum chemical approach for the calculation of Cu(II) g-tensors and its application together with our recently developed methodology for copper hyperfine coupling constants to the structural assignment of EPR signals in the LPMO from Photorhabdus luminescens, PlAA10. Based on the benchmarked accuracy of our applied methodology in the prediction of the relevant EPR parameters, we can safely assign the species observed at pH = 6.5 and 8.5. The low pH structure is a distorted square-pyramidal copper center with two water molecules displaying rhombic EPR parameters, whereas the higher pH species is a distorted square-planar copper bound to a single hydroxy ligand leading to axial EPR parameters. This assignment is fully consistent with all aspects of the quantum chemical calculations and with the spectroscopic observations, confidently excluding structural alternatives on the basis of copper EPR parameters and relative energetics of isomeric forms, as well as on the assignment of the number and magnitude of nitrogen superhyperfine parameters. The successful application of our approach provides a solid framework to establish correlations between structural features and spectroscopic data on LPMOs. By combining EPR spectroscopy and quantum chemistry, this work contributes to an accurate description of the biologically relevant species of AA10 LPMOs and allows us to reach a better understanding of the enzyme, while the computational strategy is broadly applicable for the reliable prediction of EPR parameters of intermediates involved in the catalytic mechanism of LPMO and in copper proteins in general.