Heteronuclear NMR simulations

Heteronuclear liquid-state NMR simulations in Spinach are set up in the same way as homonuclear simulations, but the pulse sequence parameters must specify multiple spin channels. This tutorial follows two examples from the Spinach example set: hsqc_strychnine.m, which simulates a phase-sensitive natural-abundance ¹³C HSQC spectrum of strychnine, and clip_hsqc_camphor.m, which simulates a natural abundance CLIP-HSQC spectrum of camphor.

Strychnine spin system

Chemical structure of strychnine

Spinach examples are functions rather than scripts: each run starts from a clean Matlab workspace. We must therefore open a new file and begin with with a function declaration:

    function hsqc_strychnine()

Unlike the previous tutorial where we had specified the spin system manually, here we use one of the standard spin systems supplied with Spinach - the strychnine spin system:

    [sys,inter]=strychnine({'1H','13C'});

The strychnine helper function returns the sys and inter input structures for the selected isotopes. In this example, proton and carbon-13 spins are imported. The strychnine data file contains proton shifts and couplings, carbon and nitrogen shifts and coordinates, and one-bond ¹³C-¹H J-couplings. The magnetic field is then set in Tesla:

    sys.magnet=5.9;

The next line requests the "greedy" parallelisation option:

    sys.enable={'greedy'};

When Spinach starts a new parallel pool, this option means that every worker process is allowed to use every CPU core.

Basis set specification

The basis set is the same connectivity-adaptive Liouville-space basis used by many liquid-state examples:

    bas.formalism='sphten-liouv';
    bas.approximation='IK-2';
    bas.connectivity='scalar_couplings';
    bas.space_level=1;

The IK-2 approximation builds a restricted state space rather than the complete Liouville space, and scalar_couplings instructs the basis generator to use the scalar-coupling graph for connectivity. The space_level setting controls how far spatially local correlations are allowed into the basis: since HSQC does not rely on spatial proximity, this parameter is set to 1.

HSQC sequence parameters

The following parameters are required by the HSQC pulse sequence function:

    parameters.J=140;
    parameters.sweep=[10000 3000];
    parameters.offset=[4000 1000];
    parameters.npoints=[128 128];
    parameters.zerofill=[512 512];
    parameters.spins={'13C','1H'};
    parameters.decouple_f1={'1H'};
    parameters.decouple_f2={'13C'};
    parameters.axis_units='ppm';

Here parameters.J is the working heteronuclear scalar coupling in Hz, used to set the transfer delay. The two entries in parameters.sweep are the F1 and F2 sweep widths in Hz, and the two entries in parameters.offset are the corresponding transmitter or receiver offsets in Hz. The point counts in parameters.npoints are the acquired F1 and F2 dimensions, and parameters.zerofill gives the Fourier transform sizes used during processing. The spin list {'13C','1H'} means that the indirect F1 dimension is carbon-13 and the directly detected F2 dimension is proton. The decouple_f1 field lists spins that receive midpoint 180-degree refocusing pulses during F1 evolution, and decouple_f2 lists spins decoupled during F2 acquisition.

Spin system construction

We now create the spin system object from the information we have supplied above:

    spin_system=create(sys,inter);

The create function performs Spinach input processing and prints a report that should be checked for warnings. Natural-abundance carbon is handled by isotope dilution:

    subsystems=dilute(spin_system,'13C');

The dilute function returns a cell array of spin systems. With the default tuple size, each returned subsystem contains a single instance of the specified dilute isotope and all other carbon-13 spins have been removed. This makes a natural-abundance ¹³C HSQC simulation practical: the final spectrum is obtained as a sum over all single-carbon isotopomers.

We are going to be summing spectra up and we must therefore preallocate the sum:

    spectrum=zeros(parameters.zerofill(2),...
                   parameters.zerofill(1),'like',1i);

The first dimension of the matrix corresponds to the direct F2 dimension, and the second dimension corresponds to the indirect F1 dimension. The use of 'like',1i requests a complex array.

HSQC isotopomer loop

Each isotopomer is processed independently:

    parfor n=1:numel(subsystems)

The basis set is built inside the loop:

        subsystem=basis(subsystems{n},bas);

This is necessary because dilute removes basis-set information when it changes the spin system. The HSQC simulation is then run in the liquid-state context with NMR assumptions:

        fid=liquid(subsystem,@hsqc,parameters,'nmr');

The liquid context builds the isotropic Liouvillian, applies offsets and any other common sequence settings, and passes the resulting Hamiltonian, relaxation superoperator, and kinetics superoperator to hsqc.m. The HSQC sequence returns two components of the States quadrature signal:

        fid.pos=apodisation(spin_system,fid.pos,{{'sqcos'},{'sqcos'}});
        fid.neg=apodisation(spin_system,fid.neg,{{'sqcos'},{'sqcos'}});

Both components are multiplied by square-cosine windows in both dimensions. The direct dimension is transformed first:

        f1_pos=fftshift(fft(fid.pos,parameters.zerofill(2),1),1);
        f1_neg=fftshift(fft(fid.neg,parameters.zerofill(2),1),1);

The States signal is then formed:

        fid=f1_pos+conj(f1_neg);

and the indirect dimension is transformed and added to the accumulated spectrum:

        spectrum=spectrum+fftshift(fft(fid,parameters.zerofill(1),2),2);

The loop is then closed:

    end

Because each subsystem is independent, the loop is a natural use case for Matlab parfor.

HSQC plotting

HSQC simulation for strychnine.

The final stage opens a figure, scales it, and plots positive contours:

    kfigure(); scale_figure([1.5 2.0]);
    plot_2d(spin_system,real(spectrum),parameters,...
            20,[0.05 1.0 0.05 1.0],2,256,6,'positive');

The plot_2d arguments are the spin system, the real part of the spectrum, the experiment parameters, the number of contours, positive and negative contour elevation ranges, the non-linear contour-spacing curvature, the colour-map size, the colour-map curvature, and the contour sign selection. In this example the last argument requests positive contours only.

CLIP-HSQC spectrum of camphor

The clip_hsqc_camphor.m example also begins with a function declaration:

    function clip_hsqc_camphor()

The spin system is imported from a Gaussian log file:

    options.min_j=3.0; options.no_xyz=0;
    [sys,inter]=g2spinach(gparse('../standard_systems/camphor.log'),...
                     {{'H','1H'},{'C','13C'}},[31.8 182.1],options);

The gparse function reads the Gaussian output file. The g2spinach function converts the parsed electronic-structure data into Spinach input structures. The particle list imports hydrogen atoms as ¹H and carbon atoms as ¹³C. The reference vector [31.8 182.1] supplies the absolute shielding references used by the conversion to chemical shifts. The options.min_j=3.0 setting discards scalar couplings smaller than 3 Hz, and options.no_xyz=0 keeps the coordinate information.

The magnetic field is then specified:

    sys.magnet=14.1;

The next block replaces the isotropic parts of the shielding tensors with experimental chemical shifts:

    inter.zeeman.matrix=shift_iso(inter.zeeman.matrix,1:26,...
                                  [ 29.70 26.80 44.20 59.70 42.80 ...
                                   218.1  49.70 17.80 17.30  7.80 ...
                                     1.35  1.67  1.97  1.31  1.96 ...
                                     0.85  0.85  0.85  0.99  0.99 ...
                                     0.99  0.90  0.90  0.90  2.33 1.76]);

The shift_iso function preserves the anisotropic rank-1 and rank-2 parts of each tensor and replaces only the isotropic component with the supplied value. In this example, the source comment states that coordinates, shielding anisotropies, and J-couplings come from DFT, while isotropic chemical shifts come from experimental data.

CLIP-HSQC basis set and options

The basis set specification is:

    bas.formalism='sphten-liouv';
    bas.approximation='IK-2';
    bas.connectivity='scalar_couplings';
    bas.space_level=1;

The algorithmic options are:

    sys.enable={'greedy'};
    sys.tols.prox_cutoff=4.0;

As in the strychnine example, this requests Spinach's greedy parallelisation option and sets the proximity cut-off to 4 Angstrom.

CLIP-HSQC sequence parameters

The CLIP-HSQC example uses the following sequence parameters:

    parameters.J=140;
    parameters.sweep=[8000 1500];
    parameters.offset=[4000 1000];
    parameters.npoints=[128 128];
    parameters.zerofill=[512 512];
    parameters.spins={'13C','1H'};
    parameters.axis_units='ppm';

The clip_hsqc.m header requires parameters.sweep, parameters.npoints, parameters.spins, and parameters.J. As above, parameters.J is the working scalar coupling in Hz, parameters.sweep gives the F1 and F2 sweep widths, parameters.offset gives the corresponding offsets, parameters.npoints gives the acquired point counts, parameters.zerofill gives the Fourier transform sizes, parameters.spins={'13C','1H'} specifies carbon-13 in F1 and proton in F2, and parameters.axis_units='ppm' requests ppm axes.

CLIP-HSQC simulation

The spin-system construction and isotopomer generation steps are the same as in the HSQC example:

    spin_system=create(sys,inter);
    subsystems=dilute(spin_system,'13C');

The answer is again preallocated as a complex matrix:

    spectrum=zeros(parameters.zerofill(2),...
                   parameters.zerofill(1),'like',1i);

The isotopomer loop starts with basis construction:

    parfor n=1:numel(subsystems)
        
        subsystem=basis(subsystems{n},bas);

The simulation call differs only in the pulse sequence function handle:

        fid=liquid(subsystem,@clip_hsqc,parameters,'nmr');

The clip_hsqc.m sequence returns fid.pos and fid.neg as the two States quadrature components. The data processing matches the HSQC example:

        fid.pos=apodisation(spin_system,fid.pos,{{'sqcos'},{'sqcos'}});
        fid.neg=apodisation(spin_system,fid.neg,{{'sqcos'},{'sqcos'}});
        
        f1_pos=fftshift(fft(fid.pos,parameters.zerofill(2),1),1);
        f1_neg=fftshift(fft(fid.neg,parameters.zerofill(2),1),1);
        
        fid=f1_pos+conj(f1_neg);
        
        spectrum=spectrum+fftshift(fft(fid,parameters.zerofill(1),2),2);
        
    end

The summation over the isotopomer loop produces the natural-abundance ¹³C CLIP-HSQC spectrum.

CLIP-HSQC plotting

CLIP-HSQC simulation for camphor.

The example plots negative contours:

    kfigure(); scale_figure([1.5 2.0]);
    plot_2d(spin_system,real(spectrum),parameters,...
            20,[0.05 0.5 0.05 0.5],2,256,6,'negative');

The contour-level vector differs from the strychnine example, and the last argument requests negative contours only.

Version 2.12, authors: Ilya Kuprov, Luke Edwards