A partial differential equation for pseudocontact shift

It is demonstrated that pseudocontact shift (PCS), viewed as a scalar or a tensor field in three dimensions, obeys an elliptic partial differential equation with a source term that depends on the Hessian of the unpaired electron probability density. The equation enables straightforward PCS prediction as well as analysis of experimental PCS data in systems with multiple and / or distributed unpaired electron centres.


Introduction
Pseudocontact shift (PCS) is an effective additional contribution to the nuclear chemical shift that arises in open-shell chemical systems due to partial polarization of the electron spin by the applied magnetic field [1]. The primary application of PCS is in structural biology, where it provides additional distance restraints for molecular structure determination [2,3]. Pseudocontact shift is different from the contact shift in that it does not require electron-nuclear overlap and propagates instead through the dipolar coupling [1].
It can be verified by direct inspection that the commonly used point electron dipole expression for the PCS [4] has a zero Laplacian everywhere except the origin: 2 where ax   and rh   are axiality and rhombicity of the electron magnetic susceptibility tensor χ and the nucleus is located at   , , x y z relative to the electron. This is to be expected -all classical electromagnetic phenomena must obey Maxwell's equations -but the singularity at the origin also suggests the possibility of an elliptic partial differential equation existing for the harder case of the PCS generated by a non-point electron probability density   in which the source term   , r  χ  is unknown and has so far resisted all derivation attempts: a direct calculation of the Laplacian of the convolution of Equation (1) with a finite unpaired electron density comes across singular integrals that cannot be regularized [5]. Yet the prize is tempting -elliptic partial differential equations are a classical topic in mathematics and a simple enough PDE would generalize all PCS analysis problems, improve data interpretation close to the unpaired electron location and also provide a way of measuring spin label probability distributions in biomolecular EPR experiments [6]. It would be convenient too -numerical PDE solvers have been available in standard software packages for a long time. In this communication we derive the equation and comment on some of its properties.

Hyperfine shift as a total energy derivative
To facilitate subsequent mathematics, and also for the sake of completeness, we provide in this section succinct derivations, using the relatively modern total energy derivative formalism [7], of the classical expressions for the various components of the hyperfine shift tensor [1,4]. For his- Placed in a magnetic field 0 B  , a point electron at the origin with a magnetic susceptibility tensor 1 χ  would acquire an average magnetic moment: The additional magnetic field 1 B  generated by this dipole at the point r  is: The change in energy E for a nuclear magnetic dipole N   located at that point would be: and therefore, the additional chemical shift tensor experienced by the nucleus, measured relative to the unperturbed conditions, would be: where the minus appears because of the relationship between chemical shielding and chemical shift [8,9]. The isotropic part of this tensor is: For a point electron located at e r  and a point nucleus located at N r  the final expression is: A neat derivation for the hyperfine shift can also be given in terms of the hyperfine coupling tensor A . The spin Hamiltonian for the hyperfine interaction is and therefore the corresponding energy change for the k-th nucleus in the system is: After using the same derivative expression for the chemical shielding [7], we obtain: A significant advantage of these equations over other physically equivalent formulations is that the excitation structure of the quantum chemistry part of the problem is hidden from the userboth the hyperfine tensor and the susceptibility tensor are effective quantities that already incorporate all of the formidable complexities of the electronic structure theory [7]. The derivations given above are classical, but the only assumption in Equation (10) is that 1 χ  (meaning also that magnetic hyperpolarizability terms in the electronic structure theory are negligible) and that the electron relaxes sufficiently rapidly for Equation (3) to always remain valid.
Equation (8), in its various forms, has done considerable service to the NMR community over the last forty years -naturally occurring calcium, magnesium and other metals in biological systems can often be substituted with lanthanides and pseudocontact shift then used to obtain distance restraints for structure determination purposes [2,3]. The point electron dipole assumption does, however, have its validity range -Equation (8) is not applicable in close proximity of the metal centre, most notably in lanthanide spin labels, and also in the cases where the electron probability density is broadly distributed within the molecular structure.

Derivation of the elliptic PDE
The source term   , r  χ  in Equation (2) for pseudocontact shift induced by a distributed electron is not currently known. The most straightforward derivation is to notice that: 3 T e e N e N e N e 1 1 3 r r r r r r r r r r r r and to note that the Laplacian of the reciprocal distance is: With these observations in place, we can conclude that: The convolution of this expression with a finite electron probability density distribution   e r   then yields (after dropping the N subscript on the nuclear coordinate vector N r  ): After abbreviating the Hessian of   we arrive at a neat final result: A similar derivation for the full 3×3 dipolar shift tensor in Equation (6) yields: If chemical shielding is considered instead of the chemical shift, the minus sign on the right hand side disappears. Simplicity of Equations (16) and (17) stands in sharp contrast with the unfathomable spherical harmonic expansions [10] generated by ab initio treatments using the ligand field theory.

Analytical and numerical solutions
General solution strategies for Equations (16) and (17) are identical to those of Poisson's equation [11]. The special case of the point electron in Equation (7) is recovered by observing that the Green's function of the Laplacian is: and taking its convolution with the right hand side of Equation (16) Analytical treatments of more general unpaired electron probability density distributions may be simplified significantly by noting that the Laplacian in Equations (16) and (17) where   , r  χ  is either the right hand side of Equation (16), or each of the nine matrix elements on the right hand side of Equation (17). This is the general analytical solution for Equation (16) and for each matrix element of Equation (17), but in practice the interpretation of multipole moments with 2 l  is difficult. We would argue here that numerical treatment is easier computationally and also more interpretable from the physical point of view because the fundamental quantity there is the unpaired electron probability density   r   . Sparse matrix representations of 3D Laplacians with Dirichlet, Neumann or periodic boundary conditions are readily available [13], and the numerical solution to Equations (16) and (17) requires a single sparse matrix-inverse-times-vector operation, for example: where L is a matrix representation of the 3D Laplacian with appropriate boundary conditions [13] and H χ denotes the index transformation that stretches   Tr   H χ , which is a three-dimensional array, into a column vector. Equation (21) has been implemented into the latest version of our Spinach library [14]. On a contemporary workstation using Matlab the example solution given in Figure 1 takes a few minutes to compute. The conclusion from Figure 1 is that, for lanthanide-containing spin labels, the accuracy of the point PCS model is very high, presumably due to the localized nature of the f orbitals, and it is mostly the contact shift that is making their interpretation difficult in practice.

Inverse problems
The most interesting possibility offered by Equations (16) and (17) where P is the matrix that projects out pseudocontact shifts at the nuclear locations, expt   is a vector of experimental PCS data and 1,2  are Tikhonov regularization parameters [15]. The uncommon choice of regularization operators (maximum entropy and minimum Laplacian norm) is dictated by two practical considerations: 1. Electron probability distribution is expected to be localized in at least one dimensioneven extended conjugated systems, such as porphyrin and carotene radicals, have electron spin densities that closely follow the bonding network. Maximum entropy regularization is known to favour highly local solutions unencumbered by baseline noise [16] A synthetic example of computing the PCS field generated by multiple paramagnetic centres and then recovering their distribution from PCS data is given in Figure 2. Additional constraints on the probability density are non-negativity, fixed integral and zero boundary conditions: where N is the number of unpaired electrons in the system. Because the error functional in Equation (22) has two regularization parameters, a generalization of the L-curve method [18] to surfaces [19] is used here. Better regularization methods for Equation (16) that could improve the fidelity of the reconstruction in Figure 2C are undoubtedly possible, but are beyond the scope of the present work.
Attention should also be drawn again to the simple general connection between hyperfine shift and hyperfine coupling provided by Equations (10) -modern electronic structure theory packages are able to compute both hyperfine tensors and magnetic susceptibility tensors, meaning that hyperfine shift tensors may be obtained essentially for free after standard magnetic property runs in ADF [21], ORCA [22] or Gaussian [23], subject only to the electron spin relaxation time being sufficiently short for the approximation made in Equation (3) to be valid.
The source code for all examples provided above, as well as numerical infrastructure functions (3D finite difference operators, 3D interpolation operators, Tikhonov solvers, volumetric scalar field visualizer, etc.), are available in versions 1.5 and higher of Spinach library [14]. In all cases, the molecular geometry was optimized and the electron spin density estimated using DFT UB3LYP method [25,26] in vacuum with cc-pVTZ basis set [27] on light atoms and Stuttgart ECP basis set on europium [28].

Figure captions
CSGT DFT UB3LYP [29] method with the same combination of basis sets was used to estimate the magnetic susceptibility tensor. The points refer to the symmetryunique atoms in the ligand, excluding nitrogens that have significant contact shifts.
Simulation source code is available within the example set of Spinach library version 1.5 and later [14]. plot of the electron probability distribution obtained by solving the inverse problem as described in the main text. Pseudocontact shift was sampled at 500 random points emulating nuclear locations within the volume and fed into Equation (22), which was then minimized from a random initial guess. Simulation source code is available within the example set of Spinach library version 1.5 and later [14].