Models of chemical kinetics

First order chemical reactions are common in solution NMR: conformers interconvert, ligands bind and unbind, protons exchange, and tautomers move between structures. In a spin dynamics simulation these processes are included into the equation of motion as an additional term, called kinetics superoperator, whose job is to shuffle probability density between reactants and products. This tutorial shows how to set up a 1D liquid-state NMR simulation with first-order chemical kinetics in Spinach. It follows the exchange_asymmetric.m example from the Spinach example set.

First-order kinetic model

Under first order chemical kinetics, the concentration vector obeys a linear differential equation:

\(\frac{d}{dt}\mathbf{c}(t)=\mathbf{K}\mathbf{c}(t)\)

where \(\mathbf{c}\) is the vector of concentrations of reactants and products, and \(\mathbf{K}\) is the kinetic rate matrix. Spinach uses the column convention: an off-diagonal element \(K_{ij}\) is the rate from chemical subsystem \(j\) into chemical subsystem \(i\), and every column must sum to zero so that matter is conserved. For two-site chemical exchange:

\(\mathrm{A}\underset{k_{\mathrm{BA}}}{\overset{k_{\mathrm{AB}}}{\rightleftharpoons}}\mathrm{B}\)

the rate matrix is therefore

\(\mathbf{K}=\begin{pmatrix}-k_{\mathrm{AB}}&k_{\mathrm{BA}}\\k_{\mathrm{AB}}&-k_{\mathrm{BA}}\end{pmatrix}\).

The numbers are rates in \(\mathrm{s}^{-1}\), labelled Hz in Spinach reports.

Chemical subsystems

Spinach represents different chemical forms as different chemical subsystems inside one spin system. The molecule must be entered once for every chemical form that participates in the exchange. In the present tutorial each form contains one proton, and the two forms have different chemical shifts:

    function exchange_asymmetric()

    sys.magnet=14.1;
    sys.isotopes={'1H','1H'};
    inter.zeeman.scalar={0,3};
    inter.chem.parts={1,2};

The first proton is the NMR-active spin of chemical form A, and the second proton is the NMR-active spin of chemical form B. The line

    inter.chem.parts={1,2};

tells create that spin 1 and spin 2 are not two coupled spins in the same molecule; they are two chemical subsystems connected by kinetic exchange. For larger molecules, every chemical form must have the same isotope sequence in the same order, for example {[1 2 3],[4 5 6]} for two three-spin conformers. The corresponding nuclei in the two forms are then mapped onto each other by the kinetics module.

Rate matrix and concentrations

The asymmetric two-site model is specified by the rate matrix and by the starting concentrations:

    inter.chem.rates=[-5e2  2e3; 
                           5e2 -2e3];
    inter.chem.concs=[2e3 5e2];

The off-diagonal element in row 2, column 1 is 5e2, so form A converts into form B at 500 s^-1. The off-diagonal element in row 1, column 2 is 2e3, so form B converts into form A at 2000 s^-1. The diagonal elements are the corresponding losses from each source form. The initial concentrations are in arbitrary but consistent units. Here they are already at equilibrium because

\(500[\mathrm{A}]=2000[\mathrm{B}]\)

for [A,B]=[2000,500]. If the equilibrium concentrations are unknown, equilibrate can compute them from a valid column-conservative rate matrix and an initial concentration vector.

Basis set and construction

Linear chemical exchange of this type is supported in the spherical-tensor Liouville-space formalism:

    bas.formalism='sphten-liouv';
    bas.approximation='none';

For this one-spin-per-site tutorial the complete basis set is small, so no restricted basis approximation is needed. In larger systems the same chemical-subsystem logic applies, but all exchanging forms must use compatible basis states. The spin system is then constructed in the usual way:

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

The create function checks that each spin belongs to no more than one chemical subsystem, that the rate matrix is square and column-conservative, that the concentration vector has one entry per subsystem, and that all exchange partners have the same number and order of isotopes. The basis call then builds the Liouville-space basis used by the Hamiltonian, relaxation, and kinetics generators.

Initial state and detection

A 1D NMR acquisition needs an initial state and a detection state. For an exchange simulation, the initial transverse magnetisation should normally be weighted by the chemical concentrations:

    parameters.spins={'1H'};
    parameters.rho0=state(spin_system,'L+','1H','chem');
    parameters.coil=state(spin_system,'L+','1H');

The 'chem' option in state multiplies each requested spin state by the concentration of the chemical subsystem to which that spin belongs. This is what gives the two forms their correct starting magnetisation. The receiver coil is not concentration-weighted: it represents the observable proton transverse magnetisation summed over the chemical forms.

Acquisition parameters

The rest of the sequence parameters are the usual 1D liquid-state NMR acquisition settings:

    parameters.decouple={};
    parameters.offset=900;
    parameters.sweep=5000;
    parameters.npoints=512;
    parameters.zerofill=1024;
    parameters.axis_units='ppm';
    parameters.invert_axis=1;

The offset and sweep are in Hz. In this example the two proton resonances are separated by 3 ppm at 14.1 T, and the 5000 Hz sweep width comfortably covers both sites. The zero-filled Fourier transform size is larger than the acquired point count, which gives a smoother plotted spectrum without adding new physical information.

Running the simulation

The calculation is run in the liquid-state context with the standard acquisition sequence:

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

The liquid context builds the isotropic Hamiltonian, the relaxation superoperator, and the kinetics superoperator, applies the requested offset and rotating-frame settings, and passes all of them to acquire. The acquisition sequence forms the total Liouvillian from those pieces and evolves the initial state while observing the coil state. No separate user call to kinetics is required for ordinary pulse-sequence simulations: context functions include the kinetics superoperator automatically.

Processing and plotting

The simulated free induction decay is processed in the usual way. A mild exponential window improves visual smoothness:

    fid=apodisation(spin_system,fid,Template:'exp',6);

The spectrum is obtained by Fourier transformation and plotted as a 1D NMR spectrum:

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

    kfigure(); plot_1d(spin_system,real(spectrum),parameters);

The rates affect the spectrum during the simulated acquisition. Slow exchange produces separate resonances with intensities governed by the concentrations. As the exchange rates become comparable with the frequency separation, the lines broaden and move towards the exchange-averaged position. In asymmetric exchange, the two directions have different rates and the equilibrium populations are different; both effects are contained in the same rate matrix.

Changing the model

The same pattern extends directly to more chemical sites. Add one copy of the spin system for every site, put the corresponding spin indices into inter.chem.parts, provide an \(N\times N\) column-conservative inter.chem.rates matrix, and provide an \(N\)-element inter.chem.concs vector. The main restrictions are that all exchanging chemical forms must have the same number of spins, the same isotope order, and compatible basis states. For processes where magnetisation is transferred between spins within a fixed molecular topology, see the magnetisation-flux section of Chemical kinetics parameters; that is a different model from the first-order chemical-subsystem exchange used here.

Exercises

These are open-ended exploration tasks. Suggestions:

1. Decrease both off-diagonal rates by a factor of ten, keep the diagonal elements column-conservative, and compare the line widths with the original calculation.
2. Increase both off-diagonal rates by a factor of ten and observe how the two-site spectrum approaches the fast-exchange limit.
3. Change the concentration vector away from equilibrium and use equilibrate to compute the equilibrium populations for the same rate matrix.
4. Replace the asymmetric matrix with the symmetric matrix used in exchange_symmetric.m and compare the spectrum.
5. Build a three-site model by adding a third chemical subsystem and a 3 by 3 column-conservative rate matrix, then check that the columns of the matrix sum to zero before running create.

Version 2.12, authors: Ilya Kuprov