Simple liquid state NMR simulations

Chemical structure of rotenone.

COSY spectrum of rotenone is easy enough for a tutorial and large enough (22 spins, irregular J-coupling network) to require Spinach. The experimental data used in this tutorial comes from http://dx.doi.org/10.1002/jhet.5570250160.

Spin system specification

Open a new file in Matlab editor. Spinach would refuse to run in script mode; it requires a function declaration:

    function cosy90_rotenone()

This is necessary because this guarantees that Matlab starts the simulation from a clean state where no previously assigned variables exist. The first stage is to specify the isotopes, 22 protons in this case:

    sys.isotopes={'1H','1H','1H','1H','1H','1H','1H','1H','1H',...
                  '1H','1H','1H','1H','1H','1H','1H','1H','1H',...
                  '1H','1H','1H','1H'};

See the spin system specification section for details on how to specify spin systems. The next step is to specify the magnet field:

    sys.magnet=5.9;

then chemical shifts, in ppm for all protons:

    inter.zeeman.scalar={6.72 6.40 4.13 4.56 4.89 6.46 7.79 3.79 2.91...
                         3.27 5.19 4.89 5.03 1.72 1.72 1.72 3.72 3.72...
                         3.72 3.76 3.76 3.76};

then all scalar couplings, in Hz:

    inter.coupling.scalar{3,4}=12.1; 
    inter.coupling.scalar{4,5}=3.1; 
    inter.coupling.scalar{3,5}=1.0; 
    inter.coupling.scalar{3,8}=1.0; 
    inter.coupling.scalar{1,8}=1.0;
    inter.coupling.scalar{6,7}=8.6; 
    inter.coupling.scalar{5,8}=4.1; 
    inter.coupling.scalar{7,9}=0.7; 
    inter.coupling.scalar{7,10}=0.7; 
    inter.coupling.scalar{9,10}=15.8;
    inter.coupling.scalar{10,11}=9.8; 
    inter.coupling.scalar{9,11}=8.1; 
    inter.coupling.scalar{13,14}=1.5; 
    inter.coupling.scalar{12,14}=0.9; 
    inter.coupling.scalar{22,22}=0;

where the last line is necessary to tell Matlab that the array is 22 by 22 and all other elements are empty or zero. This completes the spin system specifications for situations where accurate relaxation theory treatment is not required.

Formalism and basis set

It is a good idea here to request "greedy" parallelisation, allowing each worker in the Matlab parallel pool to use all CPU cores:

    sys.enable={'greedy'};

The next stage is basis set specification. The complete basis set for a 22-spin system is too large and we must therefore rely on the restricted state space approximation. See the basis set specification section for further details on basis set selection. Here we will be using IK-2 basis set with Liouville space formalism with connectivity inferred from the J-coupling network:

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

The three methyl groups contain magnetically equivalent protons and this symmetry may optionally be used to reduce the calculation time:

    bas.sym_group={'S3','S3','S3'};
    bas.sym_spins={[14 15 16],[17 18 19],[20 21 22]};

This completes the basis set specification.

Simulation parameters

The specific parameters required by the COSY sequence are in the corresponding section of this manual as well as in the function header of the COSY sequence file. The parameters required by COSY are:

    parameters.angle=pi/2;
    parameters.offset=1200;
    parameters.sweep=2000;
    parameters.npoints=[512 512];
    parameters.zerofill=[2048 2048];
    parameters.spins={'1H'};
    parameters.axis_units='ppm';

where parameters.angle is the second COSY pulse angle in radians; parameters.offset is the transmitter and receiver offset in Hz; parameters.sweep is the spectral width in Hz; parameters.npoints gives the acquired point counts in the two dimensions; parameters.zerofill gives the Fourier transform sizes after zero-filling; parameters.spins selects the working isotope channel; and parameters.axis_units selects the units used on the plotted axes. This completes the specification of the spin system, of the basis set and of the experiment parameters.

Running the simulation

The next stage is to give all that information to Spinach. This is accomplished by running the two housekeeping functions:

    spin_system=create(sys,inter);
    spin_system=basis(spin_system,bas);

Both print copious output to the console which you should always read carefully because it might contain warning messages. The next stage is simulation, which we will carry out in liquid state (hence the liquid context function) with the kernel assumptions set to 'nmr', indicating common high-field NMR spectroscopy:

    fid=liquid(spin_system,@cosy,parameters,'nmr');

Data processing and plotting

COSY simulation for rotenone.

The simulation returns the two-dimensional free induction decay that should undergo apodisation (cosine bell in both dimensions is a good choice here):

    fid=apodisation(spin_system,fid,{{'cos'},{'cos'}});

and Fourier transform:

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

Finally, the real part of the spectrum is plotted with positive and negative contours:

    kfigure(); scale_figure([1.5 2.0]);
    plot_2d(spin_system,real(spectrum),parameters,...
            20,[0.01 0.1 0.01 0.1],2,256,6,'both');

The whole simulation should take a few minutes on any modern PC. Note that Matlab auto-starts the parallelization engine when it runs for the first time, that stage only happens once per Matlab session.

Exercises

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

1. Replace the COSY pulse sequence call with the E.COSY one (use ecosy.m) and change the processing part to phase-sensitive (see other E.COSY examples in Spinach). Run the simulation to obtain a phase-sensitive E.COSY spectrum of rotenone.
2. Add rapid T1/T2 relaxation (instructions are here) and see how the resolution degrades when T2 becomes shorter. Note that Spinach requires you to specify relaxation rates rather than times.
3. Switch the basis set to IK-1, set the spatial proximity level (bas.space_level) to 1 and increase the interaction level (bas.level) starting from 1 to explore the trade-off between simulation speed and accuracy in Spinach.

Version 2.12, authors: Ilya Kuprov