Model-free extraction of spin label position distributions from pseudocontact shift data† †Electronic supplementary information (ESI) available. See DOI: 10.1039/c6sc03736d Click here for additional data file.

Not a point, but a cloud: advanced PCS data analysis using 3D probability density reconstruction provides more information.


Introduction
Pseudocontact shi (PCS) is an additional chemical shi caused by the presence of a rapidly relaxing paramagnetic centre near the nucleus. 1,2 PCS is well understood theoretically [3][4][5][6][7] and is widely employed as a source of structural restraints in metalloproteins, [8][9][10] where commonly occurring Ca 2+ , Mg 2+ , Mn 2+ and Zn 2+ binding sites can oen coordinate a lanthanide ion instead. 11 A paramagnetic centre may also be introduced articially by attaching a lanthanide ligand tag to the protein surface. 12,13 The subject has a long-standing problemlanthanide-containing protein tags have signicant conformational mobility. 14 Even DOTA-M8, 15 which uses a sterically overcrowdedand therefore rigidmetal cage, 16 still has a exible linker. The conformational mobility of lanthanide tags is visible in the distance distributions measured by double electron resonance, 17 and in molecular dynamics simulations. 18 In this situation the commonly used point paramagnetic centre approximation 3,19 for PCS is not expected to be valid, 14,20,21 but quantum chemical calculations 6,7 are prohibitively expensive.
The problem is solved by the recently discovered partial differential equation that treats the probability density of the metal ion and the resulting pseudocontact shi as scalar elds in three dimensions. 20 We demonstrate here that it may be used to recover the spin label position distribution from the experimental PCS data. This creates a new window into protein structure and dynamics.

Extracting spin label distributions from PCS data
The point paramagnetic centre approximation for the pseudocontact shi 3,19 s(r) experienced by a nucleus at the position r relative to the metal has the following form: where c is the traceless part of the magnetic susceptibility tensor. If the paramagnetic centre is distributed with some probability density r(r), the convolution of this density with eqn (1) obeys 21 the partial differential equation that we have recently derived: 20 where V is the gradient operator, H r(r) is the Hessian of r(r) and the denominator is to be understood as the inverse Laplacian. The best way to solve this equation is to use fast Fourier transforms: 21 The task of recovering r(r) from point measurements of s(r) at the nuclei may be formulated as nding the paramagnetic centre probability density that minimises the following functional: where E[r] is the least squares error relative to the experimentally measured pseudocontact shis at the nuclear locations sampled by the operator P N , and T[r] is a Tikhonov regularisation term emphasizing smooth solutions with l being the regularisation parameter selected using the L-curve method. 22 State-of-the art numerical optimisation algorithms 23 require rst and second variations of the error functional with respect to the probability density. The rst variations are: and the actions by the second variations on a probe function h(r) are: The trust region reective Newton-Raphson minimiser, as implemented in the Optimisation Toolbox 24 supplied with Matlab, was used to obtain the optimum paramagnetic centre probability density on a nite grid, subject to the non-negativity constraint. Numerical implementation details are discussed in our recent paper 21 and the associated Matlab source code is available in the paramagnetic NMR module supplied with versions 1.8 and later of the Spinach library. 25

Protein preparation
The pACA plasmid used for the production of human carbonic anhydrase II (hCA-II) mutants was a generous gi from Carol A. Fierke (University of Michigan). 26 Double and triple mutants were prepared by sequential site-directed mutagenesis (performed as recommended by Zheng et al. 27 ). All constructs were expressed in uniformly 15 N labelled form and in selectively 15 N-Leu labelled form; all expressions were carried out in BL21(DE3) pLysS competent cells using standard methods. Conjugation with Lu 3+ , Tm 3+ or Gd 3+ containing DOTA-M8 tags was performed as previously described. 15
Pseudocontact shis were obtained by comparing 1 H-15 N HSQC spectra of diamagnetic Lu 3+ and paramagnetic Tm 3+ DOTA-M8 tagged mutants. PCS assignment was performed in two stages. At the rst stage, at least eight of the 26 leucine peaks in the selectively 15 N-Leu labelled mutants were identied manually and supplied to NUMBAT, 28 along with 15 N positions from the X-ray structure (PDB:3KS3) of hCA-II. 29 This enabled the identication of all leucine 15 N atoms and the extraction of the approximate point model parameters (metal position and the anisotropic part of the magnetic susceptibility tensor). At the second stage, these parameters were used in another round of NUMBAT calculations on the uniformly 15 N labelled mutants to assist in the identication of the rest of the shied signals, yielding a total of 364 (S217C), 366 (S50C, S220C) and 397 (S166C) unambiguous 1 H and 15 N PCS assignments. The relevant data is included into the ESI. †

Electron paramagnetic resonance
Pulsed EPR measurements were performed using a home-built Q-band EPR spectrometer 30 equipped with a rectangular broadband resonator that can accommodate oversized samples. 31,32 EPR samples of the spin-labelled hCA-II were prepared as a mixture of equal volumes of glycerol and the protein solution in an aqueous phosphate buffer at pH 6.8; the resulting protein concentration was approximately 100 mM. 50 mL of this mixture was transferred into a thin-walled quartz tube with 2.9 mm outer diameter, ash-frozen by immersion into liquid nitrogen and stored at À80 C.
Spin label distance distributions were measured at 10 K using the four-pulse DEER sequence. 33 All pulses were 12 ns long; the frequency offset between the pump pulse and the detection pulse was 300 MHz. The duration of all DEER traces was at least 3.0 mslong enough to sample and subtract the intermolecular background. Distance distributions were extracted using the DeerAnalysis package. 34 The optimal values of the Tikhonov regularization parameter were found using the L-curve method. 22

Rotamer library details
Position distributions for the lanthanide ion enclosed in the DOTA-M8 chelating tag 15 were also estimated by constructing, for each of the four tagging sites (S50C, S166C, S217C, S220C), a complete set of distinct energetically accessible rotamers. The set was built by Monte-Carlo sampling of the dihedral angles present in the tag, followed by agglomerative clustering using the following distance metric: where c denotes a set of dihedral angles, m and n indices enumerate Monte-Carlo ensemble members, and k runs over the elements of c. An UPGMA agglomerative hierarchical clustering tree 35 was built with this distance metric using Matlab Statistics and Machine Learning Toolbox. 36 The initial structure for the DOTA-M8 tag attached to a cysteine side chain by a disulphide bond was generated using ZORA DFT B3LYP/SVP energy minimisation in ORCA. 37 The bond lengths were then xed and only the dihedral angles were varied, with the energies computed using the UFF model 38 and populations estimated using the Boltzmann distribution. 39 We followed the same approach as Joseph et al., 40 with soened Lennard-Jones potentials that account for the librational motion in the protein environment 41the equilibrium inter-atomic distances were scaled down until convergence was achieved in the distance distributions. Direct inspection of the energy proles along each of the dihedral angles in the DOTA-M8 linker identies 2

Results and discussion
The rst step in the paramagnetic centre probability density reconstruction is to nd its approximate location. The process is illustrated graphically in Fig. 1. The initial guess for the probability density is a uniform distribution over all points in space that are realistically accessible to the tagoutside the van der Waals radius of the protein and no further than 12Å from its surface. The initial guess for the magnetic susceptibility tensor is obtained from the point model t. The initial localisation of the paramagnetic centre distribution is typically achieved in 50 Newton-Raphson iterations or fewerabout a minute on the wall clock when a Tesla K40 coprocessor card is used to run the FFTs on a 128 Â 128 Â 128 point grid.
Once the approximate location of the paramagnetic centre becomes clear, the renement of its distribution on a ner grid can proceed with a much reduced variational volume that only involves the region of space immediately adjacent to the approximate location of the tag. An example of such a volume is given in Fig. 2 (red cube in the right panel); a 20 Â 20 Â 20Å cube is in practice sufficient.
The optimisation is then performed repeatedly for different values of the regularisation parameter l in eqn (4). The resulting L-curve is shown as an inset in the le panel of Fig. 2; the corresponding curvature plot is in the middle panel. The optimum value of the regularisation parameter (indicated with a red circle) is calculated and the optimisation is performed again with that value. This yields the paramagnetic centre probability density (red cloud in the right panel) on a ne grid, as well as the plot of the back-calculated pseudocontact shis against the experimental ones (Fig. 2, le panel). A 256 Â 256 Â 256 point grid is in practice sufficient; the calculation takes a few hours on a Tesla K40 card.
Once the dra probability density is obtained, a different least squares optimisation is run, this time with respect to the ve independent elements of the effective magnetic susceptibility tensor. The tensor is updated and the probability density reconstruction procedure described above is repeated. The whole procedure is performed multiple times until self-consistency is achieved between c and r(r).
The procedure described above relies on two signicant assumptions. Firstly, the protein structure is treated as rigid and only the paramagnetic centre is assumed to be delocalised. This is an approximationin a real protein structure the Fig. 1 Evolution of the probability density of the Tm 3+ ion attached to C220 of the S220C mutant of human carbonic anhydrase II with a DOTA-M8 tag during the error functional optimisation process. The initial guess is a uniform distribution within the volume that is at least 2.0Å from all atoms of the protein itself and at most 12Å from any of its atoms, corresponding to the region of the space realistically accessible by the Tm 3+ ion in a tag attached anywhere on the protein surface. As the optimisation proceeds, the probability density gradually becomes zero in the locations that are not consistent with the experimental PCS data. At the end of the optimisation, the probability density is localised, subject to the standard accuracy conditions associated with Tikhonov regularisation, 22 in the region of space actually accessible to the Tm 3+ ion.
pseudocontact shis are also averaged over the distributions in the nuclear positions. From NMR data, in well-dened structures these have position distributions within about 0.4Å for backbone atoms and 1.0Å for all heavy atoms. 42 It is therefore to be expected that the paramagnetic centre distribution obtained from the PCS data would be broader than the real one by approximately that amount. Secondly, eqn (2)-(6) rely on the magnetic susceptibility tensor being the same at each point in the tag distribution. This is not necessarily true because the orientation of the tag can vary. This matter has recently been studied in detail by Shishmarev and Otting; 14 their conclusion was that a single effective c tensor can describe the PCS eld reasonably well, even in the presence of signicant tag mobility. A recent experimental study by Abdelkader et al. has also concluded that using an effective magnetic susceptibility tensor to mask its orientational distribution is a good approximation. 43 The algebraic structure of eqn (2) suggests that local variations in c can be compensated by local variations in the probability densitythe practical consequences of the constant effective magnetic susceptibility tensor assumption are therefore minor ripples in the probability density. A technical analysis of the accuracy of this approximation is given in the ESI; † the conclusion is that the resulting uncertainty is multiplicativeit would never generate probability density where there was none; it can only scale the true density by a factor related to the norms of the susceptibility tensors involved.
The results of the paramagnetic centre probability density reconstructions for S50C, S166C, S217C and S220C mutants of hCA-II with a Tm 3+ containing DOTA-M8 tag attached to the corresponding cysteines are presented in Fig. 3 and 4. As could be expected, the paramagnetic centre locations predicted by the point model ts (dark grey three-dimensional crosses in Fig. 3) are located close to the centroids of the probability density distributions computed by regularization (coloured translucent bubbles). The distributions also overlap signicantly with the Tm 3+ ion positions predicted by the rotamer library (swarms of coloured spheres), providing an independent experimental conrmation of the validity of the rotamer library approach. 40,41 An important secondary contribution to the chemical shi in paramagnetic systems arises from the residual anisotropic chemical shis (RACS) that are caused by the weak alignment of the magnetic susceptibility tensor by the applied magnetic eld. This effect was recently studied in detail by Otting et al., who estimated the RACS correction magnitude for backbone 15 N nuclei to be about 0.1 ppm for a dysprosium ion rigidly coordinated inside a protein structure. 44 In the context of this work, Fig. 2 Diagnostic information and the outcome of a typical paramagnetic centre probability density reconstruction run. After the initial localisation stage (Fig. 1), the region of space in which the probability density is allowed to vary is chosen (right panel, red cube). Multiple reconstruction runs with different values of the regularisation parameter are performed to obtain the L-curve (left panel, cut-in). The optimum regularisation parameter is extracted as the maximum curvature point on the L-curve (middle panel). The final reconstruction is performed to obtain the probability density (right panel, red cloud) and the fitting plot (left panel). Blue circles in the left panel correspond to the point model fit and the red dots to the probability density fit. the RACS correction was assumed to be negligiblea thulium ion (weaker PCS than dysprosium) at the end of a exible linker would generate a much smaller chemical shi correction than 0.1 ppm, which is itself smaller than the scatter observed in Fig. 2. In situations when RACS are suspected to be signicant, we recommend running the reconstruction using only proton PCS data because the effect is negligible for protons. Fig. 4 presents the comparison between the distance distributions between lanthanide ions in tags at different sites obtained using three physically different methods: DEER, 34 PCS (this work) and rotamer libraries. 40,41 For three out of four tagging site pairs, the agreement of the PCS data with the DEER data is very good and signicantly better than the agreement of the rotamer library prediction with the DEER data. In the remaining case, both the rotamer library and the PCS prediction deviate from the DEER data by the same amount. This indicates that PCS-based reconstruction of the spatial distribution of the paramagnetic centre performs better than rotamer libraries, although testing on a broader range of proteins would be necessary to make that conclusion in a denitive way. One of the possible explanations for the difference between PCS and DEER reconstructions for the 50-166 dataset is the presence of structural changes caused by the double mutationindividual S50C and S166C mutations (used for PCS) might not have inuenced the overall protein geometry in a detectable way, but in the double mutant (used for DEER) the changes could be signicant. This conjecture is supported by the fact that a related 166-220 double mutant is completely unstable and could not be expressed in a non-degraded form. The observed difference should not therefore be held against either method; it yields a useful structural insight.
The question of tag probability density reconstruction is particularly pertinent to the many ongoing efforts to characterise domain mobility in proteins. [45][46][47][48] Pseudocontact shi is a convenient parameter for those studies because the timescale of its emergence (i.e. the unpaired electron magnetisation equilibration time in a lanthanide) is in the picoseconds, and the time scale of its observation (i.e. the reciprocal frequency difference between the signals in the paramagnetic NMR spectrum) is in the milliseconds.
The former is much faster than protein domain mobility, and the latter is much slower, meaning that the probability density is well dened and simply reects structural heterogeneitythe corresponding theory is not troubled by the local dynamics effects that make nuclear spin relaxation theory so complicated.
The following would be a reasonable usage scenario for the method described above. In a multi-domain protein or nucleic acid, one of the domains should be tagged with a lanthanide. Pseudocontact shis measured in the same domain should be used to run a probability density reconstruction. The resulting cloud would be a measure of how rigidly the tag is immobilised relative to its home domain. At the second stage, pseudocontact shis measured in the other domains should be used to reconstruct the volume that is available to the tag; that volume is an indication of the volume explored by its home domain relative to other domains. The width of the tag distribution in its home domain would then be a measure of the uncertainty in the resulting conformational mobility conclusions. Rigidly immobilised and highly predictable tags 49 are therefore likely to be benecial.
On the detailed map of protein mobility analysis methods recently published by Ravera et al., 50 the PCS technique described in this paper belongs to the L-curve class, with a signicant difference that the penalty functional is not molecular energy (of which there is no notion in the probability density formalism), but the more traditional Laplacian norm. 22 There exists a possibility of introducing a contrast functional similar to the maximum entropy one, 20 but we would not recommend using it because it is hard to justify on physical grounds, and also because the distance distribution widths are already in good agreement with other methods (Fig. 4).

Conclusions and future work
A probability density of the paramagnetic centre may be extracted from PCS data using the recently discovered partial differential equation for PCS 20,21 and the Tikhonov regularisation method. The resulting technique is experimental in the same sense as DEER 33it requires regularisation at the data processing stagebut it explores a very different range of conditions: a room-temperature solution rather than a glass at cryogenic temperatures. This makes it highly complementary to DEER because the difference between static structural heterogeneity at cryogenic temperatures and the dynamic structural ensemble at ambient conditions becomes easy to observe. The PCS method also has the advantage of only requiring NMR equipment, which is easier to come by than cryogenic pulsed EPR gear. In many cases it would not even need a chemically attached tag, since about 30% of all proteins coordinate metal ions naturally and can usually accommodate a variety of PCSfriendly ions. 8,51,52 Some mobility is suspected to exist at those coordination sites; 21 it would be an interesting target for further exploration using 3D reconstructions.
Because the extracted distributions have the physical meaning of probability densities, multiple independent datasets (for example, from different metals or different structures in a bundle) may be combined by multiplication. We did not explore this matter further, but it bears notice that the possibility exists.
The probability density reconstruction technique described above is also important because it provides an independent experimental validation for the DEER methodso far, the distributions of the tag at each labelling site could only be modelled, and no experimental technique was available to check the results, except for DEER itself. The good agreement on both the centres and the widths of the distance distributions shown in Fig. 4 is a strong endorsement of the two-electron dipolar spectroscopy results.