Chemical kinetics parameters
Contents
General notes
Spinach has a very general chemical kinetics module that can handle arbitrary reaction networks, the only restriction being that the corresponding differential equations must be linear and must have the following general form:
\(\frac{d}{dt}\left( \begin{matrix} [\text{A}] \\ [\text{B}] \\ [\text{C}] \\ \vdots \\ \end{matrix} \right)=\mathbf{K}\left( \begin{matrix} [\text{A}] \\ [\text{B}] \\ [\text{C}] \\ \vdots \\ \end{matrix} \right)\)
where \(\mathbf{K}\) is the reaction rate matrix. For example:
\(\begin{matrix} \text{A }\xrightarrow{{{k}_{1+}}}\text{ B }\underset{{{k}_{2-}}}{\overset{{{k}_{2+}}}{\longleftrightarrow}}\text{ C} \\ \frac{d}{dt}\left( \begin{matrix} \text{ }\!\![\!\!\text{ A }\!\!]\!\!\text{ } \\ \text{ }\!\![\!\!\text{ B }\!\!]\!\!\text{ } \\ \text{ }\!\![\!\!\text{ C }\!\!]\!\!\text{ } \\ \end{matrix} \right)=\left( \begin{matrix} -{{k}_{1+}} & 0 & 0 \\ {{k}_{1+}} & -{{k}_{2+}} & {{k}_{2-}} \\ 0 & {{k}_{2+}} & -{{k}_{2-}} \\ \end{matrix} \right)\left( \begin{matrix} \text{ }\!\![\!\!\text{ A }\!\!]\!\!\text{ } \\ \text{ }\!\![\!\!\text{ B }\!\!]\!\!\text{ } \\ \text{ }\!\![\!\!\text{ C }\!\!]\!\!\text{ } \\ \end{matrix} \right) \\ \end{matrix}\)
Spinach expects you to supply this matrix and the initial concentrations of the molecules, as well as to say which spins in the spin system specification belong to which molecule. Spinach can also assist in determining the steady-state concentrations (see equilibrate.m) if you do not know them.
Networks of chemical reactions
Your entire spin system must be specified multiple times, and then Spinach must be told which spins in the input stream belong to which reaction endpoint. For example, for a three-spin system where conformational exchange is happening between two configurations:
% Isotopes
sys.isotopes={'1H','13C','15N','1H',13C','15N'};
% Chemical shifts
inter.zeeman.scalar={1.0, 20.0, 15.0, 1.5, 25.0, 16.0};
% Spins 1,2,3 are in the reactant; spins 4,5,6 are in the product
inter.chem.parts={[1 2 3],[4 5 6]};
% Kinetic rate matrix (Hz)
inter.chem.rates=[-0.1 0.2; 0.1 -0.2];
% Initial concentrations (arbitrary units)
inter.chem.concs=[1.0, 2.0]
In the general case, the the parameters, supplied to create.m at the spin system specification stage, must be:
| Reaction kinetics parameters | ||
|---|---|---|
| Variable name | Variable type | Content |
| inter.chem.parts | a cell array of vectors containing integers | individual vectors in the cell array must contain the numbers of spins that belong to each of the molecules in the chemical reaction, for example {[1 2],[3 4]} indicates that spins 1 and 2 belong to the first molecule and spins 3 and 4 belong to the second molecule. |
| inter.chem.rates | a matrix of real numbers | chemical reaction rate matrix between the molecules identified in inter.chem.parts |
| inter.chem.concs | a vector of non-negative real numbers | initial concentrations of the molecules identified in inter.chem.parts |
Note that all molecules must have the same number of spins, must have those spins specified in the same order, and must have the same basis set. Only the Liouville space formalism using the spherical tensor basis set is supported for this option.
Spin state fluxes
There are rare cases when the magnetization is transferred through some indirect mechanism within a topologically fixed molecule, for example, a proton swap between two N-H groups of the same molecule. Because protons are fundamentally identical, this essentially amounts to swapping the spin state (together with all associated correlations to all external spins) between the two protons. In this case, the user should supply the matrix of magnetization flux rates between spins, for example:
% Isotopes
sys.isotopes={'15N','1H','1H'};
% Chemical shift tensors
inter.zeeman.eigs{1}=[-50 -50 100];
inter.zeeman.euler{1}=[0.1 0.2 0.3];
inter.zeeman.eigs{2}=[-5 -5 10]+1;
inter.zeeman.euler{2}=[0.2 0.3 0.1];
inter.zeeman.eigs{3}=[-5 -5 10]+5;
inter.zeeman.euler{3}=[0.3 0.1 0.2];
% J-couplings
inter.coupling.scalar=cell(3,3);
inter.coupling.scalar{1,2}=50;
% Coordinates
inter.coordinates={[0.00 0.00 0.00];
[1.05 0.00 0.00];
[0.00 3.50 0.00]};
% Swapping the protons and their correlations
inter.chem.flux_rate=zeros(3,3);
inter.chem.flux_rate(2,3)=1000;
inter.chem.flux_rate(3,2)=1000;
inter.chem.flux_type='intramolecular';
In the general case, the the parameters, supplied to create.m at the spin system specification stage, must be:
| Magnetisation flux parameters | ||
|---|---|---|
| Variable name | Variable type | Content |
| inter.chem.flux_rate | nspins × nspins matrix of real numbers | magnetization transport rates, in Hz |
| inter.chem.flux_type | string: 'intermolecular', 'intramolecular' | controls the fate of correlations with external spins, which are damped in the 'intermolecular' case and transported in the 'intramolecular' case. |
The flux type parameter requires further explanation. When we swap a pair of spins in a molecule, their coherences with any observer spins should in general be preserved. This is what the 'intramolecular' option does. There are cases, however, when the exchange happens through some external medium, such as the solvent, or the end point of the jump is a different molecule. This would lead to the emergence of inter-molecular coherences between identical molecules – a situation that the density matrix formalism is ill-equipped to accommodate. Because such coherences are non-observable, and the spin is very unlikely to jump back to the exact same molecule that it had originally come from, the coherences in question are counted as lost. The corresponding contribution to the kinetics superoperator effectively causes relaxation. This is what happens when the 'intermolecular' option is specified.
Spin-selective radical pair recombination
This exotic process occurs in radical pairs involved in CIDNP, where the two-electron singlet state recombines (and forms diamagnetic products) faster than the two-electron triplet state. Because the rate of the singlet-triplet interconversion depends on the hyperfine couplings present, the result is the sorting of nuclei by orientation (some nuclear spin states accumulate in diamagnetic products, others are left behind in the unreacted radicals).
The kinetic modelling of the process is reviewed in the recent paper by Peter Hore and co-authors (http://dx.doi.org/10.1016/j.cplett.2011.03.082). The exponential model, the Haberkorn model and the Jones-Hore model are currently supported in Spinach. An example of the specification appears below:
% Particles present
sys.isotopes ={'E','E','1H'};
% g-factors and chemical shifts
inter.zeeman.scalar={2.0023 2.0024 1.0};
% Hyperfine couplings
inter.coupling.scalar=cell(3,3);
inter.coupling.scalar{2,3}=1e7;
% Radical pair kinetics model
inter.chem.rp_theory='haberkorn';
% Electron identification
inter.chem.rp_electrons=[1 2];
% Singlet and triplet recombination rates
inter.chem.rp_rates=[1e7 0];
In the general case, the the parameters, supplied to create.m at the spin system specification stage, must be:
| Magnetisation flux parameters | ||
|---|---|---|
| Variable name | Variable type | Content |
| inter.chem.rp_theory | string: 'exponential', 'haberkorn', 'jones-hore' | radical pair recombination theory |
| inter.chem.rp_electrons | 1 x 2 vector of integers | rlectron spins involved in the recombination reaction. |
| inter.chem.rp_rates | 1 x 2 vector of reals | singlet and triplet recombination rates in a radical pair, in Hz |
The theory is only applicable to systems with exactly two unpaired electrons in Liouville space formalisms.
Markov state models
An MSM is essentially a multi-site chemical exchange model involving multiple configurations identified by clustering molecular dynamics simulations. To specify an MSM, do the following:
- Get your molecular dynamics package to output the configurations as separate PDB files.
- Generate the corresponding BMRB files with chemical shifts.
- Import individual configurations using protein.m function.
- Merge individual sys and inter structures that protein.m returns into single large sys and inter structures.
- Specify which spin belongs to which configuration as described in the networks of chemical reactions section above.
- The kinetic matrix described in networks of chemical reactions is exactly the same as your MSM transition rate matrix, just copy it over from your MD package.
That's it - all built-in pulse sequences (for example, noesy.m) will automatically generate and use your kinetics superoperator. To obtain it for some purposes of your own, call kinetics.m function.
Non-linear kinetics
As of version 2.8, Spinach supports non-linear kinetics. The functionality has not been automated or polished yet, but do take a look at react_gen.m and the nonlinear kinetics examples folder if you are interested.
Version 2.8, authors: Ilya Kuprov, Hannah Hogben, Peter Hore
