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.prox_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 prox_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.

Strychnine HSQC simulation and isotopomer loop

HSQC simulation for strychnine.

Each isotopomer is to be simulated independently in a parallel loop:

    parfor n=1:numel(subsystems)

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

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

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

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

        % Form States signal
        fid=f1_pos+conj(f1_neg);

        % F1 Fourier transform
        spectrum=spectrum+fftshift(fft(fid,parameters.zerofill(1),2),2);

    end

Here, the basis set is built inside the loop; 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. 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. Both States components are multiplied by square-cosine windows in both dimensions. The direct dimension is transformed first, the States signal is then formed and the indirect dimension is transformed and added to the accumulated spectrum. Because each subsystem is independent, the loop is a natural use case for Matlab parfor.

Strychnine HSQC plotting

The final stage is to open a figure, scale it, and plot 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.

Camphor spin system import

Chemical structure of camphor.

Create a new function file. The difference with the previous example is that here we will import spin system information 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);

CLIP-HSQC simulation for camphor.

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. We then specify the magnetic field:

    sys.magnet=14.1;

and replace 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 parts of each tensor and replaces only the isotropic component with the supplied value.

Camphor CLIP-HSQC simulation

Use exactly the same basis set and other options as we have used above for strychnine, and the same parallel loop over isotopomers. The only difference now should be that CLIP-HSQC pulse sequence is called instead of HSQC:

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

and the plot must use 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');

Exercises

These are open-ended exploration tasks, see the built-in pulse sequence list for what is available. Suggestions:

1. Run other pulse sequences in the same family, for example, hmqc.m
2. Set the working J-coupling to a deliberately incorrect value and observe the effect.
3. Increase and decrease the number of acquired points in parameters.npoints parameter.
4. Turn the decoupling off in one or both dimensions and observe the changes in the simulated spectrum.
5. Decrease the J-coupling drop threshold in the spin system import command and observe the effect on the simulation time.

Version 2.12, authors: Ilya Kuprov