Pseudocontact shift analysis
This is a tutorial page that goes through the different types of pseudocontact shift analysis that are available in Spinach. It also summarises basic facts from the theory of pseudocontact shift. Please refer to our recent publications on the subject for the technical details of the three-dimensional PCS analysis process, and to standard reviews (from Florence and Otting groups) for the point PCS model description.
Contents
Introduction
The presence of rapidly relaxing unpaired electrons generates extra nuclear chemical shifts (called paramagnetic shifts) because the unpaired electrons respond nearly instantaneoulsy to the changes in the applied magnetic field and therefore provide additional nuclear shielding that is much stronger then what would be available from a closed-shell electronic structure. Paramagnetic chemical shift has two contributions:
- Contact shift, which arises from the isotropic part of the hyperfine coupling between the electrons and the nuclei located in their close vicinity (typically less than 5 Angstrom).
- Pseudocontact shift, which is generated by the long-range (up to 50 Angstroms) electron-nuclear dipolar interaction.
Both contributions may be obtained from classical physics expressions if the electron-nuclear hyperfine coupling tensor and the electron magnetic susceptibiity tensors are known. The contact shift depends on the isotropic part of the magnetic susceptibility tensor, and the pseudocontact shift on the anisotropic part.
Point dipole approximation
This is historically the first and the simplest approach. It is based on the point dipole approximation for the electron. The expression for the pseudocontact shift for a nucleus at position [math]\mathbf{r}[/math] relative to the electron with a magnetic susceptibility tensor [math]\mathbf{\chi }[/math] is:
This approach is valid for the distances at which the paramagnetic centre may be viewed as a point - typically at least 10 Angstroms away. The PCS of any nuclei located at closer distances would not fit this model. The Spinach function that calculates pseudocontact shifts from user-specified coordinates and electron susceptibility tensors is ppcs.m. The function that extracts the paramagnetic centre position and the susceptibility tensor from user-specified nuclear coordinates and pseudocontact shifts using the point electron model is ippcs.m.
Note that the magnetic susceptibility tensor may be estimated from the g-tensor using the following expression:
The function that performs this estimation for an arbitary g-tensor is g2chi.m. The resulting estimate is only valid at high temperature.
Example 1a. PCS calculation using the point dipole approximation
In this example we compute the PCS on the protons of a Co(II)/Cu(II) porphyrin and see how magnetic anisotropy affects PCS. The following steps should be taken:
- Specify the coordinates of all protons:
% Porphyrin ring proton coordinates
nxyz=[ 4.551635888 2.658552774 0.000000000
2.658552774 4.551635889 0.000000000
-2.658552774 4.551635888 0.000000000
-4.551635889 2.658552774 0.000000000
-4.551635888 -2.658552774 0.000000000
-2.658552774 -4.551635889 0.000000000
2.658552774 -4.551635888 0.000000000
4.551635889 -2.658552774 0.000000000
4.533874147 0.000000000 0.000000000
0.000000000 -4.533874147 0.000000000
-4.533874147 0.000000000 0.000000000
0.000000000 4.533874147 0.000000000];
- Specify (or compute from the g-tensor) the susceptibility tensor (cubic Angstrom) in the same coordinate system as the one used for the proton coordinates above:
% Co(II) g-tensor
g_co=diag([3.0 3.0 2.0]);
% Cu(II) g-tensor
g_cu=diag([2.0 2.0 2.2]);
% Curie susceptibility tensors
chi_co=g2chi(g_co,298);
chi_cu=g2chi(g_cu,298);
- Specify the metal position
% Metal position
mxyz=[0 0 0];
- Call ppcs.m function
% PCS calculation
point_pcs_co=ppcs(nxyz,mxyz,chi_co);
point_pcs_cu=ppcs(nxyz,mxyz,chi_cu);
- Print the output to the console
% Output
disp('Pseudocontact shifts [Co, Cu], ppm');
disp([point_pcs_co point_pcs_cu]);
In the output you see that there are two groups of protons: four with larger absolute shift and eight with smaller shift. These correspond to two symmetry unique protons in the D4h symmetric molecule:
PCS [Co, Cu], ppm
5.9471 -0.9991
5.9471 -0.9991
5.9471 -0.9991
5.9471 -0.9991
5.9471 -0.9991
5.9471 -0.9991
5.9471 -0.9991
5.9471 -0.9991
9.3458 -1.5701
9.3458 -1.5701
9.3458 -1.5701
9.3458 -1.5701
Example 1b. PCS fitting using the point dipole approximation
In this example we will process the experimental PCS data from [Balayssac, S.; Jiménez, B.; Piccioli, M. Journal of Biomolecular NMR 2006, 34, 63] for calbindin D9K, in which one of the two calcium ions has been replaced by a thulium ion. The corresponding example file is located in nmr_paramag/calbindin directory of the example set supplied with Spinach.
The calculation proceeds in the following stages:
- Import the data. It may be pyped in explicitly, as above. In this case, atomic coordinates and pseudocontact shifts are supplied in a Matlab file:
% Load experimental data
load('tm_1igv_pcs.mat');
- Call the point model fit function ippcs.m with the initial guess of the thulium metal position set to [-5 5 -15]:
% Solve the inverse problem
[mxyz,chi,pred_pcs]=ippcs([x y z],[-5 5 -15],expt_pcs);
The other variables are contained in the file loaded in the previous line. In practice, several attempts with different paramagnetic centre location guesses may be necessary.
- Plot the results
% Plot experimental vs predicted PCS
figure(); plot(expt_pcs,pred_pcs,'bo'); hold on; grid on;
plot([min(expt_pcs) max(expt_pcs)],[min(expt_pcs) max(expt_pcs)],'r-');
xlabel('Experimental PCS, ppm'); ylabel('Predicted PCS, ppm');
- Print the position of Tm and the traceless part of the magnetic susceptibility tensor in molecular frame
% Report the parameters
disp('Susceptibility tensor:'); disp(chi);
disp('Point electron location:'); disp(mxyz);
The output appears below:
Susceptibility tensor:
0.0585 -0.0878 0.0729
-0.0878 -0.0773 0.0341
0.0729 0.0341 0.0189
Point electron location:
3.6610 17.7074 15.4140
It is clear from the figure that the point dipole model does fit the data, but does not fit it particularly well.
Magnetic multipole approximation
The next level of theory that improves on the point dipole approximation is based on the spherical harmonic expansion of the electron probability density:
It works for nuclei positioned outside the bounding sphere of the electron probability density. PCS calculation is performed using lpcs.m function and the fitting is performed using ilpcs.m function. In both cases, the list of [math]l[/math] ranks to be included into the expansion must be supplied by the user.
Example 2a. Computing PCS using the multipole model
The spin density in the copper porphyrin from Example 1a is clearly not a point object; some density is located on the nitrogen atoms because the corresponding isotropic hyperfine couplings are clearly visible in ESR spectra. If the spin populations on all relevant atoms are known (for example, from electronic structure theory calculations, or as fitting parameters), the multipole moments of the overall distribution may be estimated using points2mult.m function. If the spin density is given on a 3D grid, the same function would calculate its multipole moments by numerical integration. A typical calculation would proceed in the following stages:
- Specify nuclear coordinates and the corresponding Mulliken spin populations:
% Nuclear coordinates
xyz=[ 0.000000000000 0.000000000000 0.000000000000
-1.452280156504 1.452280156496 0.000000000000
1.452280156504 -1.452280156496 0.000000000000
-1.452280156496 -1.452280156504 0.000000000000
1.452280156496 1.452280156504 0.000000000000];
% Multipole centre coordinates
mxyz=[0 0 0];
% Mulliken spin populations
rho=[0.6 0.1 0.1 0.1 0.1]';
- Compute multipole moments relative to the centre, in this case up to spherical harmonic rank 14
% Compute multipole moments
Ilm=points2mult(xyz,mxyz,rho,0:14,'points');
- Get the susceptibility tensor from the g-tensor
% Cu(II) g-tensor
g_cu=diag([2.0 2.0 2.2]);
% Curie susceptibility tensor
chi_cu=g2chi(g_cu,298);
- Compute PCS using the multipole moments
% Get the pseudocontact shifts
theor_pcs=lpcs(nxyz,mxyz,0:14,Ilm,chi_cu)
A visual inspection of the resulting PCS values shows that the effect of the spin delocalisation on the PCS is significant: the point dipole model predicts shifts at -1.00 and -1.57 ppm, and the multipole model at -1.21 and -1.81 ppm, respectively. This example is located in nmr_paramag directory of the example set supplied with Spinach.
Example 2b. Fitting PCS to the multipole model
Fitting experimental data to the multipole model can account for the possible mobility of the lanthanide centre; this is often the case in proteins. If we repeat Example 1b identically, except for the solver line:
% Solve the inverse problem
[mxyz,chi,~,pred_pcs]=ilpcs([x y z],expt_pcs,[0 1 2],[-5 5 -15]);
where now the fit is to the multipole model with up to the second spherical rank in the probability density distribution, we would see a much better fit - the residual drops by almost a factor of three. Specifying [math]l=\left[ \begin{matrix} 0 & 1 & 2 \end{matrix} \right][/math] is typically sufficient; [math]l=0[/math] makes the calculation equivalent to the point dipole approximation. Spherical ranks higher than 4 are not well constrained by fitting.
Using hyperfine tensors from electronic structure theory
Another way to account for the paramagnetic centre delocalisation is to use hyperfine tensors obtained from electronic structure theory calculations. Those are integrals over the unpaired electron probability density and therefore include all delocalisation effects. Only anisotropic part of the hyperfine tensor affects the PCS; for light nuclei it is just the dipolar hyperfine interaction, for heavy elements the orbital part also contributes.
Spinach has two electronic structure theory package parsers: gparse.m for outputs from Gaussian and oparse.m for outputs from ORCA. Both extract hyperfine tensors from the logs produced by the program. The function hfc2pcs.m then computes PCS by multiplying the traceless part of the hyperfine tensor by the traceless part of the susceptibility tensor and taking the trace of the result. To keep the traces (and therefore the contact part of the full paramagnetic shift), use hfc2pms.m function. The general equation is:
See http://dx.doi.org/10.1039/C4CP03106G for further information. Note that this method essentially requires the molecule to be rigid; realistically, it is only applicable to small molecules, for which it is strongly recommended.
Example 3a. Computing PCS using hyperfine tensors from DFT
You should refer to ORCA and Gaussian manuals to find out how to calculate hyperfine tensors. An example ORCA output file is located in the nmr_paramag folder of the example set supplied with Spinach. The calculation proceeds in the following stages:
- ORCA log file is parsed by Spinach:
% Parse the ORCA log
props=oparse('cu_porph_hfc.out');
- Relevant hyperfine tensors are extracted:
% Extract hyperfine tensors hfcs=props.hfc.full.matrix(26:37);
- PCS values are computed
% Compute DFT PCS
for n=1:numel(hfcs)
hfc_pcs_cu(n)=hfc2pcs(gauss2mhz(hfcs{n}),chi_cu,'1H');
end
The same susceptibility tensor as in Examples 1a and 2a is used here. The point model in Example 1a was predicting PCS at -1.00 and -1.57 ppm; the multipole model in Example 2a put the PCS at -1.21 and -1.81 ppm; DFT hyperfines give even larger PCS values: -1.53 and -1.90 ppm.
Distributed paramagnetic centre model
The most general approach that works for any paramagnetic centre probability distribution is based on the partial differential equation obeyed by the pseudocontact shift and described in (http://dx.doi.org/10.1039/C4CP03106G) and implemented in the functions kpcs.m for the direct problems and ipcs.m for the inverse problems. This approach is fully numeric and operates with the probability density defined on three-dimensional grids.
Example 4a. Computing PCS using DFT spin density
In this example the protons PCS of Copper(II) porphyrin are computed from the spin density generated by DFT using kpcs.m. ORCA_PLOT offers several options to generate the spin density file. In this example we use the "3D simple format" file "cu_porph_sd80.spindens.3d". To import the density into the format that kpcs.m needs we do:
% Import spin density from ORCA in "3D simple fromat"
A=importdata('cu_porph_sd80.spindens.3d');
% read number of points along [x y z]
dim=str2num(A.textdata{2});
% read coordinates of the corner
min=str2num(A.textdata{3});
% read step size [dx dy dz]
d3r=str2num(A.textdata{4});
% form the absolute density cube
density=reshape(abs(A.data),dim(2),dim(3),dim(1));
density=permute(density,[2 3 1]);
% normalize the density
density=density/(simps(simps(simps(density)))*prod(d3r));
% compute extents
ext=[min(1) (min(1)+(dim(1)-1)*d3r(1))...
min(2) (min(2)+(dim(2)-1)*d3r(2))...
min(3) (min(3)+(dim(3)-1)*d3r(3))];
The numerical costs for kpcs.m are defined by the grid size. New version of kpcs.m uses the fast Fourier transform (FFT) method that improves the scaling significantly (for details see http://arxiv.org/abs/1607.08869). FFT has the periodic boundary conditions (PBC), therefore we have to make sure that borders are far enough from the molecule by padding the density array with zeros:
% Pad the density with zeros to avoid PBC effects pad_size=2; pad_density=padarray(density,pad_size*size(density),0,'both'); pad_ext=(2*pad_size+1)*ext;
The actual calculation takes place here:
% Solve PDE using FFT [pde_pcs_cu,pcs_cube]=kpcs(pad_density,chi_cu,pad_ext,nxyz,'fft');
Comparison of the PCS from the DFT hyperfine and from the partial differential equation (PDE) with imported DFT spin density defined on 80x80x80 grid
% Compare HFC PCS with PDE PCS
disp('Pseudocontact shifts [hfc, pde], ppm');
disp([hfc_pcs_cu' pde_pcs_cu]);
shows the good agreement (-1.53 -1.90 ppm from DFT and -1.43 -1.84 ppm from PDE). Using this method one can compute DFT spin density for a small subsystem (i.e. protein active cite) and recompute PCS in full system (entire protein) using kpcs.m without a need to do DFT for the full system.
Example 4b. Reconstructing lanthanide tag distributions in proteins
This example shows how to extract the density distribution of paramagnetic center from PCS data. The input file for this example is located in Spinach\examples\nmr_paramag\calbindin/dft_density.m together with all supporting files.
