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 and the reaction rate matrix

Spinach input must specify different chemical forms as subsystems inside one spin system: we specify everythig in one input structure and then tell Spinach which spins belong to which chemical species. In the present tutorial each species contains one proton with it own chemical shft:

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

The last line tells create that spin 1 and spin 2 belong to different chemical subsystems. The asymmetric two-site chemical exchange process is then 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 loss rates from each form. The initial concentrations are in arbitrary concentration units. Here they are already at equilibrium. If the equilibrium concentrations are unknown, equilibrate can compute them from a valid column-conservative rate matrix and an initial concentration vector.

Formalism and basis set

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

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

For this tutorial the complete basis set is small, so no restricted state space approximation is needed. The spin system is 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.

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.

Simulation, data processing, and plotting

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.

If you reproduce this Liouvillian assembly manually, remember that relaxation and kinetics superoperators are dissipative. Use L=H+1i*R+1i*K; do not use H+R+K.

NMR spectrum of a single spin-system undergoing asymmetric chemical exchange between positions with different chemical shifts.

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

    fid=apodisation(spin_system,fid,{{'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);

Forward abd backward reaction rates affect the appearance of the spectrum. Slow exchange produces separate resonances with intensities proportional to concentrations. As the exchange rates become comparable to the frequency separation, the lines broaden and move towards the exchange-averaged position.

Exercises

The same pattern extends to more chemical sites. An extra copy of the spin system is needed 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. All exchanging chemical forms must have the same number of spins and the same atom order. See the Chemical kinetics parameters section of the Spinach manual for further information. Suggestions for open-ended exploration:

1. Decrease both reaction rates by a factor of ten and compare the line widths with the original calculation.
2. Increase both reaction rates by a factor of ten and observe how the two-site spectrum approaches the fast-exchange limit.
3. Change the initial concentration vector away from equilibrium and observe the chenges inthe specttrum resulting from transient concentration dynamics.
4. Replace the asymmetric matrix with the symmetric matrix and compare the spectrum.
5. Build a three-site model by adding a third chemical subsystem and a 3 by 3 rate matrix.

Version 2.12, authors: Ilya Kuprov