Pseudocontact shifts from mobile spin labels

This paper presents a detailed analysis of the pseudocontact shift (PCS) field induced by a mobile spin label that is viewed as a probability density distribution with an associated effective magnetic susceptibility anisotropy. It is demonstrated that non-spherically-symmetric density can lead to significant deviations from the commonly used point dipole approximation for PCS. Analytical and numerical solutions are presented for the general partial differential equation that describes the non-point case. It is also demonstrated that it is possible, with some reasonable approximations, to reconstruct paramagnetic centre probability distributions from the experimental PCS data.


Introduction
Pseudocontact shift (PCS) is an additional contribution to the nuclear chemical shift caused by the presence of a paramagnetic centre in close proximity to the nucleus in question. 1,2 PCS is very well researched and is widely used as a source of structural restraints for paramagnetic metalloproteins. [3][4][5][6][7] Even when a protein is not naturally paramagnetic, the commonly occurring calcium, magnesium and zinc binding sites would usually coordinate a lanthanide well. 8 In combination with artificially introduced lanthanide-containing tags PCS is also used to determine relative orientations of protein domains. 9 In magnetic resonance imaging it is useful as a reporter for local pH and oxidation potential. 10,11 General equations describing chemical shielding, obtained by an assiduously systematic application of perturbation theory, are due to Ramsey. 12,13 The first paper dealing with a point dipole approximation for chemical shift was published by McConnell 14 in 1957 -he noted that shielding by sufficiently distant electrons could be expressed via an effective magnetic susceptibility tensor; for paramagnetic molecules this tensor is a function of the spin Hamiltonian parameters of the paramagnetic centre. 15 Analytical treatments for specific classes of d-and f-transition metal complexes using ligand field theory have been reported by Bleaney, 16 Golding, 17 and Stiles. [18][19][20][21] Modern quantum chemistry defines paramagnetic shift as the Frobenius inner product between the nuclear hyperfine coupling tensor and the magnetic susceptibility tensor. 22,23 Both parameters are difficult to compute because they have contributions from spin-orbit coupling and often require non-perturbative treatment of relativistic effects within multi-configurational ab initio methods, 24,25 as well as conformational averaging. Accurate quantum chemical calculations are therefore limited to a few dozen atoms.
Structural biologists couldn't care less -most nuclei in macromolecules are far enough away from the paramagnetic centre for McConnell's' approximation to be accurate 2,26 and the resulting formula for the pseudocontact shift produced by a point source with a given magnetic susceptibility anisotropy has been of great service to protein structure and dynamics research over the last 30 years. 8,27 Excellent software packages exist that make PCS analysis in proteins and nucleic acids straightforward and informative. 5,28 There remains one important unsolved problem: pseudocontact shift prediction and analysis in large systems that feature fast conformational mobility of the paramagnetic centre. For such systems the point dipole approximation is no longer valid at short distances and quantum chemical calculations are prohibitively expensive. For this reason, PCS measurements in close proximity to the tag are often excluded from the analysis because they are not expected to fit the point dipole formula. 29 In this communication we introduce an analytical approach based on the recently published partial differential equation for PCS 30 that views the paramagnetic centre as a probability density distribution in three dimensions. This approach clarifies the key features of that density that affect PCS. It may also be used to recover the distribution itself from the experimental PCS data.

Pseudocontact shift from a point paramagnetic centre
The point source formula for PCS is best derived using a classical physics argument. Placed in an external magnetic field 0 B , a paramagnetic centre would acquire the following magnetic dipole moment: where 0  is the vacuum permeability and χ is the magnetic susceptibility tensor. 26 This linear response assumption is valid for an ensemble of non-interacting paramagnetic centres when T e 0 kT  μ B  -true for most metalorganic systems at room temperature. For a point centre, the induced dipole creates the following magnetic field at the relative position r : where the dipolar matrix is: It is easy to demonstrate that the dipolar matrix only depends on the direction of the position vector   ,   r = and does not depend on its length r  r -this provides a clean separation of coordinates that will be useful below.
The change in the energy of a nuclear magnetic moment n μ produced by placing it at the position r relative to the paramagnetic centre would be The associated chemical shift tensor is the second derivative of E  with respect to n μ and 0 B : 12,31 The isotropic average of this tensor is the familiar point dipole expression for the pseudocontact shift: 14,26       T  point  point  3  3  T   1  1  1  1  Tr  Tr  Tr  3  4  4 Dr χ χ r r (6) that may also be rewritten via second-rank spherical harmonics: where m  are the irreducible spherical components of χ : It is important to note that the isotropic part and the first spherical rank component of the magnetic susceptibility tensor do not enter the equation for PCS. The five irreducible spherical tensor parameters in Eq. (8) may also be expressed as axiality and rhombicity, along with the three parameters (e.g. Euler angles) specifying the orientation of the principal axis frame.

Pseudocontact shift from a distributed paramagnetic centre
A less well explored situation is when the paramagnetic centre is distributed with some probability density   e  r within the molecular structure. In such a situation, the magnetic susceptibility tensor would also in general be position-dependent. Integration of the point PCS expression in Eq. (5) over the probability density produces the following expression for the effective dipolar shift tensor at position r : This integral is a convolution of the dipolar matrix D divided by the cube of the distance with the product of susceptibility tensor and probability density: The simplest way to proceed is to use the Fourier transform because convolution is equivalent to multiplication in the k-space, and the Fourier transform of the dipolar matrix is very simple: (11) where k is the angular part of the k-space vector k . Another useful property of the dipolar matrix is that the inverse Fourier transform of its product with a function in k-space can be expressed as an action of the differential operator in real space: where   g k is the Fourier transform of       g   r r r  and the differential operator S has the fol- This is a reciprocal equation to Eq. (3) -multiplication by a vector in Fourier space is equivalent to taking a gradient in real space and division by 2 r in real space is equivalent to the inverse Laplacian in k-space.
In the case where both the probability density and the magnetic susceptibility tensor are position-dependent, this operator acts on their product and the following general expression is obtained for the matrix elements of the paramagnetic shift tensor: The second term in the right hand side effectively subtracts the Fermi contact part of the full paramagnetic shift, ensuring that PCS does not depend on

 
Tr  in the same way as it happens in the point model.
Eq. (14) simplifies significantly under the assumption that the susceptibility tensor is the same at all locations. After the isotropic average is taken, the result is: In the Fourier space this equation can be written as where    k is the Fourier transform of    r and    k is the Fourier transform of    r .
In the derivation presented above,    r is the statistical probability density of the paramagnetic centre, but one can make an approximate parallel here with the quantum mechanical spin density. The dipolar part of the hyperfine coupling at the nuclear position r for a given spin density meaning that the dipolar hyperfine coupling tensor field, viewed as an integral over the spin density, can also be expressed using the differential operator from Eq. (13): It must be noted that Eq. (18) only accounts for the dipolar part of the hyperfine coupling and does not include the orbital contribution, which is important in the immediate vicinity of heavy ions. 32,33

Analytical solution to the direct problem
The "direct" problem will be defined here as the task of calculating PCS from a given magnetic susceptibility tensor, and a given probability distribution of the paramagnetic centre at the specified nuclear coordinates. General case with a position-dependent magnetic susceptibility tensor is described by Eq. (14). It is clear from the form of Eq. (14) that in order to disentangle effects of two spatial functions of density and susceptibility we need to know a priory at least one of those. Below we analyse the special case where magnetic susceptibility tensor is the same at every point of the probability density as described by Eq. (15) . This approximation is reasonable for modelling mobility of lanthanide tags as it was shown by Shishmarev and Otting with "two-hinged" approximation where orientation of susceptibility tensor is the same for each rotamer. 29

General solution
The easiest way to solve Eq. (15) analytically is to expand the probability density in spherical harmonics: (20) where   m l r  are radial functions serving as expansion coefficients in this angular function series: For a spherically isotropic probability density distribution, the sum in Eq. (20) only has one term with l=0 and m=0; for a spherically anisotropic density there would also be higher terms in the spherical harmonics expansion. The Fourier transform of the density in Eq. (20) leaves the angular part the same but the radial part is integrated with spherical Bessel functions of the first kind: We shall substitute the density from Eq. (16) into Eq. (22) and expand the products of spherical harmonics using Clebsch-Gordan coefficients. The result has three terms when 2 l  and two when 2 l  : With this substitution in place, Eq. (16) acquires the following form: To get the final answer we must take the inverse Fourier transform of Eq. (24). The spherical harmonics again remain the same, and the radial part is integrated with spherical Bessel functions. The double integrals that make an appearance are all straightforward. Changing the order of integration and integrating the product of two spherical Bessel function with respect to k gives: Taking everything together we obtain the general analytical solution for Eq. (15):  (28) There are three physically different contributions associated with the three integrals of the radial probability density functions in Eqs. (25)- (27). The first one corresponds to the integral of the density weighted with a monotonically decreasing function over the region outside the sphere of radius r : (29) For any nucleus positioned outside the bounding sphere of the paramagnetic centre probability density this contribution is zero. The second contribution is proportional to the probability density at the nucleus itself: Figure 1. A schematic representation of the paramagnetic centre probability density (grey shaded area) and regions where different terms in Equation (28) contribute to the resulting pseudocontact shift. Outside the bounding sphere (grey solid line) only P IN is non-zero; inside the bounding sphere an additional term P OUT becomes important; P 0 contributes only inside the density. At a sufficient distance from the bounding sphere (grey dashed line) the point paramagnetic centre approximation becomes valid.
We are in practice unlikely to be able to measure chemical shifts of the nuclei that are directly underneath the spin density due to their fast relaxation. We are thus left with the third and the most important contribution that is associated with the part of the density that is inside the sphere of radius r : If we assume that there is no paramagnetic centre density outside that sphere, this would allow us to extend the upper integration limit to infinity and the integrals then correspond to the multipole moments of the probability density of the paramagnetic centre: A schematic diagram of the entire argument is given in Figure 1. Outside the bounding sphere the final expression for PCS is: This solution tells us that PCS is sensitive to the multipole moments of the paramagnetic centre probability density distribution. In the case of isotropic distribution, the PCS is the same as from the point source -this follows from Eq. (33) when we put l to zero; it simplifies into Eq. (7). Therefore, only the anisotropy of the density makes a difference in the PCS compared to the point source. We shall therefore proceed to explore some simple anisotropic probability densities.
More detailed derivation of Eq. (28) and Eq. (33) is provided in Supporting Information.

PCS from a Gaussian paramagnetic centre distribution
where a is the standard deviation of the Gaussian. Figure 2 demonstrates that a significant difference between the point and the isotropic Gaussian case appears only for 3 r a  , which is not particularly interesting or dangerous because nuclei at such short distances from the paramagnetic centre are not normally visible in PCS experiments due to their rapid transverse relaxation. is parallel to the Z axis of the susceptibility tensor eigenframe. The contribution to PCS from the isotropic part of the density is on the left and the contribution from the correction from anisotropic part is on the right.
A rather more worrying situation emerges when we consider a realistically anisotropic Gaussian distribution of the paramagnetic centre: The case with a c  corresponds to an oblate ellipsoid and the case with a c  to a prolate one. The anisotropy results in the emergence of the next multipole in Eq. (32)     with the amplitude that depends on the difference of the squares of axial and equatorial sizes of the Gauss-  (38) involves fourth rank spherical harmonics and the fifth power of the distance. It is illustrated graphically in

Analytical solution to the inverse problem
Eq. (33) is linear with respect to the components of the susceptibility tensor and the multipole moments of the probability density. Both sets of parameters are therefore easy to extract by fitting a sufficiently large data set comprising nuclear coordinates and pseudocontact shifts. However, because experimental PCS measurements inside the bounding sphere of the paramagnetic centre density are not usually realistic, only the multipole moments (rather than the density itself) may be extracted in a well-defined way.
An illustration for 300 nuclei randomly placed around a lump of paramagnetic centre probability density is given in Figure 3: the density was formed by four isotropic Gaussians with a standard deviation of 0.5 Å placed randomly within a 3×3×3 Å cube; nuclei were scattered within a 1 Å thick spherical layer at varying distances from the origin; the axiality and the rhombicity of the susceptibility tensor were chosen randomly from the typical range reported in the literature: 34 = -0.45 Å 3 and = -0.1 Å 3 . PCS amplitudes in this system range from ±90 ppm (at 5 Å distance) to ±1 ppm at (20 Å distance).   33) in the case when the nuclei are scattered in the spherical layer between 4 Å and 5 Å away from the origin (a schematic is given in Figure 3). (Right) fractional back fit error in the multipole moments of the paramagnetic centre probability density distribution in Eq. (33) as a function of the radius of the spherical layer containing the nuclei (a schematic is given in Figure 3).
It is clear from Figure 3 that, at distances comparable to the size of the paramagnetic centre probability distribution, the point model breaks down and further terms are required in Eq. (33). Their beneficial effect is illustrated in Figure 4 -adding terms of higher spherical rank to the expansion dramatically reduces the PCS back-fit error in the vicinity of the bounding sphere.
Due to the steep distance dependence and the wobbly angle dependence of the higher multipole terms in Eq. (33), the accuracy with which these terms may be extracted also falls steeply as the nuclei are moved further away from the bounding sphere of the paramagnetic centre probability density. This is illustrated in Figure 4 -it is clear that multipole moments of spherical rank higher than 4 cannot be reliably extracted no matter how close the nuclei are to the bounding sphere. We would not therefore recommend taking Eq. (33) beyond the fourth spherical rank.

Numerical solution to the direct problem
If the probability density of the paramagnetic centre is discretised on a finite three-dimensional grid, two general numerical avenues become available for the solution of Eq. (15) -the finite-difference method and the Fourier transform method. The latter has the advantage of being fast, but the disadvantage of requiring periodic boundary conditions. This section explores both methods and comments on their relative merits.

Finite difference methods
A general algorithm for the generation of elementary finite difference operators of a given derivative and accuracy order on a given finite grid has been published by Fornberg. 35 For a rectangular grid with N points in X direction, M points in Y direction and K points in Z direction, finite difference matrix representations of the relevant derivative operators acting on the vectorisation of the spin density cube are: (39) and similarly for the other first and second derivatives. In these expressions,   The matrices required for the solution of Eq. (15) are In terms of these matrices, the solution may be written as 1   σ L Kρ (40) where ρ is the vectorization of the probability density cube on the chosen grid and σ is the vectorisation of the pseudocontact shift field on the same grid; here and below we take only the traceless part of the χ tensor. The final step is to project out the PCS values on the nuclei, for which the interpolation matrix N P , also published by Fornberg in the same paper, 35 may be used. The final expression is 1 N N   σ P L Kρ (41) where the action by the very sparse inverse Laplacian matrix on the Kρ vector is beset computed using an iterative solver, such as GMRES. 36 Here and below the 1  L symbol should be understood in that sensethe inverse Laplacian is never computed explicitly.
This method is very easy to set up -see, for example, the kpcs.m function supplied with Spinach library. 37 Its downside is unfavourable scaling: the dimension of the finite difference matrices involved is KMN and the scaling of sparse solvers is approximately quadratic in the matrix dimension, meaning that the numerical complexity of this method scales approximately as the sixth power of the grid size. On a contemporary computer, a solution for a 256×256×256 point grid takes about an hour.

Fourier transform methods
An alternative method for solving Eq.(15) takes advantage of the Fourier domain expression in Eq. (16), where Fourier transform can be taken numerically using the fast Fourier transform (FFT) algorithm. 38 Using the same interpolation operator as in Eq. (41), we get k χ k σ P r k k (42) where χ stands for the traceless part of the magnetic susceptibility tensor. This method is about the same as the finite difference method in terms of the implementation complexity, but it is much faster -the complexity scaling of the three-dimensional fast Fourier transform is is close to cubic scaling with respect to the grid size. This takes the simulation time into the region of seconds. This method is particularly fast on modern computing hardware because fast Fourier transform routines are available for GPGPU coprocessor cards (NVidia Tesla K40 cards were used in this work). The performance gain is shown in the left panel of Figure 5. The only significant problem with the FFT method is the periodic boundary condition -care must be taken to ensure that the images do not contribute significantly to the solution. Thankfully, PCS decays cubically with distance; it is sufficient for the distance from the paramagnetic centre to the cube boundary and the distance from the nucleus to the cube boundary to be about three times larger than the distance from the paramagnetic centre to the nucleus -a numerical example is given in Figure 5.

Numerical solution to the inverse problem
The "inverse problem" refers to the task of reconstructing the paramagnetic centre probability density from the experimental values of pseudocontact shifts and atomic coordinates. The task may be formulated as finding the probability density    r that minimises the following functional: where the least squares error is (44) and the standard Tikhonov regularisation term emphasizing smooth solutions is 39 where  is a user-specified regularisation parameter. Regularisation is only needed when the inverse problem is ill-posed, i.e. the number of points in the probability density grid exceeds the number of experimental PCS points and/or all experimental points are outside the bounding sphere of the density. In that case the recovered density is only an approximation that agrees with the experimental data and has the required smoothness at the same time -this point is discussed further in Section 8.

Finite difference method
In order to proceed with minimising         r with respect to    r , the state-of-the-art Newton-Raphson minimisers require analytical expressions for its first and second variation. 40 In a matrix representation on a discrete grid, these variations become the gradient vector and the Hessian matrix. The following vector relations are useful as a starting point.
If we use ρ to denote the vectorised form of a discrete representation of    r on a finite three-dimensional grid, the expressions for the gradient and the Hessian of the least squares error part of Eq. (43) are: where advantage was taken of the fact that L and K are symmetric matrices and that Similarly, the gradient and the Hessian for the Tikhonov regularisation term are: It is in practice unnecessary to consider every point of the grid to be a variable. The location of the paramagnetic centre is usually known approximately, and it is sufficient to only consider the points in the immediate vicinity -typically a 20×20×20 Å cube around the expected location.

Fourier transform method
About a factor of a hundred in performance may be gained (at the cost of worrying about the periodic boundary as described in Section 6.2) by using the Fourier solution in Eq. (42). The overall error functional is still the same as in Eq. (43), but the terms for the least squares error and the Tikhonov regularisation now involve three-dimensional fast Fourier transforms: The first variation and the action by the second variation of r on a probe function    r are: Figure 6. Paramagnetic centre probability density reconstruction for Tm 3+ labelled calbindin D9K (PDB code: 1IGV), 34,45 for which the point approximation is known to produce a poor fit (right panel, blue circles). Numerical probability density reconstruction (left panel) reveals that the thulium ion probability density has an elongated shape, likely as a consequence of local conformational mobility or multiple competing oxygen coordination sites. The red cube is the volume in which the probability density is allowed to vary during the Tikhonov regularised reconstruction process.
An example, shown in Figure 6, deals with the paramagnetic centre density distribution in the well characterized calbindin D 9K in which one of the two calcium ions has been replaced by a thulium ion. 45 Detailed coordinate and PCS data for this protein is available in the supplementary information of the paper by the Florence group. 34 The point paramagnetic centre model produces a rather unsatisfactory fit in the vicinity of the metal ( Figure 6, right panel). A reasonable explanation (proposed to the authors by Gottfried Otting and Thomas Huber) would be that fast exchange between metal binding sites could be present. We can now confirm that there is indeed an elongated distribution in the metal position ( Figure 6, left panel).
Simple and intuitive though the reconstructions of the kind shown in Figure 6 might appear, a few words of caution are in order. An important consequence of the asymptotic behaviour of the three components of the exact solution in Eq. (28) is that the fine internal details of the paramagnetic centre probability density distribution are expected to be lost if the PCS data is measured for the nuclei that are positioned outside the bounding sphere of the density. In practice this means that the shape of the bounding surface of the probability distribution is well reproduced by the reconstruction, but its internal details are not. The consequences are illustrated in Figure 7. The left diagram shows a simple test density composed of three Gaussian functions. The middle diagram shows the reconstruction performed using PCS data from the nuclei that are distributed randomly within a 10 Å cube around the origin, including inside the density.
The right diagram shows the reconstruction performed using PCS data from the nuclei that are placed outside the 5 Å sphere around the density. The latter case is what typically happens in paramagnetic NMR; it is clear that a reasonable likeness of the true probability density is recovered, but the internal details of the distribution are lost because PCS values at those points are not sensitive to the radial functions in the probability density expansion in Eq. (20). This is both the consequence of the structure of the multipolar expansion and the price to pay for a Tikhonov solution to what is in general an ill-posed problem.

Conclusions
Delocalisation of the paramagnetic centre creates deviations from the point dipole model for the pseudocontact shift. The analytical solution to the partial differential equation for PCS indicates that the key factor in such deviations is the angular anisotropy in the probability density -PCS from spherical paramagnetic centre distributions is identical to the point dipole PCS outside the bounding sphere of the density.
Multipoles with spherical rank l in the paramagnetic centre distribution create rank 2 l  multipoles in the resulting PCS field that decay with the distance as 3 1 l r  outside the bounding sphere of the density. Due to fast decay of these higher rank components, only terms up to spherical rank four can in practice be extracted by fitting experimental data. The behaviour of the PCS field inside the bounding sphere is more complicated and may be analysed numerically; however, experimental measurements in such close proximity to the paramagnetic centre are rarely possible.
Discretising the paramagnetic centre probability density and the PCS field on finite grids makes it possible to solve the partial differential equation for PCS numerically, using either finite difference operators or fast Fourier transforms. The ill-posed inverse problem of reconstructing the probability density from experimental PCS data can then be solved using Tikhonov regularisation under the assumption that the paramagnetic centre has the same effective magnetic susceptibility anisotropy at every point in the distribution.
The methods described in this paper are implemented in versions 1.8 and later of Spinach library. 37

Acknowledgments
This work was made possible by a grant from EPSRC (EP/N006895/1). Stimulating discussions with Thomas Huber, Alexander Karabanov, Claudio Luchinat, Gottfried Otting and David Parker are gratefully acknowledged.