Pseudocontact shift analysis
This tutorial shows how to use Spinach for analysis of paramagnetic NMR data. The examples are designed in a way to familiarize users with Spinach functions and basic theory of paramagnetic shift.
Paramagnetic chemical shift consists of two contributions: Fermi-contact (FC) and pseudocontact shift (PCS). It can be expressed via trace of the product of hyperfine (A) and susceptibility (χ) tensors. There are four different approaches for PCS analysis currently presented in Spinach. Difference in the presented methods appears only at the evaluation of the hyperfine tensor while susceptibility is either defined by user or optimized. Note that PCS depends only on the traceless part of χ, whereas Fermi contact part of the paramagnetic shift depends on the Tr(χ).
Contents
Point approximation
The simplest approach is based on the point-dipole approximation where dipolar hyperfine tensor depends on the position of a nucleus relative to the paramagnetic center \(\left(\vec r\right)\) as \(\frac{1}{{4\pi {r^3}}}\left( {\frac{{\vec r \cdot {{\vec r}^{\rm{T}}}}}{{{r^2}}} - \frac{1}{3}} \right)\). Therefore PCS is defined as \(\sigma = \frac{1}{{4\pi {r^3}}}{\rm{Tr}}\left( {\left( {\frac{{\vec r \cdot {\vec r^{\rm{T}}}}}{{{r^2}}} - \frac{1}{3}} \right) \cdot {\bf{\chi }}} \right)\). This approach is valid at the distance where paramagnetic center can be viewed as a point. To solve a direct problem use the function ppcs.m. This function takes structure of a molecule and susceptibility tensor as an input and gives PCS at the nuclei as an output. To solve an inverse problem use ippcs.m. The function evaluates susceptibility tensor and position of the paramagnetic center for provided PCS data with corresponding nuclear coordinates.
Example 1a. Direct problem for PCS using the point-dipole approximation: Metal porphyrin proton PCS.
In this example we will compute PCS on protons of the Co(II)/Cu(II) porphyrin and see how magnetic anisotropy affects PCS. The structure of the molecule in proper coordinate system is given at M_porphyrin.txt.
- Enter the coordinates of protons:
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];
- Enter the susceptibility tensor in cubic Angstrom in SI. To compute susceptibility tensor from g-tensor for S=1/2 at room temperature use the formula \({\chi _i} = {\mu _0}\mu _B^2\frac[[:Template:G i^2]]{{4{k_B}T}} \approx 1.9571\frac[[:Template:G i^2]]{{T[K]}}[{{Å}^3}]\). Note that the coordinate system for tensor is defined by molecular frame.
% g-tensor
g_Co=[3.0 3.0 2.0];
g_Cu=[2.0 2.0 2.2];
% susceptibility tensor in A^3 in SI
T=298; %K
chi_Co=1.9571/T*diag(g_Co.^2);
chi_Cu=1.9571/T*diag(g_Cu.^2);
- Define position of paramagnetic center
mxyz=[0 0 0];
- Compute proton PCS with different susceptibility tensors
pcs_Co=ppcs(nxyz,mxyz,chi_Co)
pcs_Cu=ppcs(nxyz,mxyz,chi_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. The shifts (in ppm) are
Co Cu
5.95 -1.00
9.35 -1.57
In the case of easy-plane anisotropy like in Co the PCS in equatorial plane is positive, and in the case of easy-axis anisotropy - negative. In both cases gx and gy are the same so the PCS field has a shape of second rank spherical harmonic Y20 (shown on the left). In the case of non-zero rhombisity the shape of PCS field is a linear combinations of several second rank spherical harmonics (the case of limiting rhombisity is shown on the right).
Example 1b. Inverse problem for PCS using point-dipole approximation. Calbidin D9K PCS.
In this example we will process the experimental PCS data for calbindin D9K in which one of the two calcium ions has been replaced by a thulium ion from [Balayssac, S.; Jiménez, B.; Piccioli, M. Journal of Biomolecular NMR 2006, 34, 63]. The matrix of coordinates from PDB file 1IGV atoms with corresponding PCS from the paper is stored in calbidin_xyz_pcs.txt
- Import the data from calbidin_xyz_pcs.txt to x y z and expt_pcs
- Solve the inverse problem with initial guess of the Tm position [-5 5 -15]. One may need several tries with different initial guesses to converge the solution.
mguess=[-5 5 -15];
[mxyz,chi,pred_pcs]=ippcs([x y z],mguess,expt_pcs);
- Plot experimental vs predicted PCS
plot(expt_pcs,pred_pcs,'bo');
xlabel('Experimental PCS, ppm');
ylabel('Predicted PCS, ppm');
axis equal; axis square;
- Print the position of Tm and traceless part of suceptibility tensor in molecular frame
disp('Susceptibility tensor:'); disp(chi);
disp('Point electron location:'); disp(mxyz);
You can see that the fit with the point model is quite good but there are several outliers, which we will talk about later.
Corrections to PCS from non-point source
The next level of theory that corrects the point approximation is based on the multipole expansion of a spin distribution. It works for nuclei outside the spin density bounding sphere and accounts for its anisotropy that affect PCS (see http://arxiv.org/abs/1607.08869 for details). It should be noted that any isotropic distribution of the spin gives the same PCS as the point source. Direct problem is solved using lpcs.m and inverse problem is solved using ilpcs.m. This method requires from user list of the ranks of spherical harmonics (L). If L=0 the density is isotropic and result is identical to point approximation. Normally L=[0 1 2] accounts for the major effects of the anisotropy in the spin density distribution.
Example 2a. Correction to point PCS from spin density delocalization. Copper porphyrin
Spin density in copper porphyrin is clearly not a point object. In fact the delocalized spin density over nitrogen leads to the detectable superhyperfine signal in EPR. Let's assume that 10% of spin density is localized on each nitrogen atom. How that affects proton PCS? A simple way to to account for spin delocalization is to compute multipole moments of the distribution using the function points2mult.m.
- To compute multipole moments we have to specify nuclear coordinates with corresponding spin population
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];
% spin population
rho=[0.6 0.1 0.1 0.1 0.1]';
- Compute multipole moments with ranks from 0 to 14
L=[0:14];
% multipole moments
Ilm=points2mult(xyz,mxyz,rho,L,'points')
- Compute PCS with account of spin delocalization
theor_pcs=lpcs(nxyz,mxyz,L,Ilm,chi_Cu)
The resulting PCS shows that effect of the spin delocalization on PCS is significant. Point model predicts shifts at -1.00 and -1.57 ppm and non-point model - at -1.21 and -1.81 ppm.
Example 2b. Account for spin tag mobility. Calbindin PCS
To account for possible mobility of thulium ion in the coordination cites we solve the inverse problem with additional non-point contributions using ilpcs.m. Using the experimental data from example 1b we solve the inverse problem with L=[0 1 2]:
[mxyz,Ilm,pred_pcs]=ilpcs(nxyz,expt_pcs,L,mguess)
Fit of the experimental values improves significantly with account for spin delocalization. Standard deviation drops from 1.26 to 0.67 ppm.
Hyperfine tensor from DFT
Another way to account for spin delocalization is just taking the hyperfine tensors from DFT. Spinach has two parser scripts: gparse.m - for outputs from Gaussian and oparse.m - for ORCA outputs, which help to extract the computed hyperfine tensors specified in molecular frame. It should be noted that only anisotropic part of the tensor affects PCS. For light nuclei it is just dipolar hyperfine, for heavy elements the orbital part also contributes to anisotropic part of hyperfine tensor. The function hfc2pcs.m computes PCS by multiplying hyperfine tensor with susceptibility tensor and taking its isotropic part.
Example 3a. DFT hyperfine. Copper porphyrin
To get the DFT hyperfine for protons one can run a simple ORCA job
! B3LYP def2-SVP
*xyzfile 0 2 Cu_porph.xyz
%eprnmr
Nuclei = all H {adip}
end
An example output file is located in the folder examples/nmr_paramag in Spinach. To get the PCS using these hyperfine tensor use the function hfc2pcs.m
% Parse the ORCA log
props=oparse('cu_porph_hfc.out');
% Extract hyperfine tensors
hfcs=props.hfc.full.matrix(26:37);
% Compute DFT PCS
for n=1:numel(hfcs)
dft_pcs_cu(n)=hfc2pcs(gauss2mhz(hfcs{n}),chi_cu,'1H'); %#ok<AGROW>
end
PCS from the simplistic model with spin population on nitrogen Example.2a accounts for delocalization only partially. DFT hyperfines give even larger absolute shifts: -1.53 and -1.90 ppm (it was -1.21 and -1.81 ppm in Example.2a)
Distributed model
The most general approach that works for any density distribution is 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 density distribution defined on three dimensional grid.
"Example 4a." Numerical solution for direct problem. Copper(II) porphyrin
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
% 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 using a small subsystem (i.e. protein active cite) and recompute PCS in full system (entire protein) by means of kpcs.m.
"Example 4b." Numerical solution to inverse problem. Calbindin.
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.
