Accurate structure and dynamics of the metal-site of paramagnetic metalloproteins from NMR parameters using natural bond orbitals.

A natural bond orbital (NBO) analysis of unpaired electron spin density in metalloproteins is presented, which allows a fast and robust calculation of paramagnetic NMR parameters. Approximately 90% of the unpaired electron spin density occupies metal-ligand NBOs, allowing the majority of the density to be modeled by only a few NBOs that reflect the chemical bonding environment. We show that the paramagnetic relaxation rate of protons can be calculated accurately using only the metal-ligand NBOs and that these rates are in good agreement with corresponding rates measured experimentally. This holds, in particular, for protons of ligand residues where the point-dipole approximation breaks down. To describe the paramagnetic relaxation of heavy nuclei, also the electron spin density in the local orbitals must be taken into account. Geometric distance restraints for (15)N can be derived from the paramagnetic relaxation enhancement and the Fermi contact shift when local NBOs are included in the analysis. Thus, the NBO approach allows us to include experimental paramagnetic NMR parameters of (15)N nuclei as restraints in a structure optimization protocol. We performed a molecular dynamics simulation and structure determination of oxidized rubredoxin using the experimentally obtained paramagnetic NMR parameters of (15)N. The corresponding structures obtained are in good agreement with the crystal structure of rubredoxin. Thus, the NBO approach allows an accurate description of the geometric structure and the dynamics of metalloproteins, when NMR parameters are available of nuclei in the immediate vicinity of the metal-site.

Introduction

Understanding the mechanisms that lead to the biological function of metalloproteins requires detailed knowledge of the metal-site structure. So far, metal-site structures of metalloproteins have been studied experimentally mainly by techniques such as X-ray crystallography, X-ray absorption spectroscopy,[1−3] extended X-ray absorption fine structure,[4−7] and EPR.[8−10] These techniques provide information about the geometric and electronic structure of the metal-site in the crystal phase or in frozen liquids. However, proteins normally function in solution where the flexibility and dynamics of the proteins often are essential for their biological activity. Characterization of both the dynamics and the structure of the metal-site in aqueous solution at the atomic level is therefore important for the elucidation of functional mechanisms of metalloproteins.

Paramagnetic nuclear magnetic resonance (NMR) spectroscopy is a unique tool for providing information about metal-site structure and dynamics of paramagnetic metalloproteins such as iron–sulfur proteins,[11−13] copper proteins,[14−17] heme proteins,[18,19] and other metalloproteins.[20,21] Thus, the nuclear spins act as structural probes yielding information about the geometric and electronic structures of the metal-site in solution through their interaction with the unpaired electron spin.

The chemical shift of nuclei close to a paramagnetic metal-site is affected strongly by the hyperfine Fermi contact shift.[22] The Fermi contact shift arises from a through-bond electron–nuclear scalar coupling and is proportional to the unpaired electron spin density at the nuclear position. Therefore, it provides information about the delocalization of the unpaired electrons onto the ligand nuclei. Theoretical quantum chemical calculations, experimental NMR, and ENDOR studies have shown that the Fermi contact shift depends on the geometry of the protein in the immediate vicinity of the metal-site.[12,17,23,24] Therefore, it also gives information about the structure and structural changes close to the metal-site. The nuclear longitudinal paramagnetic relaxation enhancement (PRE) is complementary to the Fermi contact shift. It stems from a through-space electron–nucleus dipole interaction, which depends on the geometric position of the nucleus relative to the total unpaired electron spin distribution of the metal-site. The PRE can therefore provide long-range distance information in proteins[11,25−30] and information about spatial distribution of the electron spin in metalloproteins.[14,16,17]

Despite the large and versatile amount of information provided by the nuclear–electron interactions, the use of paramagnetic NMR in the determination of the structure and dynamics of metal-sites of metalloproteins has so far been limited. This is mainly due to the fact that a proper description of the unpaired electron spin density in proteins requires a large basis set, which makes restrained dynamics simulations and accurate structure determinations computationally highly demanding.

Here, we show how a natural bond orbital (NBO)[31] description of the unpaired electron spin density of metalloproteins can alleviate this problem. First, we show that the majority of the unpaired electron spin density of metalloproteins can be modeled by a small number of NBOs that reflect the chemical bonding environment of the protein. Second, we show how accurate NMR parameters can be calculated from a small and transferable basis set that is based on the NBOs. For 1H nuclei, the PRE can be calculated accurately from the metal–ligand NBOs, as for example the PRE of 1H nuclei of the iron–sulfur protein rubredoxin that can be calculated using the NBOs that describe the Fe–S bond. For heavy atoms, such as 13C and 15N, local NBOs that are centered on the heavy atom must be included.

The model proteins used here are the iron–sulfur protein rubredoxin and the blue copper protein plastocyanin. Paramagnetic interactions in both proteins have been studied extensively, both experimentally by NMR and theoretically by density functional theory (DFT) calculations.[11,12,32,33] Moreover, these two metalloproeins are very different both with respect to symmetry, structure, and electronic spin. Accurate paramagnetic NMR parameters are calculated using the small NBO basis set for both plastocyanin and rubredoxin, which suggests that the approach developed below is applicable to a wide range of metalloproteins. Overall, combining the NBO description of the electron spin with experimental NMR parameters allows an accurate characterization of the dynamics and structures of metalloproteins in solution, as shown below.

Theory

Paramagnetic Relaxation Enhancement

The longitudinal paramagnetic relaxation enhancement (PRE) contains contributions from the Fermi contact relaxation,[34−36] the Curie spin relaxation,[37,38] and the dipolar relaxation.[39] The Fermi contact contribution to longitudinal paramagnetic relaxation is in general small even for hyperfine shifted signals.[40] Specifically, the contact contribution is less than 0.5% of the dipolar relaxation for all of the nuclei used in this study[11,41] and is therefore neglected. The experimentally measured longitudinal PRE of the nuclei depends on the unpaired electron spin density as follows:[38,39,42,43]where ge is the electron g-factor, μB is the Bohr magneton, S is the spin quantum number of the unpaired electron, μ0 is the magnetic permeability of free space, γI is the gyromagnetic ratio of the nucleus, τc,1 is the correlation time of the electron–nucleus dipolar interaction, τR is the rotational correlation time of the protein, and ⟨S⟩ is the thermal average of the electron spin magnetization.[38] The effect of Curie spin relaxation,[38] described by the term modulated by τR in eq 1, can often be neglected. For example, for the high-spin Fe(III) protein rubredoxin (blue copper protein plastocyanin) studied here, the contribution to R1p from Curie spin relaxation is about 2.5% (0.03%). Even in studies where Curie spin relaxation cannot be neglected, the longitudinal relaxation enhancement is proportional to reff–6, where ξ is the proportionality constant. Moreover, cross-correlated relaxation between Curie-spin relaxation (N···⟨S⟩) and dipole relaxation with a directly bond nuclei (N···H) can be neglected if (1) H or N relaxes faster than the scalar coupling between H and N or (2) the relaxation rate of N is measured on a coherence that is in-phase with respect to H. Both criteria are fulfilled here, and cross-correlated relaxation is therefore neglected. The effective nucleus–electron distance is given by the unpaired spin distribution as follows:[42]where ℱ̂2ν(r) = ∥r∥–3Y2ν(r/∥r∥) are the spatial components of the electron–nucleus dipolar operator centered at the nucleus at the position r′, and Y2ν(r/∥r∥), ν = −2,...,2, are the spherical harmonics. Furthermore, F̂2ν(r′) is the matrix representation of the dipolar operator in the basis of the spin density, and P is the spin-only Fock–Dirac density matrix.

Experimentally, the longitudinal PRE of a given nucleus is obtained as the difference between the relaxation rate of the nucleus in the paramagnetic and the corresponding diamagnetic species. For example, for the nuclei of the blue copper protein plastocyanin, a longitudinal PRE is obtained as R1p = R1(Cu(II)) – R1(Cu(I)), R1(Cu(II)) and R1(Cu(I)) being the longitudinal relaxation rates measured in the paramagnetic copper(II) protein, and in the copper(I) analogue diamagnetic protein, respectively. Subsequently, the effective distance, reff, can be calculated from the PRE once the correlation time, τc,1, or the proportionality constant, ξ, is known. Overall, the effective distance depends on the distribution of the total unpaired electron spin in the molecule and therefore provides information on the overall electron spin distribution via ρ(r) = ρα(r) – ρβ(r), as well as information on the position of the nucleus relative to the distribution via ℱ̂2ν(r – r′).

Electron Spin Density and Hyperfine Shift

The unpaired electron spin density causing the paramagnetic relaxation enhancement and given in eq 3 is located either in the metal–ligand orbitals or in the local orbitals of the relaxing nucleus. Even though only a small fraction of the spin density is in the local orbitals of the relaxing nucleus, this density is spatially very close to the nucleus and can therefore give a significant contribution to the relaxation. It is our goal to estimate the spin density of the local orbitals from the hyperfine chemical shift of heavy atoms and subsequently “correct” the experimentally derived relaxation enhancements such that these enhancements report exclusively on the position of the nucleus relative to the metal–ligand center.

The chemical shift observed for nuclei of a paramagnetic protein consists of contributions from the hyperfine interaction, δhyp, and a diamagnetic contribution. The hyperfine shift, δhyp, has in general two contributions, that is, the Fermi contact shift, δcon, and the dipolar (through space) pseudocontact shift, δpcs:For nuclei in the immediate vicinity of a paramagnetic metal ion, the Fermi contact shift is usually much larger than the pseudocontact shift. This holds particularly for high-spin Fe(III) systems investigated here.[12,22] We find that δpcs/δhyp < 0.085 for all of the nitrogen nuclei of rubredoxin included in the analysis below with the average ⟨δpcs/δhyp⟩ = 0.025 ± 0.025. An upper limit for δpcs was calculated from the susceptibility tensor of Fe(III) rubredoxin (Δχax = 5.3 × 10–28 cm3, Δχrh = 2.1 × 10–28 cm3),[44] and δhyp was obtained experimentally.[45]

The Fermi contact shift, δcon, relates to the unpaired electron spin density through:[22,46]where k is the Boltzmann constant, T is the absolute temperature, and δ̂(r′) is the matrix representation of the Dirac-delta function centered at the nuclear position. The nuclear position, r′, is given in the coordinate frame of {|ϕ⟩}, which is the basis set in which the electron spin density is described.

From eq 5, it is seen that the peak position in an NMR spectrum of hyperfine-shifted nuclei is a direct measure of the electron spin density at the position of the nucleus, if the pseudocontact shift δpcs is much smaller than δcon or if δpcs can be calculated. Below, we use the Fermi contact shift to calculate the spin density of the local orbitals to take the relaxation effect of the local orbitals into account. In general, the contribution to the PRE from the local orbital can be neglected, when the pseudocontact shift is similar to or larger than the Fermi contact shift. For the nuclei of the ligand binding residues of rubredoxin, the Fermi contact shift is obtained as δcon = δFe(III) – δdia, where δFe(III) is the chemical shift in the paramagnetic high-spin Fe(III) protein[45] and δdia is the chemical shift in the analogue diamagnetic protein, which in this case was estimated from the random coil value.[47]

Spin Population Analysis of Natural Bonding Orbitals

Natural bond orbital (NBO) analysis is based on local block eigenvectors of the one-particle density matrix and has emerged as a technique for studying hybridization and covalency effects in polyatomic wave functions.[48−50] In short, the NBO between two atoms A and B is derived by first transforming the one-particle density matrix to the natural atomic orbital basis and subsequently only considering the part of the density matrix, PAB, that involves natural atomic orbitals centered at A or B. The matrix PAB is then depleted for orbitals such as core orbitals and lone pairs that are only centered at one of the two atoms. An eigenvector of the depleted PAB density matrix with an occupation number (eigenvalues) over a certain threshold constitutes a pre-NBO (pNBO). The calculated eigenvectors can be used to transform the original set of basis functions (often Gaussian type atomic orbitals) into the basis set of pNBOs, core orbitals, and lone-pairs. The set of NBOs corresponds closely to the picture of localized bonds and lone pairs, and ab initio wave functions transformed into pNBOs are found to be in good agreement with the Lewis structure concepts of bond hybridization.[31] The NBOs are the set of orbitals obtained by orthogonalizing the two-center pNBOs.

Results and Discussion

Unpaired Electron Spin Density and Metal–Ligand pNBOs

The pNBO basis set will serve as a starting point for a robust and transferable modeling of the unpaired electron spin density in metalloproteins. Whenever possible, the two-center preorthogonalized natural bond orbitals (pNBO) will be used instead of the fully orthogonalized NBOs. The genuine two-center pNBOs are easily transferred from one molecule to another because they are linear combinations of atomic orbitals centered only at the two nuclei that are connected with a bond. Thus, as shown below, the metal–ligand pNBO description of a metal-site model complex can be transferred from the small model complex to a full metalloprotein and thereby used to describe the unpaired electron spin density in the protein. Most of the unpaired electron spin density of a paramagnetic metalloprotein occupies metal and ligand atomic orbitals. Intuitively, when the unpaired electron spin density is described in the pNBO basis set of orbitals, the metal–ligand antibond orbitals and metal and ligand lone-pairs contain the majority of the unpaired electron spin density. To verify that the metal–ligand pNBOs and lone-pairs alone form a good basis for the unpaired electron spin density in metalloproteins, we consider two small model complexes: (1) the spin 5/2 iron–sulfur model complex, Fe(SCH3)4–, and (2) the spin 1/2 copper model complex, Cu(Im)2(SCH3)S(CH3)2+. The structures of the models were based on X-ray crystal structures of the iron–sulfur protein rubredoxin (PDB code 4RXN(51)) and the blue copper protein plastocyanin (PDB code 1PLC(52)); these model complexes are shown in Figure 1.

Figure 1

The two small model complexes used to study the unpaired electron spin density of rubredoxin and plastocyanin: (a) the Fe(SCH3)4– model complex, which is used as a model of rubredoxin, and (b) Cu(SCH3)(Im)2S(CH3)2+, which is a model of plastocyanin.

Below, we calculate the spin occupation numbers for the pNBOs of the two small model complexes and demonstrate (1) that the majority of the unpaired electron spin density occupies metal–ligand pNBOs and lone-pairs and (2) that only very few pNBOs have significant occupancy. Finally, we use the pNBOs that are highly spin-occupied to calculate NMR parameters.

Identifying pNBOs with Unpaired Electron Spin Density

Electron densities, ρα and ρβ, of the two small model complexes were obtained by quantum chemical calculations at the UB3LYP/6-31G* level. Subsequently, an NBO analysis was performed,[49] and the Fock-Dirac spin density matrix in the atomic orbital (AO) basis, PAO, was transformed to the two-center pNBO basis (β-spin) by the transformation matrices provided by the NBO analysis program as implemented in Gaussian 03 and described in Material and Methods. Gross occupation numbers, n, for the unpaired electron spin density of the individual basis functions, n = (Sℬ–1Pℬ), were calculated in the bases ℬ = AO and ℬ = pNBO, respectively. Here, Sℬ is the overlap matrix of the basis set in question, while Pℬ is the Fock–Dirac spin density matrix in the same basis.

Figures 2 and 3 show the gross spin occupation numbers, n = (Sℬ–1Pℬ), of the unpaired electron spin density in the AO and the pNBO bases and for the two model complexes. As it appears from the figures, only a small number of the two-center pNBO basis functions are significantly spin-populated; that is, the unpaired electron spin density can be described, to a good approximation, by a few basis functions, which relate to the chemical bonding environment.

Figure 2

The gross spin occupation numbers, n of the unpaired electron spin, (Sℬ–1Pℬ), for the Fe(SCH3)4– model complex in two different basis sets: (a) the two-center pNBO basis and (b) the atomic orbital basis, PAO (in the 6-31G* basis). Sigma antibond orbitals are indicated as σ*, lone pairs as LP, and Rydberg orbitals as RY. Only a few orbitals have significant spin occupation in the pNBO basis, and the occupied orbitals are strongly correlated to the chemical bonding environment.

Figure 3

The gross spin occupation numbers of the unpaired spin, n = (Sℬ–1Pℬ), for the Cu(SCH3)(Im)2S(CH3)2+ model complex in two different basis sets: (a) the two-center pNBO basis and (b) the atomic orbital basis, PAO (in the 6-31G* basis). Sigma and π bond orbitals are shown as σ and π, respectively, while sigma and pi antibond orbitals are shown as σ* and π*, respectively, and LP indicates a lone pair. Only a few orbitals in the pNBO basis are occupied by unpaired electron spin.

The gross occupation numbers for the unpaired electron spin of the metal–ligand pNBOs are summarized for the two model complexes in Table 1. The two-center metal–ligand pNBOs account for approximately 90% of the total unpaired electron spin density for both the rubredoxin and the plastocyanin complex. Also, primarily valence antibond orbitals (σ*, π*) and lone pairs have significant population of unpaired electron spin, whereas only a small amount of the unpaired electron spin occupies the high-energy Rydberg (RY) orbitals.

Table 1

Gross Occupancies of Unpaired Electron Spin in the Metal–Ligand pNBOs

rubredoxina				plastocyaninb
pNBOc	(S–1P)ii	∑(S–1P)iid	type	pNBOc	(S–1P)ii	∑(S–1P)iid	type
62	0.87	0.87	LP(Fe)e	261	0.55f	0.55	σ*,π*(Cu–S)
61	0.82	1.70	LP(Fe)e	262	0.28f	0.83	σ*Cu–S)
179	0.62	2.31	σ*(Fe–S)	12	0.03f	0.85	σ,π(Cu–S)
177	0.62	2.93	σ*(Fe–S)	11	0.02f	0.87	σ,π(Cu–S)
176	0.60	3.54	σ*(Fe–S)	75	0.02	0.89	LP(N)
178	0.60	4.14	σ*(Fe–S)	68	0.02	0.91	LP(N)
56	0.09	4.23	LP(S)e
60	0.09	4.31	LP(S)e	261, 262	0.83	0.83	π*(Cu–S)
54	0.07	4.39	LP(S)e	11, 12	0.04	0.87	π(Cu–S)
58	0.07	4.46	LP(S)e	75	0.02	0.89	LP(N)
				68	0.02	0.91	LP(N)

Rubredoxin is represented by the Fe(SCH3)4– model complex.

Plastocyanin is represented by the Cu(SCH3)(Im)2S(CH3)2+ model complex.

The numbers refer to the x-axis of Figures 2 and 3.

The accumulative gross spin density.

The pNBOs number 54, 56, 58, 60, 61, and 62 of the Fe(SCH3)4– model complex correspond to π antibonds between the iron and the four sulfur atoms.

The mixed π and σ antibond orbitals number 261 and 262 can be transformed into σ and π antibond orbitals with spin occupancies of −1.1 × 10–3 and 0.83, respectively. The same transformation of the corresponding bond orbitals, that is, 11 and 12, results in σ and π bond orbitals with spin occupancies of 3.7 × 10–3 and 0.040, respectively (see text). The populations of these transformed orbitals are shown in the last four rows of columns 5−8.

From Figures 2 and 3, it is seen that the unpaired electron spin population of the bond orbitals is small as compared to the population of the antibond natural orbitals. Thus, the tight connection between the chemical bonding environment and the pNBOs is evident from the types of pNBOs that are populated, because one would expect that the bonding orbital will be fully occupied with one α and one β electron, while the antibonding orbital is only partly occupied. For the rubredoxin model complex, the spin-populated pNBOs fall in two groups; four correspond to Fe–S σ antibond (σ*), while six pNBOs correspond to π Fe–S antibond formed by lone pairs at the iron and the sulfur atoms.

For the plastocyanin model complex, the two most spin-populated pNBOs correspond to a mixture of Cu–S π and σ antibond orbitals with spin occupancy of 0.55 and 0.28, respectively. The submatrix of the spin density matrix that corresponds to the two Cu–S orbitals can be transformed into two orbitals with σ and π antibond character. In this “refined basis”, the gross spin occupation of the π antibond orbital is 0.83, while the occupation of the corresponding σ antibond orbital is −1.1 × 10–3. Thus, the unpaired electron spin density occupies only the antibond π orbital and does not occupy the σ* orbital. This is in agreement with the UV–vis spectra of plastocyanin, which show that the unpaired electron spin predominantly is in the π antibond Cu–S orbital.

The calculations above show that the unpaired electron spin density can be modeled by the pNBOs that reflect the chemical bonding environment. In Figure 4, the spin-populated pNBOs of the two model complexes are shown schematically. The β-LUMO accounts for the total spin density of the plastocyanin model complex and is shown in Figure 4h. It is evident that most of the β-LUMO spin density (given by the 286 atomic basis functions) can be represented accurately by a linear combination of the three pNBOs sketched for plastocyanin in the figure. The nitrogen lone-pairs shown in Figure 4f and g are very similar in shape, which is expected due to the high transferability of the pNBOs.

Figure 4

Schematic representation of the pNBOs of Figures 2 and 3 with high spin occupation numbers. (a) Lone pair at the sulfur, pNBO number 54 in Figure 2; (b) lone pair at the ferric, pNBO number 61 in Figure 2; (c) lone pair at the ferric, pNBO number 62 in Figures 2 and 3; (d) the σ* antibond orbital between the sulfur and the ferric ion, pNBO number 176 in Figure 2; (e) the π* orbital obtained from pNBO numbers 261 and 262 in Figure 3 as described in the text; (f) lone pair at the nitrogen, number 68 in Figure 3; (g) lone pair at the nitrogen, number 75 in Figure 3; and (h) the β-LUMO of the plastocyanin model complex (286 basis functions). For blue copper proteins, the square of the β-LUMO is nearly identical to the spin density. As indicated, the β-LUMO can be obtained by a linear combination of the pNBOs shown in (e), (f), and (g). The density at the protons of SCH3– in (h) corresponds to the pNBO 265 in Figure 3, which is a C–H antibond orbital.

Considering the large difference between the two model complexes studied here, the metal–ligand pNBOs are likely, in general, to account for the majority of the unpaired electron spin density of metalloproteins. This implies that properties that depend on the overall spin distribution can be calculated accurately from the few metal–ligand pNBOs instead of using the full basis set, thus speeding up calculations of such parameters significantly. Below, we will show that the PRE of protons, which primarily depend on the overall unpaired electron distribution, can be calculated accurately from the metal–ligand pNBOs.

Paramagnetic Relaxation of the Protons

The longitudinal dipolar PRE, eq 1, of protons is primarily caused by the unpaired electron spin density located in the metal–ligand orbitals. Although unpaired electron spin may also occupy local orbitals of the proton (see Figure 4h), these orbitals are predominantly s-type orbitals that do not contribute to dipolar relaxation because they are spherical symmetric. Overall, the longitudinal paramagnetic relaxation of protons in metalloproteins is the simplest case to analyze in terms of unpaired electron spin density and can, intuitively, be calculated using the metal–ligand pNBOs shown in Table 1. Therefore, paramagnetic proton relaxation is our starting point in the search for a general, small, and transferable basis for the unpaired electron spin density, which allows accurate calculations of paramagnetic NMR parameters.

Two large model complexes are considered to evaluate in detail the dependence of the PRE of protons on the unpaired electron spin density. The structures of these model complex are shown in Figure S1. As compared to the model complexes shown in Figure 1, the larger model complexes include all of the amino acids in the vicinity of the metal-site and the hydrogen bonds that stabilize the metal charge. Thus, a 144-atom complex was used to model the metal-site environment of plastocyanin,[53] while a 104-atom complex was used to model the metal-site environment of rubredoxin[11] (see Materials and Methods). DFT calculations and NBO analyses of the spin density form the basis for the evaluation of the PRE. The large model complexes serve as “true” models for the proteins, and the NMR parameters calculated from the full DFT spin densities of the complexes serve as “true” values. Moreover, a comparison of the derived NMR parameters with the experimentally measured parameters serves as a test of the use of metal–ligand pNBOs for calculation of NMR parameters.

In Figure 5, we compare PREs calculated from different models of the unpaired electron spin density using the corresponding effective distances (eqs 1 and 2) of the nuclei in the two proteins rubredoxin and plastocyanin as proxies for the PREs. The expected “true” values of the PREs are those calculated from the total unpaired electron spin density obtained from a DFT calculation (1380 and 1223 basis function for rubredoxin and plastocyanin, respectively). Thus, the effective distances, reff,DFT, calculated by integrating the contributions from the total unpaired electron spin density, ρα(r) – ρβ(r), obtained from the DFT calculation are compared to the effective distances calculated from the unpaired electron spin density of the metal–ligand pNBOs. The latter distances were obtained from the spin density matrix, PpNBO, corresponding to the metal–ligand pNBOs and lone-pairs of Table 1, using eq 3. In addition, the effective distances, reff,DFT, are compared to geometric metal–proton distances (Euclidean distances), rpd = ∥rH–Me∥, obtained from the X-ray structure (PDB code 1PLC(52) and PDB code 4RXN(51)). Finally, Figure 5b compares the effective distances calculated from the experimental R1p for the protons of spinach plastocyanin using eq 1 and a correlation time[54] τc,1 = (R1e,Cu(II) + τR–1)−1 = (4.1 ± 0.6) × 10–10 s at 18.8 T, with the pNBO and the geometric distances.

Figure 5

Evaluation of the effective distances, reff, between the protons and the unpaired electron spin of (a) the 144-atom model complex of plastocyanin,[53] and (c) the 104-atom model complex of rubredoxin.[11] In (a) and (c), the distances calculated from the total unpaired electron spin density obtained from the DFT calculation, reff,DFT, are compared to the corresponding effective distances calculated from the metal–ligand pNBOs of Table 1, (red ⊙) reff,pNBO, and the geometric distances obtained from the X-ray structures of the model complexes (blue ⊡), rpd = ∥rH–Me∥. In (b), the experimentally derived effective distances, reff,EXP, calculated from the experimentally measured R1p of spinach plastocyanin using eq 1 and a correlation time[54] τc,1 = (R1e,Cu(II) + τR–1)−1 = (4.1 ± 0.6) × 10–10 s, are compared to distances calculated from the pNBO density, (red ⊙) reff,pNBO, and with the geometric distances, rpd, of spinach plastocyanin[52] (blue ⊡). Vertical lines show the uncertainty calculated from the experimentally measured PREs.

Figure 5a,c shows that the metal–ligand pNBOs account well for the longitudinal PRE of all of the protons in both plastocyanin and rubredoxin. It is also noteworthy that, in the case of plastocyanin, the effective distances calculated from the four pNBOs of Table 1 agree extremely well with the experimentally measured distances, while discrepancies of up to 1 Å are seen for the point-dipole approximation (rpd vs reff,EXP = (R1p,EXP/ξ)−1/6). The main reason for the breakdown of the point-dipole approximation is the highly anisotropic electron spin delocalization in plastocyanin, where about 40% of the unpaired electron spin is located at the strongly bound cysteine sulfur;[9] yet the pNBOs account for this delocalization.

For rubredoxin, the effective distances derived on the basis of the point-dipole description deviate only about 0.5 Å from those derived from the DFT density (Figure 5c). This is because the spin delocalization in rubredoxin is smaller and more isotropic than in plastocyanin. Thus, only 25% of the unpaired electron spin is located at the ligand sulfur atoms in rubredoxin, and the sulfur ligands are arranged in a symmetrical tetrahedron around the ferric ion. Still, the pNBO description of the unpaired electron spin provides a significantly better agreement with the full DFT description than the point-dipole description.

An important property of the metal–ligand pNBOs is the transferability; that is, the metal–ligand pNBOs of the large complexes are very similar to the metal–ligand pNBOs calculated for the small complexes. As a comparison of the transferability, we have transferred the metal–ligand pNBOs from the small model complexes onto the large model complexes and subsequently calculated effective proton–electron distances in the large complexes. A comparison of the resulting effective distances, which is detailed in Supporting Information, shows that effective proton–electron distances that are within 0.25 Å are obtained by using the pNBOs from a small complex.

In conclusion, the metal–ligand pNBOs electron density allows accurate calculation of the longitudinal PRE for all protons of the two metalloproteins, rubredoxin and plastocyanin. Thus, effective distances are obtained that are within 0.25 Å (0.084 Å rmsd) of the full DFT description for plastocyanin and within 0.10 Å (0.068 Å rmsd) for rubredoxin. This implies that the longitudinal PREs of all protons in a metalloprotein can be included as constraints in a molecular dynamics simulation, provided the unpaired electron is described by the metal–ligand pNBOs or by a more sophisticated model.

Paramagnetic Relaxation of 15N Nitrogen

The paramagnetic relaxation enhancement of heavy atoms such as 13C and 15N is more complicated to interpret in terms of unpaired electron spin density and geometric distances than the relaxation of 1H nuclei discussed above. Even so, the 15N data can often provide complementary information in addition to that obtained from the 1H relaxation. The smaller gyromagnetic ratio of 15N as compared to 1H results in sharper 15N NMR signals that can give information about paramagnetic metal-sites in cases where proton data are unavailable due to extensive line broadening and signal overlap. Thus, PREs of heavy atoms are particularly valuable in studies of high-spin metalloproteins, where the PREs and line broadenings are large and the NMR signals of protons in the vicinity of the metal–site are often broadened beyond detection. Below, it is demonstrated how structural information can be derived from the relaxation rates and the chemical shifts of 15N nuclei using NBOs to describe the electron spin density.

In contrast to the proton PREs (Figure 5), the amide 15N nuclei PREs cannot be calculated from the metal–ligand orbitals alone, as shown for the rubredoxin model in Figure 6a. This is in accordance with previous observations that the relaxation of 15N nuclei is affected also by unpaired electron spin density in the local 2p orbitals.[11,42,43,55,56] In general, a very small amount of the total spin density resides in the local 15N orbitals (less than 0.03% per amide nitrogen of the rubredoxin complex). Yet, these small spin densities have substantial contributions to the paramagnetic relaxation because the local spin density is located in the immediate vicinity of the 15N nucleus. Therefore, to analyze 15N PREs of 15N nuclei in the immediate vicinity of the metal-site, the contribution from the local orbitals must be included.

Figure 6

Evaluation of the effective electron–15N distances, reff, obtained from the paramagnetic relaxation enhancement of 15N nuclei of the 104-atom rubredoxin complex calculated from (a) the Fe–S pNBOs specified in Table 1 and (b) the Fe–S pNBOs and eight local 15N orbitals. In both cases, the distances are compared to the effective distances, reff,DFT, calculated from the full DFT spin density (y-axis) taken as the “true” density.

In Figure 6b, the effective distances, reff,DFT = R1p–1/6/ξ1/6, calculated from the full DFT density are compared to effective distances calculated from a density that includes contributions from eight local orbitals as well as from the metal–ligand orbitals. For nonproline amide nitrogens, the eight local NBOs are LNBOs = {NBO(N–H), NBO*(N–H), NBO(N–Cα), NBO*(N–Cα), NBO(N–CO), NBO*(N–CO), LP(N), CR(N)}. The local antibonding NBOs, NBO*, are shown in Figure 7. For other nitrogen atoms (proline, histidine side-chain, etc.), the set of eight local orbitals consists of the lone-pair (LP), the core orbital (CR), and the NBO and NBO* to the three neighboring atoms. Figure 6b shows a good agreement between the calculated distances, which implies that the PRE of 15N can be derived from a small set of spin densities and basis functions that consists of the metal–ligand pNBOs and a few local orbitals. A similar comparison for the plastocyanin model (Figure S3) leads to the same conclusion. Overall, it is clear from the two widely different model complexes studied here that the unpaired electron spin density, which causes the relaxation of the 15N nuclei, must include the local 15N orbitals in addition to the metal–ligand orbitals.

Figure 7

Representation of local 15N orbitals, (a) the antibonding N–Cα natural bond orbital, (b) the combined core (CR) and lone-pair (LP) orbital (see text; eq 6), and (c) the antibonding N–H natural bond orbital.

Contribution from the Local Orbitals to the Fermi Contact Shift

Two experimental paramagnetic quantities can be obtained for each nucleus near a paramagnetic site, the PRE, which translates into the effective distance, reff, and the Fermi contact shift that reports on the spin density at the position of the nucleus. For protons, it is common to use the experimental relaxation parameters (R1p, reff) to derive either geometric distance constraints[28,30,41,57−61] or information about the electronic structure.[14,53] For 15N, the situation is more complicated as shown above because the 15N paramagnetic relaxation depends on the electron spin density in both the metal–ligand orbitals and the local orbitals.

As shown below, by using the 15N PRE and the 15N Fermi contact shift in combination, information on the geometric structure and geometric (metal-site)–15N distances can be obtained. In particular, we seek a single local 15N orbital, ψlocal, and its spin density, ρlocal, that together with the PRE permit the calculation of (metal-site)–15N distance restraints. To that end, we first show that the Fermi contact shift can be calculated from the density of the local natural bond orbitals, LNBOs, combined with the metal–ligand pNBOs. In particular, we show later that all of the local spin density can be projected onto one local orbital that allows distance constraints to be obtained from the PREs.

In Figure 8, the Fermi contact shift for 15N nuclei of plastocyanin and rubredoxin is calculated using two different spin densities;: the full DFT density (x-axis) and the density formed by the LNBOs and the metal–ligand pNBOs (y-axis). The δDFT shifts was obtained from eq 5 by summation over the complete DFT spin density. As above, we use the δDFT shift as the “true” experimental shift. The δNBO{Me–Ligand,local} shifts (y-axis) were obtained by including only the metal–ligand pNBOs and the LNBOs. The good correlations of Figure 8 show that the small set of orbitals consisting of the LNBOs and the metal–ligand pNBOs together with their spin densities describe accurately the unpaired electron spin density at the nitrogen position. Conversely, the spin density of the local orbitals can be obtained if the Fermi contact shift is known experimentally.

Figure 8

Fermi contact shift for 15N nuclei of (a) the 144 atom plastocyanin model complex and (b) the 104 atom rubredoxin model complex calculated from the full DFT density (x-axis) and from the electron spin density formed by the metal–ligand NBOs and eight local 15N orbitals. The excellent correlations show that the experimentally observed Fermi contact shifts (represented by δDFT, see text) can be calculated from a small number of orbitals, that is, the metal–ligand NBOs and a few local orbitals.

Describing the Local Spin Density by NBOs

Previous studies have shown that the natural bond orbitals resemble the chemical bonding environment.[31] Also, as shown above and illustrated in Figure 6b, Figure S3b, and Figure 8, the unpaired electron spin density of a paramagnetic metal-site can be represented accurately by a few (mainly antibond) natural bond orbitals. It is therefore tempting to hypothesize that the unpaired electron spin density of the local 15N orbitals can be described by one antibonding local NBO, ψlocal, that relates to the spin-polarization/spin-delocalization pathway from the metal-site to the 15N nucleus.

Many of the amide protons of rubredoxin are hydrogen bonded to one of the ligand sulfurs including the amide protons of Val8, Cys9, Tyr11, Leu41, Cys42, and Val44.[62] We therefore propose that the NBO*(N–H) antibonding orbital is the main carrier of the local unpaired electron spin in these cases. Similarly, the amide nitrogen of the ligand-residue Cys84 of plastocyanin is only three bonds from the sulfur that carries approximately 40% of the unpaired electron spin density. A very likely spin-polarization pathway is therefore through the α-carbon (Cu → Sγ→ Cβ→ Cα→ 15N) in which case the NBO*(N–Cα) is the main carrier of electron spin density into the local nitrogen orbitals. For some nitrogens, there could be multiple options for the spin-polarization path, for example, His37 Nε of plastocyanin or the Cys42 NH of rubredoxin. For His37 Nε, one could imagine a spin-polarization path both through Cε or through Cδ, while for Cys42 NH of rubredoxin, the spin-polarization can come from the hydrogen bond or through the Cα. In all cases, we find that the dominant spin-polarization path is the one that involves the lowest number of bonds between the nucleus and the metal-site, thus, Cε→ Nε for His37 and HN→ NH for Cys42. Please see Table S1 for all of the assigned spin-polarization pathways.

Deriving Effective (Metal-Site)–15N Distances from the Fermi Contact Shift and the PRE

A key property of the local 15N orbitals, ψlocal, is the ratio of p-orbitals (only relaxation effect) and s-orbitals (only Fermi contact shift). For residues where a distinct spin polarization pathway can be identified, the inherent property of the NBOs ensures that the ratio of p-orbitals and s-orbitals is correct. In these cases, the local NBO will “bridge” the Fermi contact shift and the local relaxation effect. For 15N, where no distinct “through-bond” spin polarization pathway can be identified, we propose that the unpaired spin density migrates to the 15N nucleus through its lone-pair (LP), which, in turn, polarizes an s-orbital to give a Fermi contact shift.[40] Thus, we propose that if no “through-bond” spin polarization pathway can be identified, which will often be the case when the 15N nuclei that are not hydrogen-bonded to the metal-site or in a ligand residue, the local unpaired electron spin density is given by a mixture of the lone-pair (LP; p-orbitals) and a core-orbital (CR; s-orbitals):where θ is a coefficient that describes the mixing between the CR and the LP orbitals. Moreover, it is shown in the Supporting Information that θ of eq 6 is to a very good approximation 86°. This value of θ is valid for both plastocyanin and rubredoxin that have very different metal-sites, that is, different paramagnetic metal ion, different spin, and different coordination sphere. Therefore, we propose that the angle θ ≃ 86° is general applicable for calculating 15N relaxation of other paramagnetic systems. Knowing the local orbital and its spin density allows us to calculate the effective (metal-site)–15N distances according to eq S6:where ψM is the metal–ligand pNBOs, ψlocal is the local 15N orbital that accounts for both the Fermi contact shift and the PRE, while ρlocal is the spin-density of the local orbital, and ρM is the metal–ligand spin-density; ρM ≃ 1.

Validation of the Applied NBO Approach Using Experimental 15N PREs and Fermi Contact Shifts

We have tested the derived NBO approach by investigating how well the experimentally measured PREs and Fermi contact shift agree with geometric structure information obtained from the crystal structure of rubredoxin (5RXN). Figure 9 compares the effective distances, reff(r′,ψlocal,δcon,EXP), calculated from the position, r′, of the 15N nuclei relative to the Fe–S4 site of rubredoxin and the experimentally measured Fermi contact shift, δcon,EXP, using eq 7, with the effective distances, R1p,EXP–1/6, obtained from the experimentally measured paramagnetic relaxation enhancements. The Fermi contact shifts, δcon,EXP, were those measured previously from 15N labeled rubredoxin.[45] The local spin densities, ρlocal, were derived from the experimental shifts, δcon,EXP, as described in the Supporting Information. For residues with a through-bond spin-polarization pathway the local orbital, ψlocal, was derived from that NBO, while a well-defined mixture of a lone-pair and a core orbital, eq 6, and θ = 86° was used for residues without a spin-polarization pathway. The excellent correlation between reff(r′,ψlocal,δcon,EXP) and R1p,EXP–1/6 in Figure 9 corresponding to a straight line with slope ξ–1/6 (see eq 1) shows that distance information with an rmsd of 0.15 Å can be obtained from 15N longitudinal PREs, provided that the Fermi contact shift is obtained. In contrast, the correlation of R1p,EXP–1/6 with distances obtained from the point-dipole approximation (the geometric distance, rpd, between 15N and Fe) is quite low with a Pearson’s correlation coefficient of R2 = 0.40, which clearly demonstrates the improvement achieved by using the scheme detailed above.

Figure 9

Comparison of effective electron–15N distances and experimentally measured longitudinal PREs for oxidized rubredoxin. The red circles (○) are the effective electron–15N distances, reff(r′,ψlocal,δcon,EXP), calculated from the metal–ligand pNBOs, the experimentally measured Fermi contact shift, δcon, and the position, r′, in the crystal structure. The blue squares (□) are the electron–15N distances calculated within the point dipole approximation; that is, the distance is given by the geometric distance between the nitrogen and the iron, rpd, and compared to the longitudinal PREs. Vertical bars represent uncertainties derived from the uncertainties of the experimentally measured PREs. The red (blue) dashed line is a best fit line, y = ax, to the reff(r′Xray,ψlocal,δcon,EXP) (rpd) correlations with R1p,EXP–1/6.

Structural information is now directly available from eq 7 via r′, because the spin density ρlocal (calculated from δcon,EXP) and the form of the local orbital ψlocal can be obtained for all 15N nuclei. The allowed positions for a 15N nucleus are those r that solve eq 7, with reff–6(r,ψlocal,δcon,EXP) = R1p,EXP/ξ. This is illustrated in Figure 10a for 15N of Pro40 of rubredoxin, where the red surface represents the positions r that solve eq 7 (ξ = 3000 ms–1 Å–6 for high-spin rubredoxin).[11] Each point on the surface, that is, for each r′ of reff(r′,ψlocal,δcon,EXP), took approximately 25 s to calculate on a standard Linux computer. For comparison the solution of the simple point-dipole approximation, rpd–6 = ∥rFe–N∥–6 = R1p,EXP/ξ is shown in Figure 10b. The NBO surface and thus the allowed positions, r, obtained from the solution of eq 7 are in excellent agreement with the crystal structure because the surface goes through the 15N of Pro40. In contrast, positions derived from the point-dipole approximation shown in Figure 10b are at variance with the crystal structure, because the surface does not pass through the 15N. The distance between the “geometric surface” and the nitrogen in Figure 10b is 0.9 Å.

Figure 10

Surface representation of allowed positions for 15N of Pro40 of rubredoxin, (a) derived from the experimental Fermi contact shift, the local orbital of eq 6 and the experimentally measured longitudinal relaxation rate, R1p,EXP, using eq 7. (b) Allowed position derived using the point-dipole approximation, ∥rFe–N∥–6 = R1p,EXP/ξ. Clearly, the surface derived using the NBO approach described above is in excellent agreement with the crystal structure, while the surface derived using the point-dipole approximation is at variance with the crystal structure.

All taken together, the unpaired electron spin density of the local orbitals of the 15N nuclei can to a good approximation be described by one local orbital, ψlocal, where the density ρlocal is derived from the Fermi contact shift. Thus, the metal–ligand pNBOs together with ψlocal and δcon,EXP can be combined with the experimental PRE to derive allowed positions of the nitrogen. Subsequently, the experimental longitudinal relaxation rate and Fermi contact shift can be used as constraints in structure determination of paramagnetic metalloproteins.

Molecular Dynamics Simulations and Structure Optimization Using the NBO Approach

Restrained molecular dynamics simulations has become an important tool for determining allowed conformations of protein structures and for sampling protein dynamics. These simulations and structure optimizations have previously been restrained with experimental data, such as chemical shifts,[63−67] order parameters,[68,69] residual dipole couplings,[67,70] proton pseudocontact shifts[71−73] and proton PREs,[74−76] and has provided valuable information about protein structure and dynamics. Here, we include the PRE and the Fermi contact shift of the 15N of rubredoxin as restraints in a molecular dynamics simulation and a structure optimization using the NBO approach described above.

As detailed in the Supporting Information, we generated starting structures that were sufficiently different from the crystal structures to test our approach. These initial structures were generated by a molecular dynamics simulation at 373 K for 10 ns, and only the five metal–ligand atoms, FeS4, were restrained by a set of distance restraints obtained from the crystal structure. The backbone rmsd of residues 6–11 and 39–44 that surround the metal-site to the crystal structure was approximately 1.5 Å for these starting structures. Only the rmsd of residues 6–11 and 39–44 are considered here, because these are the only residues where paramagnetic NMR data are available. To check the ability of the MD force field to generate structures that are in agreement with the crystal structure, we performed a 500 ns simulation at 298 K where only the FeS4 atoms were restrained (Figure S5). During the 500 ns, the rmsd of residues 6–11 and 39–44 stayed at approximately 1.5 Å, thus showing that the MD force field on its own is not capable of bringing the structure to a state that is in agreement with the crystal structure. Consequently, we need to restrain the simulation with experimental data to generate structures that are in agreement with the crystal structure.

The paramagnetic NMR parameters of rubredoxin were included using eq 7 to calculate effective (metal-site)–15N distances that were converted to distance restraints (see the Supporting Information). The effective distances were approximated with distance restraints and were updated every 100 ps of the simulation. Thus, every 100 ps of the simulation, the metal–ligand pNBOs and the local 15N orbitals were transferred from the initial DFT calculation of the model complex to the current structure of the simulation. The effective distance, reff(r′,ψlocal,δcon,EXP), was calculated at three positions along the Fe–15N vector using eq 7 as shown in the Supporting Information. Finally, a distance restraint was created at the distance where the experimental effective distance, (R1p,EXP/ξ)−1/6, agrees with the calculated (see Figure S6).

The structure of the residues around the metal-site approaches that of the crystal structure (rmsd of residues 6–11 and 39–44 is 0.55 ± 0.1 Å), Figure 11, when the molecular dynamics simulation is run at 298 K and is restrained by the paramagnetic restraints. Moreover, when the temperature of the simulation is lowered to 50 K, the structures becomes even more similar to the crystal structure (rmsd of residues 6–11 and 39–44 is 0.40 ± 0.03 Å; average three runs). This shows that a well-defined and accurate structure is obtained when the paramagnetic NMR parameters are included as restraints in simulations that use the NBO approach described above. To test the strength of the NBO approach as compared to the traditional point-dipole approximation, we performed two molecular dynamics simulations where the 15N PREs were included using the point-dipole approximation (blue squares of Figure 9). During the two 10 ns simulations at 298 K, the rmsd of residues 6–11 and 39–44 to the crystal structure was 0.65 ± 0.05 Å (Figure S7). When the temperature was lowered to 50 K, the rmsd dropped to 0.63 ± 0.02 Å (average of two).

Figure 11

Restrained molecular dynamics simulations of oxidized Fe(III) rubredoxin. The 15N PREs and Fermi contact shifts of residues 6–11 and 39–44 were included as geometric restraints in the simulations. (a) Backbone rmsd calculated for residues 6–11 and 39–44 between the crystal structure and the structure during the molecular dynamics simulation. During the 373 K simulation, the rmsd reaches ca. 0.6 Å within the first few nanoseconds, cooling to 298 K results in a slightly better agreement between the structures and the crystal structure (rmsd ≃ 0.55 Å), while cooling to 50 K results in an excellent agreement (inset). (b) The structures at 10, 11, ..., 19 ns (green) overlaid with the crystal structure (magenta). Residues where paramagnetic NMR data are available (residues 6–11 and 39–44) are shown with bright colors, whereas other residues are shown with pale colors. The sphere is the Fe(III) ion. (c) The final structure, at 22 ns, obtained after cooling to 50 K and overlaid with the crystal structure (5RXN) of rubredoxin.

In summary, including the 15N PRE and Fermi contact shift as restraints in a molecular dynamics simulation ensures that the structure samples states that are in accordance with the crystal structure. Moreover, the obtained structures are in excellent agreement (rmsd ≃ 0.40 Å) with the crystal structures after lowering the temperature of the simulation. Restraining the molecular dynamics simulation with the 15N PREs in the point-dipole approximation improves the simulation as compared to the unbiased simulation; yet, the backbone rmsd of the residues in the vicinity of the metal-site is significantly larger than that of the simulation restrained with the NBO approach.

Conclusion

We have shown that the natural bond orbitals (NBOs) form a simple and natural basis for the unpaired electron spin density of paramagnetic metalloproteins, which marks a new area of the use of NBOs to describe electronic structures. The majority of the unpaired electron spin density can be described by a small number of metal–ligand NBOs that are closely related to the chemical bonding environment of the metal-site. Thus, for the model complexes of the blue copper protein plastocyanin and the iron–sulfur protein rubredoxin, about 90% of the unpaired electron spin density is described by the two-center metal–ligand natural bond orbitals. The nuclear paramagnetic relaxation rates of protons in these metalloproteins are calculated accurately by the two center metal–ligand pNBOs, even for protons within a few chemical bonds from the metal-site. This holds, because the paramagnetic relaxation of the protons is affected by the delocalization of the unpaired spin onto the ligand atoms, but is unaffected by unpaired electron spin density in local orbitals. Thus, the effect on the proton relaxation caused by the electron spin delocalization can be accounted for by the metal–ligand natural bond orbitals alone.

The paramagnetic relaxation of 15N nuclei in metalloproteins is affected also by unpaired electron spin in the local orbitals. We have shown that also the paramagnetic relaxation and Fermi contact shift of the 15N nuclei can be described within the NBO formalism, despite a more complex dependence of the relaxation on the unpaired electron. A local orbital, ψlocal, that relates to the spin-polarization/spin-delocalization pathway, together with the metal–ligand pNBOs, can account for both the nuclear relaxation and the Fermi contact shift of the 15N nuclei. Conversely, the spin density of the local orbital can be determined from the observed Fermi contact shift and subsequently used to calculate the contribution to the PRE from the local orbitals. Our method allows PREs of 15N to be directly interpreted in terms of the position of the nucleus relative to the metal–site with an uncertainty of 0.15 Å for the data shown here.

Treatment of metal ions and in particular divalent metal ions is quite challenging for molecular dynamics force fields. Describing the metal ion in a metalloprotein by a quantum mechanics calculation and the rest of the protein with a molecular mechanics force field (QM-MM method) allows simulations of metalloproteins; however, these simulations are in general slow, and only a few nanoseconds can by obtained in best-case scenarios. The hybrid NBO approach that we developed allows paramagnetic PREs and Fermi contact shifts of 15N nuclei to be included as restraints in molecular dynamics simulations. We show that structures that are in good agreement with the crystal structure of rubredoxin are obtained when molecular dynamics simulations of rubredoxin are restrained with paramagnetic NMR parameters using the NBO approach. The NBO approach, which allows paramagnetic NMR data of 15N nuclei to be used as restraints in molecular dynamics simulations, adds to a growing list of experimental restraints in molecular dynamics simulations that are significantly impacting the utility of simulations to describe dynamics and allowed conformations of proteins.

Materials and Methods

Sample Preparation

Plastocyanin from spinach prepared and purified as described previously[77,78] was kindly supplied by Jens Ulstrup and Hans E. M. Christensen, the Technical University of Denmark. The protein was dissolved in 99.99% D2O at pH 7.0 (meter reading). All samples contained 100 mM NaCl.

NMR Experiments

The NMR experiments on plastocyanin were performed at 298.2 K and a 1H frequency of 800 MHz using a 1.5 mM fully oxidized spinach plastocyanin sample. The signal eliminating relaxation filter (SERF) experiment[79] was used for measuring the relaxation rates of the fast relaxing contact shifted signals, as described previously.[14] The sweep width used in the SERF experiments was 160 kHz, and the specific SERF parameters[79] were R1° = 50 s–1, R1s = 2 s–1, ta = 90 ms, and c = 1.

All 15N resonances of oxidized rubredoxin were assigned experimentally by reference to chemically synthesized rubredoxin samples that contain 15N labeling of only one type of amino acid or one specific amino acid.[45] The experimental assignments were verified using DFT calculations.[12,45] Longitudinal relaxation measurements of the hyperfine-shifted 15N signals of oxidized rubredoxin were collected on a Bruker DMX-500 NMR spectrometer by direct observation of the 15N spectrum[11] and were recorded at 283 K. The relaxation decays were monitored by an inversion–recovery method[80] with the 180° pulse replaced by a composite inversion pulse.

DFT Calculations of Rubredoxin

The structure of the Fe(SCH3)4– model complex shown in Figure 1A was constructed from the crystal structure of oxidized C. pasteurianum rubredoxin (PDB code 4RXN(51)). The density functional calculations were performed with the Gaussian 03 program[81] using the spin-unrestricted B3LYP hybrid density functional approach.[82] The basis set used for all atoms was the 6-31G* basis[83,84] (196 basis functions), which was implemented with the following keyword in the Gaussian input: “#p ub3lyp/6-31g* SCF=Tight Prop=(efg,epr)”. The positions of the hydrogen atoms were optimized by minimizing the total DFT energy, while the positions of the ferric iron and the four sulfur atoms were fixed. The NBO analysis was performed with the Gaussian program using the keyword “pop=nboread”. Transformation matrices from the atomic orbital basis to the pNHO, to the pNBO, and to the NBO basis were obtained with the NBO keylist “$nbo archive plot $end” in the Gaussian input. Here, the NBO keyword “plot” provides the transformation matrices from the atomic orbital basis to the following bases: pNAO, NAO, pNHO, NHO, pNBO, NBO, pNLMO, NLMO, and MO. Alternatively, the transformation matrices from the atomic orbital basis to the pNHO, to the pNBO, and to the NBO basis can be provided by the NBO keylist “$nbo aopnho aopnbo aonbo $end”.

The structure of the 104-atom model complex was identical to the model complex used previously[11,12] to study the hyperfine Fermi contact shifts and paramagnetic relaxation of oxidized rubredoxin. Thus, the model complex was constructed from the X-ray structure of oxidized C. pasteurianum rubredoxin (PDB code 4RXN(51)) and consisted of the ferric iron and the two hexapeptide chains, Cys6-Thr7-Val8-Cys9-Gly10-Tyr11 and Cys39-Pro40-Leu41-Cys42-Gly43-Val44. The residues Thr7, Val8, and Leu41 were converted into glycine residues, formyl groups were added to the N-terminus, while the C-terminal residues Tyr11 and Val44 were truncated to N-methyl groups. The density functional calculations were performed with the Gaussian 03 program[81] using the spin-unrestricted B3LYP hybrid density functional approach,[82] which allows the effects of spin polarization and electron correlation to be included for large model systems. The basis set used for all 104 atoms in the model complex was the 6-311G** basis[85] (1380 basis functions). Subsequently, the NBO analysis was performed with the NBO 5.0 program.[49]

DFT Calculations of Plastocyanin

The structure of the Cu(Im)2(SCH3)S(CH3)2+ model complex shown in Figure 1B was constructed from the crystal structure of plastocyanin from poplar leaves (PDB code 1PLC(52)). The density functional calculations were performed with the Gaussian 03 program[81] using the spin-unrestricted B3LYP hybrid density functional approach,[82] and the basis set used for all atoms was the 6-31G* basis[83,84] (286 basis functions). The positions of the hydrogen atoms were optimized by minimizing the total DFT energy, while the positions of the other atoms were fixed. The NBO analysis was performed as described for the Fe(SCH3)4– model complex.

The 144-atom model complex of plastocyanin was constructed from the crystal structure of plastocyanin from poplar leaves (PDB code 1PLC(52)). The model consisted of the cupric ion and two peptide chains, Pro36-His37-Asn38 and Tyr83-Cys84-Ser85-Pro86-His87-Gln88-Gly89-Ala90-Gly91-Met92-Val93. The model was further truncated by substituting P36 and Y83 by formyl groups and N38 and V93 by N-methyl groups. Although Ser85 and Gln88 were substituted by alanine residues to reduce computational time, all charge stabilizing atoms were included in the model. The density functional calculations were performed with the program Gaussian 03 revision C.02[81] using the spin-unrestricted B3LYP hybrid density functional approach.[82] The basis set used for the sulfur and copper atoms was triple-ζ basis TZVP,[86,87] while the 6-31G* basis[83,84] was used for all other atoms (1223 basis functions). The positions of the hydrogen atoms, the carbonyl oxygen atoms, and the imidazole ring of H37 were optimized in the SCF calculation. The NBO analysis was performed as described for the Fe(SCH3)4– model complex above. The DFT calculations on the two protein model complexes were carried out using an Apple G5 Xserve cluster with 34 processors at the University of Copenhagen, and two SGI Altix 3300 servers each with 12 Intel Itanium 2 processors at the University of Wisconsin–Madison.

Calculation of Effective Distances

The effective distances were calculated from the matrix representation of the dipole operator and the Fock–Dirac spin density matrix, as described in eq 3. Here, the spin density matrix in the atomic orbital basis was obtained as a part of the output from the NBO analysis. The matrix representation of the dipole operator in the atomic orbital basis, F2ν, was calculated by in-house software. Here, the dipole tensor was calculated in the irreducible tensor representation,[42] ℱ̂2ν(r′) = ∥r′∥–3Y2ν(r′/∥r′∥) andwhere ϕ(r) and ϕ(r) are atomic orbitals and the integral is over the complete volume (V). The contribution to the integral from the singularity, r → r′, was accounted for by a Dirac-delta function, as described previously.[88]

Molecular Dynamics Simulations and Structure Optimization of Rubredoxin

All MD simulations were performed using the Gromacs 4.5.4[89] package with standard parameters and the AMBER99SB-ILDN[90] force field. The simulations were initiated from the crystal structure of rubredoxin, 5RXN, where protons were added. The protein was placed in a box (150 nm3) of explicit TIP3P water, and periodic boundary conditions were applied. Thirteen sodium ions were added randomly to neutralize the system. The structure was first energy minimized with 1000 steps of steepest decent gradient minimization, and subsequently the water and the proton positions were equilibrated by a 100 ps simulation at 298 K with positional restraints applied to all heavy atoms. For all of the simulations thereafter, the iron and sulfur atoms were restrained to the conformation of the crystal structure with a distance matrix restraint with force constant of 106 kJ/mol/nm2.

First, a 20 ns simulations at 373 K was carried out to create structures that are different from the crystal structure and structures on which we could test the performance of the paramagnetic NMR restraints for improving the structure. The total backbone rmsd after 20 ns to the crystal structure was 2.3 Å, and the backbone rmsd of the residues around the metal-site (residue 6–11 and residue 39–44) was 1.5 Å. The structures obtained after the 20 ns run serve as the starting point for the restrained molecular dynamics simulations and structure optimization.