Extended Bloch–McConnell equations for mechanistic analysis of hyperpolarized ¹³C magnetic resonance experiments on enzyme systems

We describe an approach to formulating the kinetic master equations of the time evolution of NMR signals in reacting (bio)chemical systems. Special focus is given to studies that employ signal enhancement (hyperpolarization) methods such as dissolution dynamic nuclear polarization (dDNP) and involving nuclear spin-bearing solutes that undergo reactions mediated by enzymes and membrane transport proteins. We extend the work given in a recent presentation on this topic (Kuchel and Shishmarev, 2020) to now include enzymes with two or more substrates and various enzyme reaction mechanisms as classified by Cleland, with particular reference to non-first-order processes. Using this approach, we can address some pressing questions in the field from a theoretical standpoint. For example, why does binding of a hyperpolarized substrate to an enzyme not cause an appreciable loss of the signal from the substrate or product? Why does the concentration of an unlabelled pool of substrate, for example ¹²C lactate, cause an increase in the rate of exchange of the ¹³C-labelled pool? To what extent is the equilibrium position of the reaction perturbed during administration of the substrate? The formalism gives a full mechanistic understanding of the time courses derived and is of relevance to ongoing clinical trials using these techniques.

1 Introduction

Nuclear magnetic resonance (NMR) spectroscopy and imaging (MRI) are widely employed techniques with far-reaching applications in physics, chemistry, medicine and the life sciences. NMR and MRI provide a wealth of information from structure elucidation, protein dynamics and metabolic profiling through to disease diagnostics in oncology, cardiology and neurology among others. The technique's low sensitivity is one of the primary concerns in the magnetic resonance community and is often a limiting factor in experiments from solid-state NMR to medical imaging. Recent work has shown that the sensitivity of NMR experiments can be improved by using non-equilibrium hyperpolarization techniques such as dissolution dynamic nuclear polarization (dDNP) to boost signal intensities by many orders of magnitude (Ardenkjaer-Larsen et al., 2003). Such techniques have led to new applications (Golman et al., 2003, 2006; Keshari and Wilson, 2014) and necessitated the development of acquisition strategies to exploit the hyperpolarized magnetization in a time-efficient manner (Yen et al., 2009) as well as new tools for signal processing and image reconstruction (Hu et al., 2010). A challenge with the interpretation of these recordings is that, unlike radio tracers, hyperpolarized MR is a non-tracer technique requiring the injection of physiological or even supra-physiological concentrations of substrate.

To date there have been many mathematical methods devised for analysing the kinetic time courses in dDNP NMR studies (Zierhut et al., 2010; Hill et al., 2013b; Pagès and Kuchel, 2015; Daniels et al., 2016). However, until recently there has been little consensus on the best methods for analysing and then interpreting reaction kinetics measured therein. A theoretical framework has only recently appeared to fully elucidate the underlying mechanisms (Kuchel and Shishmarev, 2020). One challenge is that the widely used Bloch–McConnell equations describe the exchange of magnetization of only the MR active nuclei, while the reaction kinetics are subject to a plethora of molecular interactions in a (bio)chemical milieu. Furthermore, in a typical hyperpolarized MR experiment the initial injection of a non-tracer concentration of substrate causes the reaction system to be perturbed from its equilibrium state, or quasi-steady state, and therefore the concentrations of the reactants are time-dependent. In this regard, challenges relate to the description of non-linear kinetics, for example second-order reactions, and the involvement of unobservable (non-labelled) metabolites to the overall kinetics, e.g. enzyme cofactors, co-substrates and natural abundance ¹²C-containing metabolites (Hill et al., 2013a), as well as explicit descriptions of enzyme mechanisms, e.g. sequential-ordered, sequential-random, double-displacement (ping–pong) reactions, and allosteric interactions that occur on an enzyme far from its active site. Enzyme activity is also influenced by inhibitors that can be competitive, non-competitive, or uncompetitive (Cleland, 1967; Cook and Cleland, 2007). Mathematical models of enzyme systems should agree with standard descriptions of (bio)chemical kinetics while remaining capable of describing the time evolution of magnetization that is described by the Bloch–McConnell equations (McConnell, 1958).

Here we address these issues in a stepwise manner by developing a mechanistic approach that combines the MR interactions with the chemical and/or enzyme-mediated reactions described by the Bloch–McConnell equations. These equations are grounded in the concept of conservation of mass of the species responsible for the hyperpolarized signal plus its non-hyperpolarized counterpart and the various products; this was recently highlighted (Kuchel and Shishmarev, 2020) where the MR-visible signal decays to produce an MR-invisible one.

1.1 Basic concepts – sensitivity

We begin addressing the problem by defining the signal-to-noise ratio (SNR) in MR. In its most basic form, sensitivity is described by the ratio of the signal amplitude divided by the root mean square of the amplitude of the noise. When a signal S(t) is detected in the NMR receiver coil that surrounds the sample, the magnitude of the induced current is a function of (i) the perturbation of nuclear spin populations from thermal equilibrium S_sample(t) plus (ii) a random contribution from the noise in the electronic circuitry S_electronics(t). Hence

$$\begin{array}{}\text{(1)}& \mathrm{S}\left(t\right)={\mathrm{S}}_{\mathrm{sample}}\left(t\right)+{\mathrm{S}}_{\mathrm{electronics}}\left(t\right).\end{array}$$

The current induced in the coil is time-dependent and proportional to the magnetization that precesses in the x,y plane. In other words, the signal S(t) is recorded until decoherence renders S_sample(t) undetectable against the noise, S_electronics(t). The latter is mainly attributed to the radiofrequency (RF) circuitry in the probe head and the preamplifier(s) (e.g. Johnson noise – Johnson, 1928) of the spectrometer. If the NMR signal (free induction decay; FID) that is detected in a subsequent experiment is indistinguishable from the first and the two are added together, then the signal amplitude (peak area) will scale linearly with the number of added FIDs, N. The noise associated with each experiment is random, and assuming its source remains fixed over time, i.e. stationary noise, then the amplitude scales with the square root of the number of FIDs, $N^{\mathrm{1}/\mathrm{2}}$. Hence signal summation enhances the SNR of an NMR experiment in proportion to the square root of the number of FIDs. In other words, to achieve an enhancement by a factor ξ requires an increase in experiment duration of ξ². Therefore, unavoidably, FID summation is a slow process, and experiments can sometimes take days or weeks to achieve a sufficient SNR from a sample of a low-sensitivity nuclide or one with a long relaxation time. The amount of attainable signal averaging is constrained when monitoring dynamic processes by NMR spectroscopy, and an inherently good SNR is required from the outset for a time course experiment.

1.2 Thermal effects

The usual way to proceed when calculating the NMR response of a spin system to RF pulse sequences is to solve the ordinary quantum mechanical master equation that describes the evolution of the spin density operator (Hore et al., 2015). This is the Liouville–von Neumann equation, which has been extended to include non-coherent interactions (predominantly relaxation phenomena) (Ernst et al., 1987):

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\mathit{\rho }=-i\widehat{H}\mathit{\rho }-\widehat{\mathrm{\Gamma }}\left(\mathit{\rho }-{\mathit{\rho }}_{\mathrm{0}}\right),\end{array}$$

where $\widehat{H}$ is the commutation superoperator of the coherent Hamiltonian H given by $\widehat{H}\mathit{\rho }=\left[H,\mathit{\rho }\right]$, which contains information on all spin–spin and field–spin interactions, while $\widehat{\mathrm{\Gamma }}$ is the relaxation superoperator that describes all longitudinal (T₁) and transverse (T₂) relaxation processes, as well as any cross-relaxation or cross-correlation interactions. Note that in the interests of reducing clutter in equations (for which the operator context should be clear), hereafter we have omitted carets denoting operators and only used them to denote superoperators.

Our aim here is to describe the kinetics of exchange between different solutes that contain hyperpolarized nuclei, e.g. A↔B, in which the relaxation times are constant. In this quest, the first simplifying assumption that is worth exploring is that all intermolecular interactions, notably scalar coupling, dipolar coupling, cross-relaxation and cross-correlation between species A and B, can be ignored. This applies to non-interacting solute molecules in solution in which motional averaging occurs, and we focus on thermal effects on the evolution of the FID.

The so-called Zeeman polarization term describes the sensitivity of S_sample(t) in Eq. (1) to temperature and magnetic field in an NMR experiment. Magnetic polarization is described by the equilibrium density operator ρ₀ that specifies the probability distribution of states. Zeeman polarization corresponds to the magnitude of normalized longitudinal spin order I_z that is contained in ρ₀. Specifically, for an ensemble of spin-$\frac{\mathrm{1}}{\mathrm{2}}$ nuclei this is given by (Ernst et al., 1987)

$$\begin{array}{}\text{(3)}& {\mathit{\rho }}_{\mathrm{0}}={\displaystyle \frac{\exp (-\mathrm{\hslash }{H}_{\mathrm{0}}/kT)}{\mathrm{Tr}\mathit{\left\{}\exp \right(-\mathrm{\hslash }{H}_{\mathrm{0}}/kT\left)\mathit{\right\}}}},\end{array}$$

where k is the Boltzmann constant and T is the temperature (Kelvin). The Zeeman Hamiltonian H₀ describes the interaction of the spins with the static magnetic field of magnitude B₀, given by H₀=ω₀I_z, where ω₀ is the Larmor frequency (rad s⁻¹). In the basis of the two eigenstates |α〉 (“spin-up”) and |β〉 (“spin-down”), the equilibrium density operator is written in matrix form as

$$\begin{array}{}\text{(4)}& {\mathit{\rho }}_{\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{Z}}\left[\begin{array}{cc}\exp (\mathrm{\hslash }{\mathit{\omega }}_{\mathrm{0}}/\mathrm{2}kT)& \mathrm{0}\\ \mathrm{0}& \exp (-\mathrm{\hslash }{\mathit{\omega }}_{\mathrm{0}}/\mathrm{2}kT)\end{array}\right],\end{array}$$

where Z is the partition function, given by $Z={\sum }_{i=\mathrm{1}}^M\exp \left(-{\mathit{\epsilon }}_i/kT\right)$, and M is the number of states (M=2 for an $I=\frac{\mathrm{1}}{\mathrm{2}}$ nucleus). In the case of a spin-$\frac{\mathrm{1}}{\mathrm{2}}$ system, the partition function is the sum of the populations $Z=\exp (\mathrm{\hslash }{\mathit{\omega }}_{\mathrm{0}}/\mathrm{2}kT)+\exp (-\mathrm{\hslash }{\mathit{\omega }}_{\mathrm{0}}/\mathrm{2}kT)\approx \mathrm{2}$ when ε_i is very small, as is typically the case at thermal equilibrium in NMR systems. The Zeeman polarization is proportional to the projection of the spin density operator onto the angular momentum operator. In other words, it is proportional to the expectation value of 〈I_z〉 and is given by (Keeler, 2010)

$$\begin{array}{}\text{(5)}& 〈I_z〉=\mathrm{Tr}\left[{\mathit{\rho }}_{\mathrm{0}}{I}_z\right]={\displaystyle \frac{\mathrm{1}}{\mathrm{2}Z}}\left[\exp \right(\mathrm{\hslash }{\mathit{\omega }}_{\mathrm{0}}/\mathrm{2}kT)-\exp (-\mathrm{\hslash }{\mathit{\omega }}_{\mathrm{0}}/\mathrm{2}kT\left)\right].\end{array}$$

Hence, the Zeeman polarization for an ensemble of nuclear spins is the normalized imbalance between the populations of the |α〉 and |β〉 states, p_α and p_β, respectively; in other words, it is the normalized net population difference that is given by

$$\begin{array}{}\text{(6)}& P={\displaystyle \frac{p_{\mathit{\alpha }}-p_{\mathit{\beta }}}{p_{\mathit{\alpha }}+p_{\mathit{\beta }}}}.\end{array}$$

This normalization is carried out with respect to the total population of the nuclear ensemble such that $p_{\mathit{\alpha }}+p_{\mathit{\beta }}=\mathrm{1}$. Therefore, the bounds on the polarization are $-\mathrm{1}<P<+\mathrm{1}$. At room temperature (∼298 K) and in a field of 11.75 T (500 MHz for ¹H nuclei), the thermal equilibrium Zeeman polarization, P_z,eq, is a mere $\sim \mathrm{4}\times {\mathrm{10}}^{-\mathrm{5}}$. Thus, there is only a tiny population difference between the spin states of a nuclear ensemble, which implies inherently weak polarization. It is this small population imbalance which is manipulated in NMR experiments under thermal equilibrium conditions. This weak polarization is a consequence of the small difference in energy (∼0.1 J mol⁻¹) between nuclear spin energy levels at room temperature (∼2.5 kJ mol⁻¹), and it implies only weak alignment of nuclear spins in the static magnetic field of all contemporary superconducting magnets.

In the usual quantum mechanical analysis of multiple spin systems, the density operator (that describes the probability density of states) is normalized to 1, meaning that the summed (total) probability density of all states is 1. This is expressed mathematically as Tr[ρ]=1, where Tr denotes the trace of the matrix (Hore et al., 2015). To describe non-equilibrium reactions in terms of solute concentrations requires a scaled density operator (Kuhne et al., 1979):

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_i=\left[{\mathrm{A}}_i\right]{\mathit{\rho }}_i,\end{array}$$

where σ_i is now proportional to concentration [A_i]. Differentiation of Eq. (7) leads to

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{d}{\mathit{\sigma }}_i}{\mathrm{d}t}}=\left[{\mathrm{A}}_i\right]{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_i}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}\left[{\mathrm{A}}_i\right]}{\mathrm{d}t}}{\mathit{\rho }}_i.\end{array}$$

Therefore, it follows that for a system at chemical equilibrium $\mathrm{d}\left[{\mathrm{A}}_i\right]/\mathrm{d}t=\mathrm{0}$, the scaled density operator is directly proportional to the normalized density operator. For non-equilibrium systems the concentrations are time-dependent viz. $\mathrm{d}\left[{\mathrm{A}}_i\right]/\mathrm{d}t\ne \mathrm{0}$, so the two no longer scale in a straightforward manner.

On the other hand, equilibrium magnetization (M_z,eq) is a bulk property that is the net magnetic dipole moment per unit volume and is proportional to 〈I_z〉, where the proportionality factor is Nℏγ. From Eq. (5) this yields the expression for the magnetization in terms of magnetic field strength, temperature and number of spins in the sample (or more specifically in the detection volume of the NMR spectrometer):

$$\begin{array}{}\text{(9)}& M_{z,\mathrm{eq}}={\displaystyle \frac{N\mathrm{\hslash }\mathit{\gamma }}{\mathrm{2}}}\tanh \left({\displaystyle \frac{\mathrm{\hslash }\mathit{\gamma }{B}_{\mathrm{0}}}{\mathrm{2}kT}}\right).\end{array}$$

In the so-called “high-temperature limit” (room temperature, in the cases addressed here), Eq. (9) simplifies to

$$\begin{array}{}\text{(10)}& M_{z,\mathrm{eq}}={\displaystyle \frac{N{\mathrm{\hslash }}^{\mathrm{2}}{\mathit{\gamma }}^{\mathrm{2}}{B}_{\mathrm{0}}}{\mathrm{4}kT}}.\end{array}$$

In words, “thermal magnetization” is proportional to the magnitude of the external magnetic field strength, B₀, and is inversely proportional to the temperature, T, while being proportional to the number of spins, N. Therefore, it is proportional to the concentration [A_i] of the solute that bears the NMR-active nucleus.

2 Equation of motion – the Bloch equations

In the absence of intermolecular binding (however transient) or scalar couplings, the motion (time evolution) of magnetizations is described by the Bloch equations. Magnetization is explicitly declared to be proportional to reactant concentrations [A] and [B], as has recently been discussed (Kuchel and Shishmarev, 2020). To explore this situation, we start with the basic Bloch equations for a single spin-$\frac{\mathrm{1}}{\mathrm{2}}$ ensemble. The equation describes the time evolution of x,y and z magnetization in the rotating frame and includes the influence of chemical shift, RF fields, and transverse (T₂) and longitudinal relaxation (T₁) time constants. The Bloch equations in their complete form are described as being inhomogeneous, and they can be written using a matrix and vectors:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}{M}_x\\ M_y\\ M_z\end{array}\right]& =-\left[\begin{array}{ccc}{R}_{\mathrm{2}}& \mathrm{\Omega }& -{\mathit{\omega }}_y\\ -\mathrm{\Omega }& R_{\mathrm{2}}& {\mathit{\omega }}_x\\ {\mathit{\omega }}_y& -{\mathit{\omega }}_x& R_{\mathrm{1}}\end{array}\right]\left[\begin{array}{c}{M}_x\\ M_y\\ M_z\end{array}\right]\\ & +\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ R_{\mathrm{1}}{M}_{z,\mathrm{eq}}\end{array}\right],\end{array}\end{array}$$

where $\mathrm{\Omega }={\mathit{\omega }}_{\mathrm{0}}-{\mathit{\omega }}_{\mathrm{RF}}$ is the “offset frequency” in the rotating frame; ω₀ (rad s⁻¹) is the Larmor frequency; ω_RF (rad s⁻¹) is the RF frequency; the x component of the RF field (rad s⁻¹) is ${\mathit{\omega }}_x=-\mathit{\gamma }{B}_{\mathrm{1}}\cos ({\mathit{\omega }}_{\mathrm{RF}}t+\mathit{\phi }$); and the y component is ${\mathit{\omega }}_y=-\mathit{\gamma }{B}_{\mathrm{1}}\sin ({\mathit{\omega }}_{\mathrm{RF}}t+\mathit{\phi }$), where the magnitude of the field strength is B₁ and the phase of the wave form relative to an internal reference source is φ. The longitudinal relaxation rate constant is denoted by $R_{\mathrm{1}}=\mathrm{1}/T_{\mathrm{1}}$, the transverse one by $R_{\mathrm{2}}=\mathrm{1}/T_{\mathrm{2}}$, and the equilibrium magnetization by M_z,eq.

Equation (11) is tedious to solve analytically, but it is readily solved numerically (Allard et al., 1998; Helgstrand et al., 2000). On the other hand, by including the identity operator in the basis set and adding a constant to the equilibrium magnetization (Levitt and Dibari, 1992), we obtain a much more compliant (to analysis) matrix equation:

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}\frac{E}{\mathrm{2}}\\ M_x\\ M_y\\ M_z\end{array}\right]=-\left[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& R_{\mathrm{2}}& \mathrm{\Omega }& -{\mathit{\omega }}_y\\ \mathrm{0}& -\mathrm{\Omega }& R_{\mathrm{2}}& {\mathit{\omega }}_x\\ -\mathrm{2}\mathrm{\Theta }& {\mathit{\omega }}_y& -{\mathit{\omega }}_x& R_{\mathrm{1}}\end{array}\right]\left[\begin{array}{c}\frac{E}{\mathrm{2}}\\ M_x\\ M_y\\ M_z\end{array}\right],\end{array}$$

where E is equal to 1 and the factor $\mathrm{\Theta }=R_{\mathrm{1}}{M}_{z,\mathrm{eq}}$ describes the equilibrium magnetization.

2.1 Chemical exchange kinetics of systems prior to and at equilibrium – the Bloch–McConnell equations

We can extend the system of equations from describing an ensemble of single spins to two or more exchanging spins. The system of equations now accounts for the magnetization interaction with the lattice and exchange via the forward and reverse chemical reactions. These are the Bloch–McConnell equations (McConnell, 1958).

First, consider the rate expressions for a simple bi-directional chemical reaction. The coupled differential equations describing first-order reaction kinetics of solute A becoming solute B and back again, A↔B, are typically expressed in terms of molar concentrations:

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\mathrm{d}\left[\mathrm{A}\left(t\right)\right]}{\mathrm{d}t}}=-k_{\mathrm{1}}\left[\mathrm{A}\left(t\right)\right]+k_{-\mathrm{1}}\left[\mathrm{B}\left(t\right)\right],\text{(14)}& {\displaystyle \frac{\mathrm{d}\left[\mathrm{B}\left(t\right)\right]}{\mathrm{d}t}}=k_{\mathrm{1}}\left[\mathrm{A}\left(t\right)\right]-k_{-\mathrm{1}}\left[\mathrm{B}\left(t\right)\right],\end{array}$$

which can be expressed in matrix form:

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}\left[\mathrm{A}\left(t\right)\right]\\ \left[\mathrm{B}\left(t\right)\right]\end{array}\right]=\left[\begin{array}{cc}-k_{\mathrm{1}}& k_{-\mathrm{1}}\\ k_{\mathrm{1}}& -k_{-\mathrm{1}}\end{array}\right]\left[\begin{array}{c}\left[\mathrm{A}\left(t\right)\right]\\ \left[\mathrm{B}\left(t\right)\right]\end{array}\right].\end{array}$$

The rate constant for the forward reaction is denoted by k₁, while for the reverse reaction it is k₋₁. The time-dependent concentrations are given by [A(t)] and [B(t)]. As required by the fact that this is a closed system, the equations must conform to the principle of conservation of mass. Specifically, the sum of the rates of change in [A(t)] and [B(t)] given by $\mathrm{d}\left[\mathrm{A}\left(t\right)\right]/\mathrm{d}t+\mathrm{d}\left[\mathrm{B}\left(t\right)\right]/\mathrm{d}t$ is zero. We return to this point below. In other words, mass is neither created nor destroyed during the reaction in such a closed system.

For the simplest case of two magnetically active solutes, each possessing a single spin-$\frac{\mathrm{1}}{\mathrm{2}}$ nuclide, in chemical exchange, A↔B, the direct product (a mathematical operation used in quantum mechanics to generate the necessary combinations of states) of the chemical (solute) space {[A],[B]} and the magnetization vector space $\left\{M_x,M_y,M_z\right\}$ for each of A and B is given by

$$\begin{array}{}\text{(16)}& \left[\begin{array}{c}\mathrm{1}\\ \mathrm{1}\end{array}\right]\otimes \left[\begin{array}{c}{M}_x\\ M_y\\ M_z\end{array}\right]=\left[\begin{array}{c}{M}_x^{\mathrm{A}}\\ M_y^{\mathrm{A}}\\ M_z^{\mathrm{A}}\\ M_x^{\mathrm{B}}\\ M_y^{\mathrm{B}}\\ M_z^{\mathrm{B}}\end{array}\right].\end{array}$$

A new exchange matrix in the basis of the new magnetization space $\left\{M_x^{\mathrm{A}},M_y^{\mathrm{A}},M_z^{\mathrm{A}},M_x^{\mathrm{B}},M_y^{\mathrm{B}},M_z^{\mathrm{B}}\right\}$ is calculated by taking the direct product of the exchange matrix with the identity operator I that is chosen to have the same dimensions as the magnetization space. The direct product is given by

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& \left[\begin{array}{cc}-k_{\mathrm{1}}& k_{-\mathrm{1}}\\ k_{\mathrm{1}}& -k_{-\mathrm{1}}\end{array}\right]\otimes \left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]=\\ & \left[\begin{array}{cccccc}-k_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -k_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& k_{-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& k_{-\mathrm{1}}\\ k_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& -k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& k_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& -k_{-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& k_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& -k_{-\mathrm{1}}\end{array}\right].\end{array}\end{array}$$

Likewise, the matrix describing coherent and incoherent magnetization interactions can be recast in a similar fashion to give

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& \left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}\end{array}\right]\otimes \left[\begin{array}{ccc}{R}_{\mathrm{2}}& \mathrm{\Omega }& -{\mathit{\omega }}_y\\ -\mathrm{\Omega }& R_{\mathrm{2}}& {\mathit{\omega }}_x\\ {\mathit{\omega }}_y& -{\mathit{\omega }}_x& R_{\mathrm{1}}\end{array}\right]\\ & =\left[\begin{array}{cccccc}{R}_{\mathrm{2}}^{\mathrm{A}}& {\mathrm{\Omega }}^{\mathrm{A}}& -{\mathit{\omega }}_y& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -{\mathrm{\Omega }}^{\mathrm{A}}& R_{\mathrm{2}}^{\mathrm{A}}& {\mathit{\omega }}_x& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\mathit{\omega }}_y& -{\mathit{\omega }}_x& R_{\mathrm{1}}^{\mathrm{A}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& R_{\mathrm{2}}^{\mathrm{B}}& {\mathrm{\Omega }}^{\mathrm{B}}& -{\mathit{\omega }}_y\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -{\mathrm{\Omega }}^{\mathrm{B}}& R_{\mathrm{2}}^{\mathrm{B}}& {\mathit{\omega }}_x\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{\omega }}_y& -{\mathit{\omega }}_x& R_{\mathrm{1}}^{\mathrm{B}}\end{array}\right].\end{array}\end{array}$$

The inhomogeneous form of the Bloch equations can now be constructed to take into account both the coherent and incoherent interactions as well as chemical exchange. This yields the inhomogeneous form of the Bloch–McConnell equations, which are written (again in matrix form) as

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}{M}_x^{\mathrm{A}}\\ M_y^{\mathrm{A}}\\ M_z^{\mathrm{A}}\\ M_x^{\mathrm{B}}\\ M_y^{\mathrm{B}}\\ M_z^{\mathrm{B}}\end{array}\right]=\\ & -\left[\begin{array}{cccccc}{R}_{\mathrm{2}}^{\mathrm{A}}+k_{\mathrm{1}}& {\mathrm{\Omega }}^{\mathrm{A}}& -{\mathit{\omega }}_y& -k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ -{\mathrm{\Omega }}^{\mathrm{A}}& R_{\mathrm{2}}^{\mathrm{A}}+k_{\mathrm{1}}& {\mathit{\omega }}_x& \mathrm{0}& -k_{-\mathrm{1}}& \mathrm{0}\\ {\mathit{\omega }}_y& -{\mathit{\omega }}_x& R_{\mathrm{1}}^{\mathrm{A}}+k_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& -k_{-\mathrm{1}}\\ -k_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& R_{\mathrm{2}}^{\mathrm{B}}+k_{-\mathrm{1}}& {\mathrm{\Omega }}^{\mathrm{B}}& -{\mathit{\omega }}_y\\ \mathrm{0}& -k_{\mathrm{1}}& \mathrm{0}& -{\mathrm{\Omega }}^{\mathrm{B}}& R_{\mathrm{2}}^{\mathrm{B}}+k_{-\mathrm{1}}& {\mathit{\omega }}_x\\ \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}& {\mathit{\omega }}_y& -{\mathit{\omega }}_x& R_{\mathrm{1}}^{\mathrm{B}}+k_{-\mathrm{1}}\end{array}\right]\\ & \left[\begin{array}{c}{M}_x^{\mathrm{A}}\\ M_y^{\mathrm{A}}\\ M_z^{\mathrm{A}}\\ M_x^{\mathrm{B}}\\ M_y^{\mathrm{B}}\\ M_z^{\mathrm{B}}\end{array}\right]+\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ R_{\mathrm{1}}^{\mathrm{A}}{M}_{z,\mathrm{eq}}^{\mathrm{A}}\\ \mathrm{0}\\ \mathrm{0}\\ R_{\mathrm{1}}^{\mathrm{B}}{M}_{z,\mathrm{eq}}^{\mathrm{B}}\end{array}\right],\end{array}\end{array}$$

where $M_{z,\mathrm{eq}}^{\mathrm{A}}$ and $M_{z,\mathrm{eq}}^{\mathrm{B}}$ denote the respective equilibrium magnetizations (hence the subscript eq).

The inhomogeneous form of the Bloch–McConnell equations can similarly be modified by incorporating the equilibrium magnetization to create a homogeneous form of this master equation:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}\frac{E}{\mathrm{2}}\\ M_x^{\mathrm{A}}\\ M_y^{\mathrm{A}}\\ M_z^{\mathrm{A}}\\ M_x^{\mathrm{B}}\\ M_y^{\mathrm{B}}\\ M_z^{\mathrm{B}}\end{array}\right]=\\ & -\left[\begin{array}{ccccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& R_{\mathrm{2}}^{\mathrm{A}}+k_{\mathrm{1}}& {\mathrm{\Omega }}^{\mathrm{A}}& -{\mathit{\omega }}_y& -k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -{\mathrm{\Omega }}^{\mathrm{A}}& R_{\mathrm{2}}^{\mathrm{A}}+k_{\mathrm{1}}& {\mathit{\omega }}_x& \mathrm{0}& -k_{-\mathrm{1}}& \mathrm{0}\\ -\mathrm{2}{\mathrm{\Theta }}^{\mathrm{A}}& {\mathit{\omega }}_y& -{\mathit{\omega }}_x& R_{\mathrm{1}}^{\mathrm{A}}+k_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& -k_{-\mathrm{1}}\\ \mathrm{0}& -k_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& R_{\mathrm{2}}^{\mathrm{B}}+k_{-\mathrm{1}}& {\mathrm{\Omega }}^{\mathrm{B}}& -{\mathit{\omega }}_y\\ \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}& \mathrm{0}& -{\mathrm{\Omega }}^{\mathrm{B}}& R_{\mathrm{2}}^{\mathrm{B}}+k_{-\mathrm{1}}& {\mathit{\omega }}_x\\ -\mathrm{2}{\mathrm{\Theta }}^{\mathrm{B}}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}& {\mathit{\omega }}_y& -{\mathit{\omega }}_x& R_{\mathrm{1}}^{\mathrm{B}}+k_{-\mathrm{1}}\end{array}\right]\\ & \left[\begin{array}{c}\frac{E}{\mathrm{2}}\\ M_x^{\mathrm{A}}\\ M_y^{\mathrm{A}}\\ M_z^{\mathrm{A}}\\ M_x^{\mathrm{B}}\\ M_y^{\mathrm{B}}\\ M_z^{\mathrm{B}}\end{array}\right].\end{array}\end{array}$$

Again, the factors ${\mathrm{\Theta }}^{\mathrm{A}}=R_{\mathrm{1}}^{\mathrm{A}}{M}_{z,\mathrm{eq}}^{\mathrm{A}}$ and ${\mathrm{\Theta }}^{\mathrm{B}}=R_{\mathrm{1}}^{\mathrm{B}}{M}_{z,\mathrm{eq}}^{\mathrm{B}}$ account for the respective equilibrium magnetizations.

2.1.1 Simulations of thermal kinetics using Eq. (19)

Next, consider Eq. (19) for simulating the evolution of the x,y, and z components of the magnetization of a “thermal magnetization” (non-hyperpolarized) sample. We seek the NMR spectrum that results from a two-site exchange reaction between solutes A and B, Fig. 1a, as conventionally observed in room temperature NMR experiments.

Figure 1 Simulated NMR spectra resulting from a two-site exchange process between thermally polarized solutes, A↔B, shown schematically in (a). Simulated FIDs S(t) are shown in (b)–(e), left panel, with corresponding spectra s(ω), middle panel, and the recovery of z magnetizations, $M_z^{\mathrm{A}}\left(t\right)$ and $M_z^{\mathrm{B}}\left(t\right)$, right panel. Spectra were simulated with rate constants, (b) $k_{\mathrm{1}}=k_{-\mathrm{1}}=\mathrm{0}$; (c) k₁=2 s⁻¹, $k_{-\mathrm{1}}=\mathrm{1}$ s⁻¹; (d) k₁=20 s⁻¹, $k_{-\mathrm{1}}=\mathrm{10}$ s⁻¹; and (e) k₁=2000 s⁻¹, $k_{-\mathrm{1}}=\mathrm{1000}$ s⁻¹, corresponding to no exchange, slow, intermediate, and fast exchange regimes, respectively.

Simulations were performed in MatLab with equilibrium z magnetizations $M_{z,\mathrm{eq}}^{\mathrm{A}}=\mathrm{1.0}$ and $M_{z,\mathrm{eq}}^{\mathrm{B}}=\mathrm{0.8}$ and an initial magnetization vector given by ${\mathit{M}}_{\mathrm{0}}=\left[\mathrm{0},\mathrm{0},\mathrm{1.0},\mathrm{0},\mathrm{0},\mathrm{0.8}\right]$. Chemical shift offsets were ${\mathrm{\Omega }}^{\mathrm{A}}=\mathrm{10}\times \mathrm{2}\mathit{\pi }$ rad s⁻¹ and ${\mathrm{\Omega }}^{\mathrm{B}}=-\mathrm{10}\times \mathrm{2}\mathit{\pi }$ rad s⁻¹. Relaxation rate constants were $R_{\mathrm{1}}^{\mathrm{A}}=R_{\mathrm{1}}^{\mathrm{B}}=\mathrm{1}$ s⁻¹ and $R_{\mathrm{2}}^{\mathrm{A}}=R_{\mathrm{2}}^{\mathrm{B}}=\mathrm{1}$ s⁻¹. The influence of an RF_y pulse was then calculated with ${\mathit{\omega }}_x=-\mathit{\gamma }{B}_{\mathrm{1}}\cos \left(\mathit{\pi }/\mathrm{2}\right)$ and ${\mathit{\omega }}_y=-\mathit{\gamma }{B}_{\mathrm{1}}\sin \left(\mathit{\pi }/\mathrm{2}\right)$ and with a field strength of 1.5 kHz, corresponding to ${\mathit{\omega }}_y=-\mathit{\gamma }{B}_{\mathrm{1}}=-\mathrm{1500}\times \mathrm{2}\mathit{\pi }$ rad s⁻¹ and ω_x=0. For a 90^∘ RF nutation (flip) angle the pulse duration is $t_{\mathrm{p}}=\mathit{\pi }/\mathrm{2}{\mathit{\omega }}_y$, which gave a transformed magnetization vector after the pulse of $\mathit{M}\left(t\right)=\left[\mathrm{0.999},\mathrm{0.007},\mathrm{0.000},\mathrm{0.800},-\mathrm{0.005},\mathrm{0.000}\right]$; this was composed mostly of $M_x^{\mathrm{A}}+M_x^{\mathrm{B}}$ with a residual contribution from $M_y^{\mathrm{A}}+M_y^{\mathrm{B}}$ arising from evolution of the chemical shift during the RF pulse and a small contribution from $M_z^{\mathrm{A}}+M_z^{\mathrm{B}}$ due to return of the magnetization to the equilibrium state.

The observable signal (the FID, which is a function of time) is proportional to the complex signal $\mathrm{S}\left(t\right)=M_x^{\mathrm{A}}\left(t\right)-i{M}_y^{\mathrm{A}}\left(t\right)+M_x^{\mathrm{B}}\left(t\right)-i{M}_y^{\mathrm{B}}\left(t\right)$. Noise was simulated by adding to the FID a normally distributed complex random vector with mean = 0 and standard deviation (SD) = 0.1. The spectrum s(ω) was then calculated by taking the Fourier transform of S(t). Simulated FIDs, S(t), are shown in Fig. 1b–e left panel, the corresponding spectra s(ω) in Fig. 1b–e middle panel, and the recovery of the z magnetizations $M_z^{\mathrm{A}}\left(t\right)$ and $M_z^{\mathrm{B}}\left(t\right)$ are shown in Fig. 1b–e, right panel. Spectra were simulated for a range of rate constants, where exchange was either absent, $k_{\mathrm{1}}=k_{-\mathrm{1}}=\mathrm{0}$, Fig. 1b, or for increasing rates of exchange. Thus, (c) k₁=2 s⁻¹, $k_{-\mathrm{1}}=\mathrm{1}$ s⁻¹; (d) k₁=20 s⁻¹, $k_{-\mathrm{1}}=\mathrm{10}$ s⁻¹; and (e) k₁=2000 s⁻¹, $k_{-\mathrm{1}}=\mathrm{1000}$ s⁻¹, corresponding to the slow, intermediate and fast regimes, respectively.

The equilibrium constant was fixed so that $K=k_{\mathrm{1}}/k_{-\mathrm{1}}=\mathrm{2}$; hence the system was not at chemical equilibrium at t=0 s. The simulations highlight an important point: in the absence of exchange the Bloch–McConnell equations predict the recovery of the z magnetizations back to their magnetic equilibrium values $M_{z,\mathrm{eq}}^{\mathrm{A}}$ and $M_{z,\mathrm{eq}}^{\mathrm{B}}$, while under conditions of fast exchange this no longer takes place during the experiment. A non-equilibrium system will rapidly recover to its chemical equilibrium but not to its initial thermal equilibrium $M_{z,\mathrm{eq}}^{\mathrm{A}}$ and $M_{z,\mathrm{eq}}^{\mathrm{B}}$; again, in other words, this does not take place within the timescale of the NMR experiment, which is typically within five T₁ values.

2.2 Describing hyperpolarized kinetics with the Bloch–McConnell equations

We now consider the predictions made by using Eq. (19) when simulating the evolution of the x, y, and z components of the magnetization of a hyperpolarized sample and the resulting spectrum for a two-site exchange reaction between solutes A and B. In the previous example the initial condition was $M_z^{\mathrm{A}}\left(\mathrm{0}\right)=\mathrm{1.0}$ and $M_z^{\mathrm{B}}\left(\mathrm{0}\right)=\mathrm{0.8}$. To extend the Bloch–McConnell formalism to be able to predict the dynamics of a hyperpolarized experiment, we recognize that for the same magnitude of noise in the receiver circuit (although this may not be true for a hyperpolarized sample) the initial hyperpolarized magnetization is given by

$$\begin{array}{}\text{(21)}& M_{z,\mathrm{hyp}}=\mathit{\eta }{M}_{z,\mathrm{eq}},\end{array}$$

where η is the enhancement factor that varies from one hyperpolarization experiment to another. In the case of dDNP experiments η≅10⁴ is typical, although this depends on the method of hyperpolarization, the solute(s) in question and a set of physicochemical parameters that are described in detail in e.g. Ardenkjaer-Larsen et al. (2015).

Simulations of hyperpolarized kinetics using Eq. (19)

These were performed with equilibrium z magnetizations $M_{z,\mathrm{eq}}^{\mathrm{A}}=\mathrm{1.0}$ and $M_{z,\mathrm{eq}}^{\mathrm{B}}=\mathrm{0.8}$, as used above, but now with an initial magnetization vector $\mathit{M}\left(\mathrm{0}\right)=\left[\mathrm{0},\mathrm{0},\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{4}},\mathrm{0},\mathrm{0},\mathrm{0}\right]$. This situation corresponds to an initial hyperpolarized magnetization $M_{z,\mathrm{hyp}}^{\mathrm{A}}\left(\mathrm{0}\right)$ of only solute A and an enhancement factor of η=10⁴. Chemical shifts were ${\mathrm{\Omega }}^{\mathrm{A}}=\mathrm{10}\times \mathrm{2}\mathit{\pi }$ rad s⁻¹ and ${\mathrm{\Omega }}^{\mathrm{B}}=-\mathrm{10}\times \mathrm{2}\mathit{\pi }$ rad s⁻¹, while relaxation times were increased to represent a hyperpolarized ¹³C substrate, $R_{\mathrm{1}}^{\mathrm{A}}=R_{\mathrm{1}}^{\mathrm{B}}=\mathrm{1}/\mathrm{60}$ s⁻¹ and $R_{\mathrm{2}}^{\mathrm{A}}=R_{\mathrm{2}}^{\mathrm{B}}=\mathrm{1}$ s⁻¹, with the rate constants representing an enzyme-mediated cell reaction $k_{\mathrm{1}}=k_{-\mathrm{1}}=\mathrm{0.005}$ s⁻¹. Figure 2a shows the time evolution of the z components of the magnetization, displaying the familiar (Day et al., 2007) bi-exponential time dependence of $M_{z,\mathrm{hyp}}^{\mathrm{A}}\left(t\right)$ and $M_{z,\mathrm{hyp}}^{\mathrm{B}}\left(t\right)$ magnetizations.

Figure 2 (a) Simulated evolution of the z components of the magnetization $M_z^{\mathrm{A}}$ and $M_z^{\mathrm{B}}$ for a hyperpolarized solute $M_z^{\mathrm{A}}\left(\mathrm{0}\right)=\mathrm{1}\times {\mathrm{10}}^{\mathrm{4}}$ undergoing a two-site exchange reaction, A↔B. Longitudinal relaxation rate constants were $R_{\mathrm{1}}^{\mathrm{A}}=R_{\mathrm{1}}^{\mathrm{B}}=\mathrm{1}/\mathrm{60}$ s⁻¹ and $R_{\mathrm{2}}^{\mathrm{A}}=R_{\mathrm{2}}^{\mathrm{B}}=\mathrm{1}$ s⁻¹. Rate constants were $k_{\mathrm{1}}=k_{-\mathrm{1}}=\mathrm{0.005}$ s⁻¹. (b) Simple pulse sequence for acquiring a time course experiment with multiple sampling of the magnetization and acquisition of an FID at each time point. (c–d) Waterfall plots of simulated spectra resulting from sequential application of the pulse sequence in (b) for an initial hyperpolarized solute A undergoing two-site exchange with solute B, calculated with a flip angle: (c) β=1^∘; and (d) β=20^∘.

We next simulate the effect of applying the pulse sequence shown in Fig. 2b corresponding to a time course type of experiment with multiple sampling of the magnetization and acquisition of an FID at each time point. This is representative of real experiments that have been presented in the literature (Gabellieri et al., 2008; Hill et al., 2013b). The time delays correspond to a pre-scan delay τ, the duration of the pulse t_p and the duration of the FID t_aq. The experiment is repeated n times to sample the entire time course where the temporal resolution is then given by the total repetition time $\mathrm{TR}=\mathit{\tau }+t_{\mathrm{p}}+t_{\mathrm{aq}}$, and the total duration of the experiment is given by nTR. In this experiment we make the assumption that the transverse magnetization from one experiment to the next is not recovered by the application of a subsequent pulse. This assumption is reasonable provided the acquisition time is much longer that the time taken for the FID to decay to zero, namely $t_{\mathrm{aq}}\gg T_{\mathrm{2}}^{\ast }$.

The influence of this pulse sequence was then calculated, accounting for multiple sampling of the magnetization. The RF pulse was again specified by ${\mathit{\omega }}_x=-\mathit{\gamma }{B}_{\mathrm{1}}\cos \left(\mathit{\pi }/\mathrm{2}\right)$ and ${\mathit{\omega }}_y=-\mathit{\gamma }{B}_{\mathrm{1}}\sin \left(\mathit{\pi }/\mathrm{2}\right)$ with a field strength of 1.5 kHz, which corresponds to ${\mathit{\omega }}_y=-\mathit{\gamma }{B}_{\mathrm{1}}=-\mathrm{1500}\times \mathrm{2}\mathit{\pi }$ rad s⁻¹. Application of an RF pulse tilts the hyperpolarized magnetization away from the z axis by an angle of β radians. The magnitude of the observable transverse magnetization is proportional to sin(β), and the remaining longitudinal magnetization is proportional to cos(β).

Simulations were performed with the same magnitude of noise as in Fig. 1. The time evolution of the magnetization was recorded for the pulse sequence shown in Fig. 2b with sequential acquisition of 64 spectra and a repetition time of TR=4.25 s. The effect of acquiring a time series of spectra with either a flip angle β=1^∘, Fig. 2c, or β=20^∘, Fig. 2d, are seen in the stack plots. The pulse length (duration) was $t_{\mathrm{p}}=\mathit{\beta }\mathit{\pi }/\mathrm{180}{\mathit{\omega }}_y$. After a single β=1^∘ pulse was applied to M(0), the magnetization vector was tilted to become $\mathit{M}\left(t\right)=\left[\mathrm{0.174},\mathrm{0.000},\mathrm{9.998},\mathrm{0.000},\mathrm{0.000},\mathrm{0.000}\right]\times {\mathrm{10}}^{\mathrm{3}}$ prior to acquisition of the FID. This was composed mostly of $M_z^{\mathrm{A}}$ with a small contribution from $M_x^{\mathrm{A}}$ that arose from excitation by the β=1^∘ pulse or, following a β=20^∘ pulse the magnetization vector was tilted to become $\mathit{M}\left(t\right)=\left[\mathrm{3.420},\mathrm{0.004},\mathrm{9.397},\mathrm{0.000},\mathrm{0.000},\mathrm{0.000}\right]\times {\mathrm{10}}^{\mathrm{3}}$, again composed mostly of $M_z^{\mathrm{A}}$ but with a greater contribution from $M_x^{\mathrm{A}}$ due to excitation by a pulse with a larger value of β. Since the magnetization relaxed to its thermal equilibrium state, the hyperpolarized magnetization was effectively destroyed during application of the RF (sampling) pulse, and it was not re-generated. This may not be the outcome when non-linear effects such as radiation damping cause recovery of the hyperpolarized signal (Weber et al., 2019).

The z magnetization after the application of a single RF pulse and delay TR is therefore given by

$$\begin{array}{}\text{(22)}& \mathrm{S}\left(\mathrm{TR}\right)=\mathrm{S}\left(\mathrm{0}\right)\cos \left(\mathit{\theta }\right)\exp \left(-R_{\mathrm{1}}\mathrm{TR}\right).\end{array}$$

Following the application of a series of n RF pulses with a total delay nTR=t, the signal is given by (Kuchel and Shishmarev, 2020)

$$\begin{array}{}\text{(23)}& \mathrm{S}\left(t\right)=\mathrm{S}\left(\mathrm{0}\right){\cos }^n\left(\mathit{\theta }\right)\exp \left(-R_{\mathrm{1}}t\right).\end{array}$$

The apparent relaxation time constant of the hyperpolarized signal, including the influence of both the intrinsic T₁ and flip angle correction, is given by Hill et al. (2013b) and Kuchel and Shishmarev (2020):

$$\begin{array}{}\text{(24)}& {\displaystyle }\exp \left(-R_{\mathrm{1},\mathrm{app}}t\right)={\cos }^n\left(\mathit{\theta }\right)\exp \left(-R_{\mathrm{1}}t\right),\text{(25)}& {\displaystyle }{R}_{\mathrm{1},\mathrm{app}}=R_{\mathrm{1}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{TR}}}\ln \cos \left(\mathit{\theta }\right).\end{array}$$

In the previous examples in Fig. 2c and d, with a typical T₁=60 s (Keshari and Wilson, 2014) corresponding to $R_{\mathrm{1}}=\mathrm{1.67}\times {\mathrm{10}}^{-\mathrm{2}}$ s⁻¹ and a TR=4.25 s, the flip angle correction for a β=1^∘ pulse was $\mathrm{3.58}\times {\mathrm{10}}^{-\mathrm{5}}$, which “for all intents and purposes” is negligible, giving $R_{\mathrm{1},\mathrm{app}}=\mathrm{1.67}\times {\mathrm{10}}^{-\mathrm{2}}$ s⁻¹ and $T_{\mathrm{1},\mathrm{app}}=\mathrm{59.87}$ s. Hence, the time dependence of the signal shown in Fig. 2c is a robust reflection of the M_z(t) seen in Fig. 2a. For β=20^∘ the flip angle correction was $\mathrm{1.46}\times {\mathrm{10}}^{-\mathrm{2}}$, giving $R_{\mathrm{1},\mathrm{app}}=\mathrm{3.13}\times {\mathrm{10}}^{-\mathrm{2}}$ s⁻¹ and $T_{\mathrm{1},\mathrm{app}}=\mathrm{31.95}$ s. Therefore, for the larger flip angle there was a tradeoff between the increased sensitivity and the corresponding reduction in T_1,app with the more rapid decay of the NMR signal. The time dependence seen in Fig. 2d is no longer a good reflection of the M_z(t) shown in Fig. 2a. We conclude that when the RF flip angle is small, <1^∘, and the magnetization is sampled many times, the flip angle correction is negligible; accordingly, it is ignored in the next sections.

3 Relaxation of hyperpolarized magnetization in ¹³C substrates

We now take a detour into relaxation theory to give an overview of the factors that determine the values of $R_{\mathrm{1}}=\mathrm{1}/T_{\mathrm{1}}$ of hyperpolarized ¹³C solutes in a (bio)chemical system taking into account the main relaxation mechanisms responsible for the decay of the nuclear magnetization in the solution state at temperatures between ∼20 and 180 ^∘C and static magnetic field strengths between 1 mT and 23.5 T. The spin interactions discussed here are relevant to the outcome of numerous dissolution-dynamic nuclear polarization (dDNP) experiments.

A master equation for spin systems far from equilibrium based on a Lindblad dissipator formalism has recently been presented and shown to correctly predict the spin dynamics of hyperpolarized systems (Bengs and Levitt, 2020). In brief, Eq. (2) is only valid for the high temperature limit and weak-order approximation of a spin system at thermal equilibrium. However, we do not pursue this line of enquiry here because for the enzyme systems studied thus far with dDNP a constant value of T₁ has been statistically satisfactory in regression analyses of the data (Pagès et al., 2013; Shishmarev et al., 2018b).

Once a sufficiently high level of nuclear spin polarization has been achieved by implementing dDNP methodologies (often for ¹³C nuclei P_C>60 %), a jet of superheated solvent (e.g. H₂O and/or D₂O at 150–180 ^∘C) is injected directly onto the hyperpolarized sample (Ardenkjaer-Larsen et al., 2003; Wolber et al., 2004). Upon contact with the warm solvent, the frozen sample rapidly dissolves and is then transferred under the pressure of helium gas (6–9 bar) to a separate NMR/MRI spectrometer for the detection of hyperpolarized MRS signals or to a collection/quality control point for use in biological applications (Comment and Merritt, 2014). There are several potential issues related to spin relaxation during these processes, and we focus on nuclear spin relaxation in solution during the sample transfer stage (i.e. subject to changes in magnetic field strength) or situations where a solute has an altered rotational correlation time (i.e. dependence on temperature or when bound to a protein). This requires an understanding of the (potentially) large variety of molecular interactions that give rise to nuclear spin relaxation.

• Dipole–dipole couplings (DD). The dominant mechanism for the relaxation of nuclear spin magnetization is often the stochastic modulation of dipole–dipole interactions (couplings) to other nuclei, either in the same molecule or other molecules, including the solvent, as the molecule re-orientates in solution by tumbling.

• Chemical shift anisotropy (CSA). Nuclear spins resonate at different frequencies depending on the chemical shielding imposed by the local electronic environment and its orientation (a tensor property). The modulation of the chemical shift tensor by molecular tumbling in solution has a quadratic dependence on the strength of the static magnetic field and therefore increases markedly with B₀ (Kowalewski and Maler, 2019).

• Paramagnetic sites. Dissolved paramagnetic solutes (often impurities, but they can be purposely added as required by the experimental design), such as radical agents that remain in the dissolution solvent, molecular oxygen, and metal ions, which can be deleterious to the nuclear-spin relaxation, particularly in regions of low magnetic field (Blumberg, 1960; Pell et al., 2019). However, all species can be easily scavenged by co-dissolving chelating agents in the dissolution medium (Mieville et al., 2010).

• Scalar relaxation of the second kind. This mechanism operates when the nuclei of interest have scalar couplings to neighbouring nuclei that also relax rapidly (Pileio, 2011; Kubica et al., 2014; Elliott et al., 2019). In dDNP NMR experiments this relaxation mechanism is often enhanced during sample transfer steps through areas of low magnetic field (Chiavazza et al., 2013; Kubica et al., 2014).

• Spin rotation. The coupling of nuclear magnetization to that of a whole molecule or to mobile parts of a molecule, e.g. methyl groups, can act as an efficient relaxation mechanism. This mechanism has an unusual dependence on temperature, with the relaxation rate usually increasing at higher temperatures (Matson, 1977).

• Quadrupolar. Many molecules of interest in dDNP experiments contain either ²H or ¹⁴N nuclei. NMR relaxation times of such nuclei are often <1 s and therefore not sufficiently long to be relevant for dDNP experiments. However, there are two notable exceptions in ⁶Li⁺ and ¹³³Cs⁺, which have small nuclear quadrupole moments and therefore have intrinsically long T₁ values (van Heeswijk et al., 2009; Kuchel et al., 2019).

Derivations of relaxation rate expressions are well established and based on plausible physical models. For simplicity, we skip the majority of these since they are comprehensively presented by several authors (Kowalewski and Maler, 2019), and instead we focus on the main results of their analyses. Assuming a two-spin system composed of a ¹³C and ¹H, equations for the ¹³C–¹H dipole–dipole and the ¹³C CSA contributions to the ¹³C longitudinal relaxation rate constant (R₁) are given by Keeler (Keeler, 2010):

$$\begin{array}{}\text{(26)}& {\displaystyle }\begin{array}{rl}{R}_{\mathrm{1},\mathrm{DD}}& =b_{\mathrm{HC}}^{\mathrm{2}}\left[{\displaystyle \frac{\mathrm{3}}{\mathrm{20}}}J\left({\mathit{\omega }}_{\mathrm{C}}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{20}}}J\left({\mathit{\omega }}_{\mathrm{H}}-{\mathit{\omega }}_{\mathrm{C}}\right)\right.\\ & \left.+{\displaystyle \frac{\mathrm{3}}{\mathrm{10}}}J\left({\mathit{\omega }}_{\mathrm{H}}+{\mathit{\omega }}_{\mathrm{C}}\right)\right],\end{array}\text{(27)}& {\displaystyle }{R}_{\mathrm{1},\mathrm{CSA}}=c^{\mathrm{2}}\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{15}}}J\left({\mathit{\omega }}_{\mathrm{C}}\right)\right],\end{array}$$

where b_HC is the dipole–dipole coupling constant, defined as

$$\begin{array}{}\text{(28)}& b_{\mathrm{HC}}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}{\mathit{\gamma }}_{\mathrm{H}}{\mathit{\gamma }}_{\mathrm{C}}\mathrm{\hslash }}{\mathrm{4}\mathit{\pi }{r}_{\mathrm{HC}}^{\mathrm{3}}}},\end{array}$$

and c is the magnitude of the CSA assuming an axially symmetrical tensor given by

$$\begin{array}{}\text{(29)}& c={\mathit{\gamma }}_{\mathrm{C}}{B}_{\mathrm{0}}\left({\mathit{\sigma }}_{\parallel }-{\mathit{\sigma }}_{\perp }\right),\end{array}$$

where γ_H and γ_C are the magnetogyric ratios of the ¹H and ¹³C spins, respectively, r_HC is the internuclear distance between the ¹H and ¹³C atoms and σ_∥ and σ_⊥ are the parallel and perpendicular components of the axially symmetrical CSA tensor, respectively.

The so-called spectral density function that is a function of the Larmor frequency, ω, is

$$\begin{array}{}\text{(30)}& J\left(\mathit{\omega }\right)={\displaystyle \frac{\mathrm{2}{\mathit{\tau }}_{\mathrm{c}}}{\mathrm{1}+{\mathit{\omega }}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{c}}^{\mathrm{2}}}},\end{array}$$

where τ_c is the rotational correlation time (tumbling motion) of the re-orientating spin-bearing molecule in solution. The overall longitudinal relaxation rate constant is the sum of these two dominant contributions and is given by

$$\begin{array}{}\text{(31)}& R_{\mathrm{1}}=R_{\mathrm{1},\mathrm{DD}}+R_{\mathrm{1},\mathrm{CSA}}.\end{array}$$

3.1 Relaxation analysis

It is important (for experimental design purposes) to note the influence that a nearby ¹H spin has on the ¹³C nuclear T₁. Figure 3a shows the calculated ¹³C T₁ for a fixed rotational correlation time of ${\mathit{\tau }}_{\mathrm{c}}=\mathrm{0.4}\times {\mathrm{10}}^{-\mathrm{11}}$ s (previously reported for glycine in saline at 310 K – Endre et al., 1983), ¹³C CSA ${\mathit{\sigma }}_{\parallel }-{\mathit{\sigma }}_{\perp }=-\mathrm{98}$ ppm (previously reported for phosphoenolpyruvate – Bechmann et al., 2004) and a magnetic field strength of B₀=7 T as a function of the ¹H–¹³C internuclear distance r_HC. Biaxiality of the CSA interaction has been ignored here. A rapid rise occurs in T₁ as the ¹H–¹³C internuclear separation increases. In the case of r_HC=1.09 Å, which is typical of a ¹H–¹³C single bond, the ¹³C nuclear T₁ is predicted to be ∼11.4 s. The ¹H–¹³C dipole–dipole coupling constant scales with $r_{\mathrm{HC}}^{-\mathrm{3}}$; consequently, the presence of a directly bonded proton significantly shortens the relaxation time constant of the ¹³C magnetization. Small molecules containing ¹³C atoms that do not have directly bonded ¹H, or at least ¹H spins located at significant internuclear distances, are required. Such moieties include the carboxyl group that is present in many low molecular weight metabolites such as pyruvate, lactate, and methylglyoxal (Shishmarev et al., 2018a). At the longer ¹H–¹³C internuclear distance of 1.45 Å, implying a ¹H–¹³C dipole–dipole coupling constant of $b_{\mathrm{HC}}/\mathrm{2}\mathit{\pi }=-\mathrm{10.2}$ kHz, a ¹³C nuclear T₁ of ∼60 s is predicted. At very long distances, the ¹³C relaxation time constant will tend to that of the CSA relaxation contribution alone.

Figure 3 (a) Simulation of the ¹³C nuclear T₁ for a two-spin ¹H–¹³C system as a function of the internuclear distance (r_HC) with a rotational correlation time ${\mathit{\tau }}_{\mathrm{c}}=\mathrm{0.4}\times {\mathrm{10}}^{-\mathrm{11}}$ s, ¹³C CSA ${\mathit{\sigma }}_{\parallel }-{\mathit{\sigma }}_{\perp }=-\mathrm{98}$ ppm and at a magnetic field strength B=7 T. (b) Dependence of the ¹³C nuclear T₁ as a function of the magnetic field B and the rotational correlation time τ_c.

The dependence of R₁ on temperature and molecular size (e.g. due to binding) scales with the rotational correlation time. Figure 3b shows the dependence of the ¹³C nuclear T₁ (1/R₁) as a function of τ_c and B₀ for this 2-spin-1/2 system with r_HC=1.45 Å and ${\mathit{\sigma }}_{\parallel }-{\mathit{\sigma }}_{\perp }=-\mathrm{98}$ ppm. In the extreme narrowing limit, i.e. ${\mathit{\omega }}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{c}}^{\mathrm{2}}\ll \mathrm{1}$, the following familiar equations describe the relaxation of ¹³C spins under the dipole–dipole and CSA relaxation mechanisms (Kowalewski and Maler, 2019):

$$\begin{array}{}\text{(32)}& {\displaystyle }{R}_{\mathrm{1},\mathrm{DD}}=b_{\mathrm{HC}}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{c}},\text{(33)}& {\displaystyle }{R}_{\mathrm{1},\mathrm{CSA}}={\displaystyle \frac{\mathrm{2}}{\mathrm{15}}}{c}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{c}}.\end{array}$$

In the extreme narrowing regime the ¹³C nuclear T₁ becomes shorter with increasing magnetic field strength due to the $B_{\mathrm{0}}^{\mathrm{2}}$ dependence of R_1,CSA. At low field strengths, the magnitude of T₁ will mostly be attributed to dipole–dipole relaxation with the nearby ¹H spin. It is also worth noting that the ¹³C T₁ follows the usual Lorentzian spectral density functional dependence on the rotational correlation time. This is clearly seen at high magnetic fields.

3.2 Molecular considerations

The majority of dDNP experiments used to study biological systems employ ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}/{\mathrm{D}}_{\mathrm{2}}\mathrm{O}$ as the dissolution solvent. Detection of hyperpolarized NMR/MRI signals typically occurs in a magnetic field range of 1.5–9.4 T; thus, Fig. 3b indicates a ¹³C nuclear T₁ of the order of ∼60 s for a carbonyl group, and this is commonly seen in practice (Shishmarev et al., 2018a). It is important to remember that Eqs. (26)–(31) provide a greatly simplified picture of the problem at hand; in reality there are many magnetic nuclei (often within the same molecule) which contribute to the relaxation of ¹³C magnetization. The additional dipole–dipole interactions are likely to be responsible for differences between predicted and measured ¹³C relaxation times, along with the other (more exotic) signal attenuation mechanisms that are described above.

In a dDNP experiment the dissolution and transfer process can take as long as 15 s; it depends on the distance to the point of use from the polarizing source, and in clinical applications an additional 30 s can easily be added for quality control processes. Such requirements place a bound on the usable time in which hyperpolarized ¹³C magnetization must be maintained, and it is typical to expect 45 s to be this limit. Given that the magnetic field strength “felt” by the hyperpolarized sample can be controlled (to a reasonable extent) throughout its voyage between the dDNP polarizer and the point of use (Milani et al., 2015), the rotational correlation time becomes the most important factor that impacts upon the ¹³C nuclear T₁. Figure 3b indicates that even for a rotational correlation time of the order of ${\mathit{\tau }}_{\mathrm{c}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{10}}$ s, such as found in proteins in solution (Wilbur et al., 1976), Eqs. (26)–(31) yield ¹³C nuclear T₁ relaxation times which are too short to allow practical use of such samples, i.e. $\mathrm{5}\times T_{\mathrm{1}}\ll \mathrm{45}$ s, in comparison to the overall time required by a dDNP experiment.

A major parameter that controls the magnitude of the rotational correlation time of a spin-bearing molecule is its molecular weight (M_w). Since τ_c∝M_w the rotational correlation time has a noticeable impact on the ¹³C nuclear T₁ with even the smallest increase in molecular weight. In order to achieve ¹³C nuclear T₁ relaxation times that are sufficiently long to enable hyperpolarized ¹³C magnetization to survive the dissolution and transfer process, the ¹³C NMR signals must be detectable above the spectral noise for ∼45 s. Hence, dDNP samples used in biological experiments are currently restricted to small molecules (or ions – van Heeswijk et al., 2009; Kuchel et al., 2019). For example, the estimate of ∼60 s for the ¹³C nuclear T₁ of the model system described above was predicted with a rotational correlation time of ${\mathit{\tau }}_{\mathrm{c}}=\mathrm{0.4}\times {\mathrm{10}}^{-\mathrm{11}}$ (Endre et al., 1983), and this is sufficiently long for dDNP experiments.

3.3 Enzyme binding

The worst-case scenario for the model system described in Fig. 3b would be a moderate rotational correlation time of the order of ${\mathit{\tau }}_{\mathrm{c}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{8}}-\mathrm{1}\times {\mathrm{10}}^{-\mathrm{10}}$ s for which ¹³C nuclear T₁ relaxation times in the millisecond regime are predicted. Such correlation times correspond to a system with a molecular weight comparable to that of an enzyme. If the small molecule (ligand) or ion becomes bound to the enzyme, then it will assume the rotational correlation time of the higher mass binding partner. In the case of ${\mathit{\tau }}_{\mathrm{c}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{9}}$ for an enzyme–ligand complex, a ¹³C substrate will have a predicted nuclear magnetic T₁ of ∼276.4 ms at a static magnetic field strength of 7 T. Such a stark variation in ¹³C nuclear T₁ values provides good contrast in relaxation-based ligand–protein binding experiments (Valensin et al., 1982).

4 Mechanistic description of reaction kinetics of hyperpolarized substrates

We now consider the interpretation of hyperpolarized dynamics for complex chemical reactions. To help tease apart the key features of the analysis, we begin with some simplifying assumptions. First, in the absence of an RF pulse Eq. (20) becomes block diagonal, since transverse and longitudinal magnetization are not interconverted. The evolution of the z magnetization is then dependent only on the initial conditions, T₁, and the rate constants that characterize the chemical exchange. Second, we assume that the z magnetization is sampled many times with an infinitesimally small flip angle (≪1^∘) so the longitudinal magnetization decays with its intrinsic T₁ value rather than an apparent T_1,app value. Finally, the hyperpolarized magnetization decays to zero; i.e. the enhancement factor η (Eq. 21) is such that M₀ is greater than M_eq by many orders of magnitude. Thus, the equilibrium magnetization at t=∞ is effectively zero, and it can be ignored in the analysis of real experimental data.

To reduce clutter in the equations, for all the discussions that now follow, we drop the subscript z since we hereafter deal only with longitudinal magnetization and denote $M_{z,\mathrm{hyp}}^{\mathrm{A}}$ and $M_{z,\mathrm{hyp}}^{\mathrm{B}}$ as A^∗(t) and B^∗(t) corresponding to hyperpolarized magnetization (identified with an asterisk ^∗).

4.1 Simple first-order exchange kinetics of hyperpolarized substrates

Confining our analysis to the physical subspace that is composed of longitudinal magnetizations, which describe first-order kinetics of a two-site exchange reaction of hyperpolarized substrates, ${\mathrm{A}}^{\ast }\leftrightarrow {\mathrm{B}}^{\ast }$, Eq. (20) simplifies to

$$\begin{array}{}\text{(34)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\end{array}\right]=\left[\begin{array}{cc}-k_{\mathrm{1}}-R_{\mathrm{1}}^{\mathrm{A}}& k_{-\mathrm{1}}\\ k_{\mathrm{1}}& -k_{-\mathrm{1}}-R_{\mathrm{1}}^{\mathrm{B}}\end{array}\right]\left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\end{array}\right].\end{array}$$

Equivalently, Eq. (34) can be expanded to give

$$\begin{array}{}\text{(35)}& {\displaystyle \frac{{\mathrm{dA}}^{\ast }\left(t\right)}{\mathrm{d}t}}=-k_{\mathrm{1}}{\mathrm{A}}^{\ast }\left(t\right)+k_{-\mathrm{1}}{\mathrm{B}}^{\ast }\left(t\right)-R_{\mathrm{1}}^{\mathrm{A}}{\mathrm{A}}^{\ast }\left(t\right),\text{(36)}& {\displaystyle \frac{{\mathrm{dB}}^{\ast }\left(t\right)}{\mathrm{d}t}}=k_{\mathrm{1}}{\mathrm{A}}^{\ast }\left(t\right)-k_{-\mathrm{1}}{\mathrm{B}}^{\ast }\left(t\right)-R_{\mathrm{1}}^{\mathrm{B}}{\mathrm{B}}^{\ast }\left(t\right),\end{array}$$

where k₁ and k₋₁ denote first-order rate constants, and $R_{\mathrm{1}}^{\mathrm{A}}=\mathrm{1}/T_{\mathrm{1}}^{\mathrm{A}}$ and $R_{\mathrm{1}}^{\mathrm{B}}=\mathrm{1}/T_{\mathrm{1}}^{\mathrm{B}}$ are the longitudinal relaxation rate constants of A and B, respectively.

Since Eqs. (35) and (36) describe the time evolution of the z magnetizations (that is proportional to concentration/mass), they do not satisfy the conservation of mass requirement because $\mathrm{d}\left[{\mathrm{A}}^{\ast }\left(t\right)+{\mathrm{B}}^{\ast }\left(t\right)\right]/\mathrm{d}t=-R_{\mathrm{1}}^{\mathrm{A}}{\mathrm{A}}^{\ast }\left(t\right)-R_{\mathrm{1}}^{\mathrm{B}}{\mathrm{B}}^{\ast }\left(t\right)$, and this tends to zero with time. However, the equations can be recast to specify that the pools of hyperpolarized substrates relax to form pools of non-polarized substrates A↔B. These pools are denoted simply by A(t) and B(t) (without the asterisks) as shown in Fig. 4a. The analogy with radioactive tracers is a useful one here. A “hot” pool of radioactive material decays with first-order kinetics (half-life) to form a “cold” pool of non-radioactive material with the sum of “hot” and “cold” being constant.

Figure 4 Simulated first-order two-site exchange kinetics of hyperpolarized solutes, ${\mathrm{A}}^{\ast }\leftrightarrow {\mathrm{B}}^{\ast }$, conforming to conservation of mass, assuming initial hyperpolarized magnetization of only solute ${\mathrm{A}}^{\ast }\left(\mathrm{0}\right)=\mathrm{1}$. Longitudinal relaxation rate constants were $R_{\mathrm{1}}^{\mathrm{A}}=R_{\mathrm{1}}^{\mathrm{B}}=\mathrm{1}/\mathrm{60}$ s⁻¹. The time dependences of A^∗(t), A(t), B^∗(t) and B(t) (left panel) were calculated numerically using Eq. (35)–(38) with rate constants (b) k₁=0.01 s⁻¹, $k_{-\mathrm{1}}=\mathrm{0}$ s⁻¹, corresponding to uni-directional kinetics, (c) k₁=0.01 s⁻¹, $k_{-\mathrm{1}}=\mathrm{0.005}$ s⁻¹ and (d) k₁=0.01 s⁻¹, $k_{-\mathrm{1}}=\mathrm{0.01}$ s${}^{-\mathrm{1},}$ corresponding to exchange kinetics. The right panels show plots of the time dependence of ${\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)=\left[\mathrm{A}\left(t\right)\right]$ and ${\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)=\left[\mathrm{B}\left(t\right)\right]$.

The kinetics of the non-polarized pools are described by

$$\begin{array}{}\text{(37)}& {\displaystyle \frac{\mathrm{dA}\left(t\right)}{\mathrm{d}t}}=-k_{\mathrm{1}}\mathrm{A}\left(t\right)+k_{-\mathrm{1}}\mathrm{B}\left(t\right)+R_{\mathrm{1}}^{\mathrm{A}}{\mathrm{A}}^{\ast }\left(t\right),\text{(38)}& {\displaystyle \frac{\mathrm{dB}\left(t\right)}{\mathrm{d}t}}=k_{\mathrm{1}}\mathrm{A}\left(t\right)-k_{-\mathrm{1}}\mathrm{B}\left(t\right)+R_{\mathrm{1}}^{\mathrm{B}}{\mathrm{B}}^{\ast }\left(t\right).\end{array}$$

Equations (35)–(38) now satisfy conservation of mass, since the rate of change $\mathrm{d}\left[{\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)+{\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)\right]/\mathrm{d}t$ is always zero. Note that A(t) and B(t) are not observed in the dDNP NMR experiment, but they are the counterparts of real concentrations of solute that would be assayable (bio)chemically.

Equations (35)–(38) can be written as

$$\begin{array}{}\text{(39)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\\ \mathrm{A}\left(t\right)\\ \mathrm{B}\left(t\right)\end{array}\right]=\\ & \left[\begin{array}{cccc}-k_{\mathrm{1}}-R_{\mathrm{1}}^{\mathrm{A}}& k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ k_{\mathrm{1}}& -k_{-\mathrm{1}}-R_{\mathrm{1}}^{\mathrm{B}}& \mathrm{0}& \mathrm{0}\\ R_{\mathrm{1}}^{\mathrm{A}}& \mathrm{0}& -k_{\mathrm{1}}& k_{-\mathrm{1}}\\ \mathrm{0}& R_{\mathrm{1}}^{\mathrm{B}}& k_{\mathrm{1}}& -k_{-\mathrm{1}}\end{array}\right]\left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\\ \mathrm{A}\left(t\right)\\ \mathrm{B}\left(t\right)\end{array}\right].\end{array}\end{array}$$

Equation (39) can be written as

$$\begin{array}{}\text{(40)}& {\displaystyle \frac{\mathrm{d}\mathit{M}\left(t\right)}{\mathrm{d}t}}=L\mathit{M}\left(t\right).\end{array}$$

We can apply a similarity transform given by

$$\begin{array}{}\text{(41)}& U=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{1}\end{array}\right],\end{array}$$

to yield an equation of motion in a transformed vector basis

$$\begin{array}{}\text{(42)}& {\displaystyle \frac{\mathrm{d}{\mathit{M}}^{\prime }\left(t\right)}{\mathrm{d}t}}=UL{U}^{-\mathrm{1}}{\mathit{M}}^{\prime }\left(t\right),\end{array}$$

given by

$$\begin{array}{}\text{(43)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\\ {\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)\end{array}\right]=\\ & \left[\begin{array}{cccc}-k_{\mathrm{1}}-R_{\mathrm{1}}^{\mathrm{A}}& k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ k_{\mathrm{1}}& -k_{-\mathrm{1}}-R_{\mathrm{1}}^{\mathrm{B}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}& k_{-\mathrm{1}}\\ \mathrm{0}& \mathrm{0}& k_{\mathrm{1}}& -k_{-\mathrm{1}}\end{array}\right]\left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\\ {\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)\end{array}\right].\end{array}\end{array}$$

We can now appreciate the equivalence between this formalism and conventional chemical reaction kinetics written in terms of molecular concentrations. For first-order reactions, the hyperpolarized magnetization evolves according to the Bloch–McConnell equations, while the concentrations given by the sum of the “hot” and “cold” pools evolve according to the conventional form of chemical reaction kinetics for a closed system. Therefore, ${\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)$ and ${\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)$ are proportional to [A(t)] and [B(t)], respectively, where the constant of proportionality depends on the initial experimental conditions viz. [A]₀ and [B]₀. In other words, provided ${\mathrm{A}}^{\ast }\left(\mathrm{0}\right)+\mathrm{A}\left(\mathrm{0}\right)={\left[\mathrm{A}\right]}_{\mathrm{0}}$ and ${\mathrm{B}}^{\ast }\left(\mathrm{0}\right)+\mathrm{B}\left(\mathrm{0}\right)={\left[\mathrm{B}\right]}_{\mathrm{0}}$, then the constant of proportionality is 1, and we can equate ${\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)=\left[\mathrm{A}\left(t\right)\right]$ and ${\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)=\left[\mathrm{B}\left(t\right)\right]$. This is a crucial concept that we return to below.

Figure 4 shows numerical simulations of the time evolution of the system described by Eq. (39) with an initial magnetization vector $\mathit{M}\left(\mathrm{0}\right)=\left[\mathrm{1},\mathrm{0},\mathrm{0},\mathrm{0}\right]$ that corresponds to only hyperpolarized ${\mathrm{A}}^{\ast }\left(\mathrm{0}\right)=\mathrm{1}$ and longitudinal relaxation rate constants $R_{\mathrm{1}}^{\mathrm{A}}=R_{\mathrm{1}}^{\mathrm{B}}=\mathrm{1}/\mathrm{60}$ s⁻¹. The time dependence of A^∗(t), A(t), B^∗(t) and B(t) were calculated numerically (left panel) for different rate constants: Fig. 4b, k₁=0.01 s⁻¹, $k_{-\mathrm{1}}=\mathrm{0}$ s⁻¹, corresponding to a uni-directional reaction; Fig. 4c, k₁=0.01 s⁻¹, $k_{-\mathrm{1}}=\mathrm{0.005}$ s⁻¹, corresponding to bi-directional exchange with an equilibrium constant K=2; and Fig. 4d, k₁=0.01 s⁻¹, $k_{-\mathrm{1}}=\mathrm{0.01}$ s⁻¹, also corresponding to bi-directional exchange with an equilibrium constant K=1. The right columns of plots show the time dependence of ${\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)$ and ${\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)$ that reproduce conventional kinetics of [A(t)] and [B(t)], as required for mathematical and physical consistency.

The approach used here (as laid out in Kuchel and Shishmarev, 2020) enables us to create systems of differential equations that satisfy conservation of mass and therefore allow a study of the influence of non-hyperpolarized pools of substrates on reaction kinetics. The approach enables more complicated reaction mechanisms to be described to allow the inclusion of MR-invisible pools of substrates such as ¹²C, which are known to affect the outcome of dDNP experiments in vivo. We consider some of these scenarios next.

4.2 Sequential reaction kinetics of hyperpolarized substrates

Equation (39) can be extended to compartmental models of arbitrary complexity: consider a reaction scheme involving three substrates ${\mathrm{A}}^{\ast }\leftrightarrow {\mathrm{B}}^{\ast }\leftrightarrow {\mathrm{C}}^{\ast }$ which relax through T₁ processes to form a pool of non-polarized substrates $\mathrm{A}\leftrightarrow \mathrm{B}\leftrightarrow \mathrm{C}$, as shown in Fig. 5a. This is analogous to a system where a solution of hyperpolarized solute A^∗ is introduced into the extracellular medium in a cell suspension, is transported into the cells where it is denoted by B^∗ and is subsequently acted upon by an enzyme to form C^∗. The system of differential equations that describe the kinetics of this scheme is

$$\begin{array}{}\text{(44)}& {\displaystyle \frac{{\mathrm{dA}}^{\ast }\left(t\right)}{\mathrm{d}t}}=-k_{\mathrm{1}}{\mathrm{A}}^{\ast }\left(t\right)+k_{-\mathrm{1}}{\mathrm{B}}^{\ast }\left(t\right)-R_{\mathrm{1}}^{\mathrm{A}}{\mathrm{A}}^{\ast }\left(t\right),\text{(45)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{{\mathrm{dB}}^{\ast }\left(t\right)}{\mathrm{d}t}}& =k_{\mathrm{1}}{\mathrm{A}}^{\ast }\left(t\right)-k_{-\mathrm{1}}{\mathrm{B}}^{\ast }\left(t\right)-k_{\mathrm{2}}{\mathrm{B}}^{\ast }\left(t\right)\\ & +k_{-\mathrm{2}}{\mathrm{C}}^{\ast }\left(t\right)-R_{\mathrm{1}}^{\mathrm{B}}{\mathrm{B}}^{\ast }\left(t\right),\end{array}\text{(46)}& {\displaystyle \frac{{\mathrm{dC}}^{\ast }\left(t\right)}{\mathrm{d}t}}=k_{\mathrm{2}}{\mathrm{B}}^{\ast }\left(t\right)-k_{-\mathrm{2}}{\mathrm{C}}^{\ast }\left(t\right)-R_{\mathrm{1}}^{\mathrm{C}}{\mathrm{C}}^{\ast }\left(t\right),\text{(47)}& {\displaystyle \frac{\mathrm{dA}\left(t\right)}{\mathrm{d}t}}=-k_{\mathrm{1}}\mathrm{A}\left(t\right)+k_{-\mathrm{1}}\mathrm{B}\left(t\right)+R_{\mathrm{1}}^{\mathrm{A}}{\mathrm{A}}^{\ast }\left(t\right),\text{(48)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dB}\left(t\right)}{\mathrm{d}t}}& =k_{\mathrm{1}}\mathrm{A}\left(t\right)-k_{-\mathrm{1}}\mathrm{B}\left(t\right)-k_{\mathrm{2}}\mathrm{B}\left(t\right)\\ & +k_{-\mathrm{2}}\mathrm{C}\left(t\right)+R_{\mathrm{1}}^{\mathrm{B}}{\mathrm{B}}^{\ast }\left(t\right),\end{array}\text{(49)}& {\displaystyle \frac{\mathrm{dC}\left(t\right)}{\mathrm{d}t}}=k_{\mathrm{2}}\mathrm{B}\left(t\right)-k_{-\mathrm{2}}\mathrm{C}\left(t\right)+R_{\mathrm{1}}^{\mathrm{C}}{\mathrm{C}}^{\ast }\left(t\right),\end{array}$$

where we have removed the square brackets that denote molar concentration to avoid some of the clutter. However, it is important to recall that there is a factor that relates magnetization to concentration, and this is estimated from the known initial experimental conditions.

Figure 5 Simulated first-order three-site exchange kinetics of hyperpolarized solutes, ${\mathrm{A}}^{\ast }\leftrightarrow {\mathrm{B}}^{\ast }\leftrightarrow {\mathrm{C}}^{\ast }$, conforming to conservation of mass, assuming initial hyperpolarized magnetization of only solute ${\mathrm{A}}^{\ast }\left(\mathrm{0}\right)=\mathrm{1}$. Longitudinal relaxation rate constants were $R_{\mathrm{1}}^{\mathrm{A}}=R_{\mathrm{1}}^{\mathrm{B}}=R_{\mathrm{1}}^{\mathrm{C}}=\mathrm{1}/\mathrm{60}$ s⁻¹. The time dependencies of A^∗(t), A(t), B^∗(t), B(t), C^∗(t) and C(t) (left column of plots) were calculated numerically using Eqs. (41)–(46) with rate constants (b) $k_{\mathrm{1}}=k_{\mathrm{2}}=\mathrm{0.01}$ s⁻¹, $k_{-\mathrm{1}}=k_{-\mathrm{2}}=\mathrm{0}$ s⁻¹, corresponding to uni-directional kinetics, (c) $k_{\mathrm{1}}=k_{\mathrm{2}}=\mathrm{0.01}$ s⁻¹, $k_{-\mathrm{1}}=k_{-\mathrm{2}}=\mathrm{0.005}$ s⁻¹ and (d) $k_{\mathrm{1}}=k_{\mathrm{2}}=k_{-\mathrm{1}}=k_{-\mathrm{2}}=\mathrm{0.01}$ s⁻¹, corresponding to exchange kinetics. The right column of plots shows the time dependence of ${\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)=\left[\mathrm{A}\left(t\right)\right]$, ${\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)=\left[\mathrm{B}\left(t\right)\right]$ and ${\mathrm{C}}^{\ast }\left(t\right)+\mathrm{C}\left(t\right)=\left[\mathrm{C}\left(t\right)\right]$.

Equations (44)–(49) can be recast in matrix form to give

$$\begin{array}{}\text{(50)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\\ {\mathrm{C}}^{\ast }\left(t\right)\\ \mathrm{A}\left(t\right)\\ \mathrm{B}\left(t\right)\\ \mathrm{C}\left(t\right)\end{array}\right]=\\ & \left[\begin{array}{cccccc}-k_{\mathrm{1}}-R_{\mathrm{1}}^{\mathrm{A}}& k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ k_{\mathrm{1}}& -k_{-\mathrm{1}}-k_{\mathrm{2}}-R_{\mathrm{1}}^{\mathrm{B}}& k_{-\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& k_{\mathrm{2}}& -k_{-\mathrm{2}}-R_{\mathrm{1}}^{\mathrm{C}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ R_{\mathrm{1}}^{\mathrm{A}}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}& k_{-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& R_{\mathrm{1}}^{\mathrm{B}}& \mathrm{0}& k_{\mathrm{1}}& -k_{-\mathrm{1}}-k_{\mathrm{2}}& k_{-\mathrm{2}}\\ \mathrm{0}& \mathrm{0}& R_{\mathrm{1}}^{\mathrm{C}}& \mathrm{0}& k_{\mathrm{2}}& -k_{-\mathrm{2}}\end{array}\right]\\ & \left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\\ {\mathrm{C}}^{\ast }\left(t\right)\\ \mathrm{A}\left(t\right)\\ \mathrm{B}\left(t\right)\\ \mathrm{C}\left(t\right)\end{array}\right].\end{array}\end{array}$$

It is readily verified that Eq. (50) satisfies conservation of mass, since the rate of change $\mathrm{d}\left({\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)+{\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)+{\mathrm{C}}^{\ast }\left(t\right)+\mathrm{C}\left(t\right)\right)/\mathrm{d}t=\mathrm{0}$.

We can apply a similarity transform given by

$$\begin{array}{}\text{(51)}& U=\left[\begin{array}{cccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right].\end{array}$$

To yield an equation of motion in the transformed basis vector given by

$$\begin{array}{}\text{(52)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\\ {\mathrm{C}}^{\ast }\left(t\right)\\ {\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)\\ {\mathrm{C}}^{\ast }\left(t\right)+\mathrm{C}\left(t\right)\end{array}\right]=\\ & \left[\begin{array}{cccccc}-k_{\mathrm{1}}-R_{\mathrm{1}}^{\mathrm{A}}& k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ k_{\mathrm{1}}& -k_{-\mathrm{1}}-k_{\mathrm{2}}-R_{\mathrm{1}}^{\mathrm{B}}& k_{-\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& k_{\mathrm{2}}& -k_{-\mathrm{2}}-R_{\mathrm{1}}^{\mathrm{C}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}& k_{-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{1}}& -k_{-\mathrm{1}}-k_{\mathrm{2}}& k_{-\mathrm{2}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{2}}& -k_{-\mathrm{2}}\end{array}\right]\\ & \left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\\ {\mathrm{C}}^{\ast }\left(t\right)\\ {\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)\\ {\mathrm{C}}^{\ast }\left(t\right)+\mathrm{C}\left(t\right)\end{array}\right],\end{array}\end{array}$$

the hyperpolarized magnetization evolves according to the Bloch–McConnell equations, while the concentrations given by the sum of the “hot” and “cold” pools evolve according to the conventional form of chemical reaction kinetics for a closed system. Therefore, provided ${\mathrm{A}}^{\ast }\left(\mathrm{0}\right)+\mathrm{A}\left(\mathrm{0}\right)={\left[\mathrm{A}\right]}_{\mathrm{0}}$, ${\mathrm{B}}^{\ast }\left(\mathrm{0}\right)+\mathrm{B}\left(\mathrm{0}\right)={\left[\mathrm{B}\right]}_{\mathrm{0}}$ and ${\mathrm{C}}^{\ast }\left(\mathrm{0}\right)+\mathrm{C}\left(\mathrm{0}\right)={\left[\mathrm{C}\right]}_{\mathrm{0}}$, then ${\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)=\left[\mathrm{A}\left(t\right)\right]$, ${\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)=\left[\mathrm{B}\left(t\right)\right]$ and ${\mathrm{C}}^{\ast }\left(t\right)+\mathrm{C}\left(t\right)=\left[\mathrm{C}\left(t\right)\right]$, respectively.

Figure 5 shows the results of numerical integration of Eq. (50) with the initial magnetization vector $\mathit{M}\left(\mathrm{0}\right)=\left[\mathrm{1},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0}\right]$ that corresponds to having only hyperpolarized ${\mathrm{A}}^{\ast }\left(\mathrm{0}\right)=\mathrm{1}$ and longitudinal relaxation rate constants $R_{\mathrm{1}}^{\mathrm{A}}=R_{\mathrm{1}}^{\mathrm{B}}=R_{\mathrm{1}}^{\mathrm{C}}=\mathrm{1}/\mathrm{60}$ s⁻¹. The time dependence of A^∗(t), A(t), B^∗(t), B(t), C^∗(t) and C(t) were calculated (left panel) for different rate constants: Fig. 5b, $k_{\mathrm{1}}=k_{\mathrm{2}}=\mathrm{0.01}$ s⁻¹, $k_{-\mathrm{1}}=k_{-\mathrm{2}}=\mathrm{0}$ s⁻¹, corresponding to uni-directional kinetics; Fig. 5c, $k_{\mathrm{1}}=k_{\mathrm{2}}=\mathrm{0.01}$ s⁻¹, $k_{-\mathrm{1}}=k_{-\mathrm{2}}=\mathrm{0.005}$ s⁻¹, corresponding to bi-directional exchange kinetics; and Fig. 5d, $k_{\mathrm{1}}=k_{\mathrm{2}}=k_{-\mathrm{1}}=k_{-\mathrm{2}}=\mathrm{0.01}$ s⁻¹, also corresponding to bi-directional exchange kinetics. The right column shows plots of the time dependence of ${\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)$, ${\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)$ and ${\mathrm{C}}^{\ast }\left(t\right)+\mathrm{C}\left(t\right)$, which reproduce the conventional chemical kinetics of [A(t)], [B(t)] and [C(t)], as required for mathematical and physical consistency.

4.3 Second-order kinetics of hyperpolarized substrates

We now describe hyperpolarized substrates A^∗(t) and B^∗(t) reacting with non-hyperpolarized substrates [C(t)] and [D(t)]. The system of differential equations describes the second-order kinetics of ${\mathrm{A}}^{\ast }+\mathrm{C}\leftrightarrow {\mathrm{B}}^{\ast }+\mathrm{D}$ with only the hyperpolarized pools relaxing through T₁ processes to form a pool of non-polarized substrates $A+\mathrm{C}\leftrightarrow \mathrm{B}+\mathrm{D}$. The reactant concentrations [C(t)] and [D(t)] are common to both pools, as shown in Fig. 6a. The relevant system of differential equations (again omitting the square brackets that denote concentration) is

$$\begin{array}{}\text{(53)}& {\displaystyle \frac{{\mathrm{dA}}^{\ast }\left(t\right)}{\mathrm{d}t}}=-k_{\mathrm{1}}\mathrm{C}\left(t\right){\mathrm{A}}^{\ast }\left(t\right)+k_{-\mathrm{1}}\mathrm{D}\left(t\right){\mathrm{B}}^{\ast }\left(t\right)-R_{\mathrm{1}}^{\mathrm{A}}{\mathrm{A}}^{\ast }\left(t\right),\text{(54)}& {\displaystyle \frac{{\mathrm{dB}}^{\ast }\left(t\right)}{\mathrm{d}t}}=k_{\mathrm{1}}\mathrm{C}\left(t\right){\mathrm{A}}^{\ast }\left(t\right)-k_{-\mathrm{1}}\mathrm{D}\left(t\right){\mathrm{B}}^{\ast }\left(t\right)-R_{\mathrm{1}}^{\mathrm{B}}{\mathrm{B}}^{\ast }\left(t\right),\text{(55)}& {\displaystyle \frac{\mathrm{dA}\left(t\right)}{\mathrm{d}t}}=-k_{\mathrm{1}}\mathrm{C}\left(t\right)\mathrm{A}\left(t\right)+k_{-\mathrm{1}}\mathrm{D}\left(t\right)\mathrm{B}\left(t\right)+R_{\mathrm{1}}^{\mathrm{A}}{\mathrm{A}}^{\ast }\left(t\right),\text{(56)}& {\displaystyle \frac{\mathrm{dB}\left(t\right)}{\mathrm{d}t}}=k_{\mathrm{1}}\mathrm{C}\left(t\right)\mathrm{A}\left(t\right)-k_{-\mathrm{1}}\mathrm{D}\left(t\right)\mathrm{B}\left(t\right)+R_{\mathrm{1}}^{\mathrm{B}}{\mathrm{B}}^{\ast }\left(t\right),\text{(57)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dC}\left(t\right)}{\mathrm{d}t}}& =-k_{\mathrm{1}}\left({\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)\right)\mathrm{C}\left(t\right)\\ & +k_{-\mathrm{1}}\left({\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)\right)\mathrm{D}\left(t\right),\end{array}\text{(58)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dD}\left(t\right)}{\mathrm{d}t}}& =k_{\mathrm{1}}\left({\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)\right)\mathrm{C}\left(t\right)\\ & -k_{-\mathrm{1}}\left({\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)\right)\mathrm{D}\left(t\right).\end{array}\end{array}$$

Figure 6 Simulated second-order exchange kinetics of hyperpolarized solutes, ${\mathrm{A}}^{\ast }+\mathrm{C}\leftrightarrow {\mathrm{B}}^{\ast }+\mathrm{D}$, conforming to conservation of mass, assuming initial hyperpolarized magnetization of only solute ${\mathrm{A}}^{\ast }\left(\mathrm{0}\right)=\mathrm{1}$. Longitudinal relaxation rate constants were $R_{\mathrm{1}}^{\mathrm{A}}=R_{\mathrm{1}}^{\mathrm{B}}=\mathrm{1}/\mathrm{60}$ s⁻¹. The time dependences of A^∗(t), A(t), B^∗(t) and B(t) were simulated (left column of plots) using Eqs. (53)–(58) with rate constants (b) k₁=0.01 M⁻¹ s⁻¹, $k_{-\mathrm{1}}=\mathrm{0}$ M⁻¹ s⁻¹, corresponding to uni-directional kinetics (c) k₁=0.01 M⁻¹ s⁻¹, $k_{-\mathrm{1}}=\mathrm{0.005}$ M⁻¹ s⁻¹ and (d) $k_{\mathrm{1}}=k_{-\mathrm{1}}=\mathrm{0.01}$ M⁻¹ s⁻¹, corresponding to exchange kinetics. The right column of plots shows the time dependence of ${\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)=\left[\mathrm{A}\left(t\right)\right]$, ${\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)=\left[\mathrm{B}\left(t\right)\right]$ and non-polarized reactants [C(t)] and [D(t)].

Again, mass is conserved, as seen by the fact that $\mathrm{d}\left(\right({\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)+{\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right))/\mathrm{d}t=\mathrm{0}$ and $\mathrm{d}\left(\mathrm{C}\left(t\right)+\mathrm{D}\left(t\right)\right)/\mathrm{d}t=\mathrm{0}$. Also, recall that, provided ${\mathrm{A}}^{\ast }\left(\mathrm{0}\right)+\mathrm{A}\left(\mathrm{0}\right)={\left[\mathrm{A}\right]}_{\mathrm{0}}$, ${\mathrm{B}}^{\ast }\left(\mathrm{0}\right)+\mathrm{B}\left(\mathrm{0}\right)={\left[\mathrm{B}\right]}_{\mathrm{0}}$, C(0)=[C]₀ and D(0)=[D]₀, then we can make use of the equalities ${\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)=\left[\mathrm{A}\left(t\right)\right]$, ${\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)=\left[\mathrm{B}\left(t\right)\right]$, C(t)=[C(t)] and D(t)=[D(t)], respectively. It is now very evident why we must equate the initial signal with the concentration via an experimentally estimated scaling factor.

Equations (53)–(58) can be written in matrix vector form as

$$\begin{array}{}\text{(59)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\\ \mathrm{A}\left(t\right)\\ \mathrm{B}\left(t\right)\\ \mathrm{C}\left(t\right)\\ \mathrm{D}\left(t\right)\end{array}\right]=\\ & \left[\begin{array}{cccccc}-k_{\mathrm{1}}\mathrm{C}\left(t\right)-R_{\mathrm{1}}^{\mathrm{A}}& k_{-\mathrm{1}}\mathrm{D}\left(t\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ k_{\mathrm{1}}\mathrm{C}\left(t\right)& -k_{-\mathrm{1}}\mathrm{D}\left(t\right)-R_{\mathrm{1}}^{\mathrm{B}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ R_{\mathrm{1}}^{\mathrm{A}}& \mathrm{0}& -k_{\mathrm{1}}\mathrm{C}\left(t\right)& k_{-\mathrm{1}}\mathrm{D}\left(t\right)& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& R_{\mathrm{1}}^{\mathrm{B}}& k_{\mathrm{1}}\mathrm{C}\left(t\right)& -k_{-\mathrm{1}}\mathrm{D}\left(t\right)& \mathrm{0}& \mathrm{0}\\ -k_{\mathrm{1}}\mathrm{C}\left(t\right)& k_{-\mathrm{1}}\mathrm{D}\left(t\right)& -k_{\mathrm{1}}\mathrm{C}\left(t\right)& k_{-\mathrm{1}}\mathrm{D}\left(t\right)& \mathrm{0}& \mathrm{0}\\ k_{\mathrm{1}}\mathrm{C}\left(t\right)& -k_{-\mathrm{1}}\mathrm{D}\left(t\right)& k_{\mathrm{1}}\mathrm{C}\left(t\right)& -k_{-\mathrm{1}}\mathrm{D}\left(t\right)& \mathrm{0}& \mathrm{0}\end{array}\right]\\ & \left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\\ \mathrm{A}\left(t\right)\\ \mathrm{B}\left(t\right)\\ \mathrm{C}\left(t\right)\\ \mathrm{D}\left(t\right)\end{array}\right].\end{array}\end{array}$$

We can apply a similarity transform given by

$$\begin{array}{}\text{(60)}& U=\left[\begin{array}{cccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],\end{array}$$

to yield an equation of motion in the transformed basis vector given by

$$\begin{array}{}\text{(61)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\\ {\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)\\ \mathrm{C}\left(t\right)\\ \mathrm{D}\left(t\right)\end{array}\right]=\\ & \left[\begin{array}{cccccc}-k_{\mathrm{1}}\mathrm{C}\left(t\right)-R_{\mathrm{1}}^{\mathrm{A}}& k_{-\mathrm{1}}\mathrm{D}\left(t\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ k_{\mathrm{1}}\mathrm{C}\left(t\right)& k_{-\mathrm{1}}\mathrm{D}\left(t\right)-R_{\mathrm{1}}^{\mathrm{B}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}\mathrm{C}\left(t\right)& k_{-\mathrm{1}}\mathrm{D}\left(t\right)& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& k_{\mathrm{1}}\mathrm{C}\left(t\right)& -k_{-\mathrm{1}}\mathrm{D}\left(t\right)& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}\mathrm{C}\left(t\right)& k_{-\mathrm{1}}\mathrm{D}\left(t\right)& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& k_{\mathrm{1}}\mathrm{C}\left(t\right)& -k_{-\mathrm{1}}\mathrm{D}\left(t\right)& \mathrm{0}& \mathrm{0}\end{array}\right]\\ & \left[\begin{array}{c}{\mathrm{A}}^{\ast }\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)\\ {\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)\\ {\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)\\ \mathrm{C}\left(t\right)\\ \mathrm{D}\left(t\right)\end{array}\right].\end{array}\end{array}$$

Figure 6 shows numerical simulations of the time evolution of the system of Eqs. (53)–(58) with initial magnetization corresponding to the hyperpolarized signal ${\mathrm{A}}^{\ast }\left(\mathrm{0}\right)=\mathrm{1}$ and non-polarized substrates C(0)=0.95 and D(0)=0.05. The longitudinal relaxation rate constants were $R_{\mathrm{1}}^{\mathrm{A}}=R_{\mathrm{1}}^{\mathrm{B}}=\mathrm{1}/\mathrm{60}$ s⁻¹. The time dependences of A^∗(t), A(t), B^∗(t) and B(t) are subject to second-order kinetics and were calculated numerically (left column of plots) for different rate constants: Fig. 6b, k₁=0.01 M⁻¹ s⁻¹, $k_{-\mathrm{1}}=\mathrm{0}$ M⁻¹ s⁻¹, corresponding to uni-directional kinetics; Fig. 6c, k₁=0.01 M⁻¹ s⁻¹, $k_{-\mathrm{1}}=\mathrm{0.005}$ M⁻¹ s⁻¹, corresponding to bi-directional exchange kinetics with an equilibrium constant K=2; and Fig. 6d $k_{\mathrm{1}}=k_{-\mathrm{1}}=\mathrm{0.01}$ M⁻¹ s⁻¹, with an equilibrium constant K=1, also corresponding to bi-directional exchange kinetics. The right column of plots shows the time dependence of ${\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)$, ${\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)$, which capture conventional chemical kinetics of the concentrations of [A(t)] and [B(t)], as expected, as well as the kinetics of the non-polarized reactants [C(t)] and [D(t)].

4.3.1 An ersatz solution

The system of differential equations in Eq. (59) describing a second-order reaction can be reduced to one with pseudo first-order kinetics by introducing time-dependent rate constants $k_{\mathrm{1}}^{\prime }\left(t\right)=k_{\mathrm{1}}\mathrm{C}\left(t\right)$ and $k_{-\mathrm{1}}^{\prime }\left(t\right)=k_{-\mathrm{1}}\mathrm{D}\left(t\right)$. Importantly, the pseudo first-order rate constants $k_{\mathrm{1}}^{\prime }\left(t\right)$ and $k_{-\mathrm{1}}^{\prime }\left(t\right)$ are now time-dependent. This approach has been used previously (Mariotti et al., 2016), but it constitutes a special case of the more general method described here, which we now advocate.

However, we now encounter a problem. The pseudo first-order rate constants for the reactions of [C(t)] and [D(t)] are now given by $k_{\mathrm{1}}^{\prime }\left(t\right)=k_{\mathrm{1}}\left({\mathrm{A}}^{\ast }\left(t\right)+\mathrm{A}\left(t\right)\right)$ and $k_{-\mathrm{1}}^{\prime }\left(t\right)=k_{-\mathrm{1}}\left({\mathrm{B}}^{\ast }\left(t\right)+\mathrm{B}\left(t\right)\right)$, respectively. The time-dependent pseudo first-order rate constants are dependent on the concentrations of both “hot” and “cold” pools. In turn, the pseudo first-order rate constants for A^∗(t) and B^∗(t) are $k_{\mathrm{1}}^{\prime }\left(t\right)=k_{\mathrm{1}}\mathrm{C}\left(t\right)$ and $k_{-\mathrm{1}}^{\prime }\left(t\right)=k_{-\mathrm{1}}\mathrm{D}\left(t\right)$. Thus, the kinetics of the “hot” pools A^∗(t) and B^∗(t) become dependent on the kinetics of the “cold” pools A(t) and B(t). This is of particular relevance (as highlighted by Kuchel and Shishmarev, 2020) when extending the equations to describe enzyme kinetics. It is this that we turn our attention to next.

5 Michaelis–Menten equation for a hyperpolarized substrate

Next consider an enzyme-catalysed reaction with a hyperpolarized substrate. The simplest model involves a hyperpolarized substrate S^∗(t) that is in equilibrium with a free enzyme of concentration [E]₀ to form a hyperpolarized enzyme–substrate complex ES^∗(t), which then reacts to form a hyperpolarized product P^∗(t). This is followed by release of the free enzyme that is then available for further reactions: $\mathrm{E}+{\mathrm{S}}^{\ast }\leftrightarrow {\mathrm{ES}}^{\ast }\leftrightarrow {\mathrm{P}}^{\ast }+\mathrm{E}$. All hyperpolarized substrates relax through T₁ processes to form non-polarized pools of substrates $\mathrm{E}+\mathrm{S}\leftrightarrow \mathrm{ES}\leftrightarrow \mathrm{P}+\mathrm{E}$ as shown in Fig. 7a. The differential equations (again omitting the square brackets denoting concentration) that describe the reaction kinetics are

$$\begin{array}{}\text{(62)}& {\displaystyle \frac{{\mathrm{dS}}^{\ast }\left(t\right)}{\mathrm{d}t}}=-k_{\mathrm{1}}\mathrm{E}\left(t\right){\mathrm{S}}^{\ast }\left(t\right)+k_{-\mathrm{1}}{\mathrm{ES}}^{\ast }\left(t\right)-R_{\mathrm{1}}^{\mathrm{S}}{\mathrm{S}}^{\ast }\left(t\right),\text{(63)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{{\mathrm{dES}}^{\ast }\left(t\right)}{\mathrm{d}t}}& =k_{\mathrm{1}}\mathrm{E}\left(t\right){\mathrm{S}}^{\ast }\left(t\right)-k_{-\mathrm{1}}{\mathrm{ES}}^{\ast }\left(t\right)-k_{\mathrm{2}}{\mathrm{ES}}^{\ast }\left(t\right)\\ & +k_{-\mathrm{2}}\mathrm{E}\left(t\right){\mathrm{P}}^{\ast }\left(t\right)-R_{\mathrm{1}}^{\mathrm{ES}}{\mathrm{ES}}^{\ast }\left(t\right),\end{array}\text{(64)}& {\displaystyle \frac{{\mathrm{dP}}^{\ast }\left(t\right)}{\mathrm{d}t}}=k_{\mathrm{2}}{\mathrm{ES}}^{\ast }\left(t\right)-k_{-\mathrm{2}}\mathrm{E}\left(t\right){\mathrm{P}}^{\ast }\left(t\right)-R_{\mathrm{1}}^{\mathrm{P}}{\mathrm{P}}^{\ast }\left(t\right),\text{(65)}& {\displaystyle \frac{\mathrm{dS}\left(t\right)}{\mathrm{d}t}}=-k_{\mathrm{1}}\mathrm{E}\left(t\right)\mathrm{S}\left(t\right)+k_{-\mathrm{1}}\mathrm{ES}\left(t\right)+R_{\mathrm{1}}^{\mathrm{S}}{\mathrm{S}}^{\ast }\left(t\right),\text{(66)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dES}\left(t\right)}{\mathrm{d}t}}& =k_{\mathrm{1}}\mathrm{E}\left(t\right)\mathrm{S}\left(t\right)-k_{-\mathrm{1}}\mathrm{ES}\left(t\right)-k_{\mathrm{2}}\mathrm{ES}\left(t\right)\\ & +k_{-\mathrm{2}}\mathrm{E}\left(t\right)\mathrm{P}\left(t\right)+R_{\mathrm{1}}^{\mathrm{ES}}{\mathrm{ES}}^{\ast }\left(t\right),\end{array}\text{(67)}& {\displaystyle \frac{\mathrm{dP}\left(t\right)}{\mathrm{d}t}}=k_{\mathrm{2}}\mathrm{ES}\left(t\right)-k_{-\mathrm{2}}\mathrm{E}\left(t\right)\mathrm{P}\left(t\right)+R_{\mathrm{1}}^{\mathrm{P}}{\mathrm{P}}^{\ast }\left(t\right),\text{(68)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dE}\left(t\right)}{\mathrm{d}t}}& =-k_{\mathrm{1}}\mathrm{E}\left(t\right)\left({\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)\right)+\left(k_{-\mathrm{1}}+k_{\mathrm{2}}\right)\left({\mathrm{ES}}^{\ast }\left(t\right)\right.\\ & \left.+\mathrm{ES}\left(t\right)\right)-k_{-\mathrm{2}}\mathrm{E}\left(t\right)\left({\mathrm{P}}^{\ast }\left(t\right)+\mathrm{P}\left(t\right)\right),\end{array}\end{array}$$

where E(t) is the free enzyme, ES(t) is the enzyme–substrate complex, S(t) is the free substrate and P(t) is the free product, with relaxation rate constants $R_{\mathrm{1}}^{\mathrm{S}}$, $R_{\mathrm{1}}^{\mathrm{ES}}$ and $R_{\mathrm{1}}^{\mathrm{P}}$, respectively. Note the appearance of the free enzyme E(t) as both a reactant and product; it is regenerated through the reactions that are characterized by the rate constants k₁ and k₋₁ and also k₂ and k₋₂, thereby being recycled.

Figure 7 Simulated Michaelis–Menten kinetics for exchange of hyperpolarized solutes $\mathrm{E}+{\mathrm{S}}^{\ast }\leftrightarrow {\mathrm{ES}}^{\ast }\leftrightarrow {\mathrm{P}}^{\ast }+\mathrm{E}$ conforming to conservation of mass, assuming initial hyperpolarized magnetization of only solute ${\mathrm{S}}^{\ast }\left(\mathrm{0}\right)=\mathrm{0.001}$ and ${\left[\mathrm{E}\right]}_{\mathrm{0}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{9}}$ M. Longitudinal relaxation rate constants were $R_{\mathrm{1}}^{\mathrm{S}}=R_{\mathrm{1}}^{\mathrm{ES}}=R_{\mathrm{1}}^{\mathrm{P}}=\mathrm{1}/\mathrm{60}$ s⁻¹. The reaction rate constants were $k_{\mathrm{1}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{7}}$ M⁻¹ s⁻¹, $k_{-\mathrm{1}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{2}}$ s⁻¹, $k_{\mathrm{2}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ s⁻¹ and $k_{-\mathrm{2}}=\mathrm{0}$ M⁻¹ s⁻¹, such that $K_{\mathrm{M}}=\mathrm{5.1}\times {\mathrm{10}}^{-\mathrm{4}}$ M and $V_{\max }=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$ M s⁻¹. Left panels: (b) simulated time dependence of S^∗(t), S(t), P^∗(t) and P(t), and (c) simulated time dependence of ES^∗(t) and ES(t). Right panels: (b) simulated time dependence of ${\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)=\left[\mathrm{S}\left(t\right)\right]$ and ${\mathrm{P}}^{\ast }\left(t\right)+\mathrm{P}\left(t\right)=\left[\mathrm{P}\left(t\right)\right]$, and (c) ${\mathrm{ES}}^{\ast }\left(t\right)+\mathrm{ES}\left(t\right)=\left[\mathrm{ES}\left(t\right)\right]$ and [E(t)].

Mass is conserved as confirmed by the fact that $\mathrm{d}\left({\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)+{\mathrm{ES}}^{\ast }\left(t\right)+\mathrm{ES}\left(t\right)+{\mathrm{P}}^{\ast }\left(t\right)+\mathrm{P}\left(t\right)\right)/\mathrm{d}t=\mathrm{0}$ and $\mathrm{d}\left({\mathrm{ES}}^{\ast }\left(t\right)+\mathrm{ES}\left(t\right)+\mathrm{E}\left(t\right)\right)/\mathrm{d}t=\mathrm{0}$. Therefore, provided ${\mathrm{S}}^{\ast }\left(\mathrm{0}\right)+\mathrm{S}\left(\mathrm{0}\right)={\left[\mathrm{S}\right]}_{\mathrm{0}}$, ${\mathrm{ES}}^{\ast }\left(\mathrm{0}\right)+\mathrm{ES}\left(\mathrm{0}\right)={\left[\mathrm{ES}\right]}_{\mathrm{0}}$ and ${\mathrm{P}}^{\ast }\left(\mathrm{0}\right)+\mathrm{P}\left(\mathrm{0}\right)={\left[\mathrm{P}\right]}_{\mathrm{0}}$, then ${\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)=\left[\mathrm{S}\left(t\right)\right]$, ${\mathrm{ES}}^{\ast }\left(t\right)+\mathrm{ES}\left(t\right)=\left[\mathrm{ES}\left(t\right)\right]$ and ${\mathrm{P}}^{\ast }\left(t\right)+\mathrm{P}\left(t\right)=\left[\mathrm{P}\left(t\right)\right]$, respectively.

Equations (62)–(68) can be written in matrix vector form given by Eq. (A69); see Appendix. We can apply a similarity transform, given by Eq. (A70) (see Appendix), to yield an equation of motion in the transformed basis vector given by Eq. (A71); see Appendix.

5.1 Steady state of ES complex

A simplified uni-directional enzyme-catalysed reaction is described by setting the reverse rate constant $k_{-\mathrm{2}}=\mathrm{0}$ (see Fig. 7a). If it is assumed that a steady state of [ES] is attained very rapidly, then $\mathrm{d}\left({\mathrm{ES}}^{\ast }\left(t\right)+\mathrm{ES}\left(t\right)\right)/\mathrm{d}t=\mathrm{0}$, and we obtain (reverting to using square brackets to denote molar concentration)

$$\begin{array}{}\text{(72)}& k_{\mathrm{1}}\left[\mathrm{E}\left(t\right)\right]\left[{\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)\right]=\left(k_{-\mathrm{1}}+k_{\mathrm{2}}\right)\left[{\mathrm{ES}}^{\ast }\left(t\right)+\mathrm{ES}\left(t\right)\right].\end{array}$$

Rearranging Eq. (72) yields the Michaelis constant in terms of hyperpolarized and non-polarized pools of substrate:

$$\begin{array}{}\text{(73)}& K_{\mathrm{M}}={\displaystyle \frac{\left(k_{-\mathrm{1}}+k_{\mathrm{2}}\right)}{k_{\mathrm{1}}}}={\displaystyle \frac{\left[\mathrm{E}\left(t\right)\right]\left[{\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)\right]}{\left[{\mathrm{ES}}^{\ast }\left(t\right)+\mathrm{ES}\left(t\right)\right]}}.\end{array}$$

Calibrating the signals to molar concentrations is important since the signals now relate to a real parameter (K_M) of the enzyme that has units of concentration (typically mM).

Thus, using conservation of enzyme mass, the free enzyme concentration is given by

$$\begin{array}{}\text{(74)}& \left[\mathrm{E}\right(t\left)\right]=[\mathrm{E}{]}_{\mathrm{0}}-\left[{\mathrm{ES}}^{\ast }\left(t\right)+\mathrm{ES}\left(t\right)\right].\end{array}$$

Then

$$\begin{array}{}\text{(75)}& {\displaystyle \frac{\mathrm{d}\left(\left[{\mathrm{P}}^{\ast }\left(t\right)+\mathrm{P}\left(t\right)\right]\right)}{\mathrm{d}t}}={\displaystyle \frac{k_{\mathrm{2}}\left[\mathrm{E}{]}_{\mathrm{0}}\right[{\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)]}{K_{\mathrm{M}}+\left[{\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)\right]}}.\end{array}$$

In other words, this is the standard form of the Michaelis–Menten equation written as a function of both polarized and unpolarized pools of substrate.

5.2 Simulations of the Michaelis–Menten reaction

Figure 7b–c show the results of numerical integration of Eqs. (62)–(68) with an initial hyperpolarized signal ${\mathrm{S}}^{\ast }\left(\mathrm{0}\right)=\mathrm{0.001}$ (corresponding to a concentration [S]₀=1 mM via the experimentally determined scaling factor, which here was set to 1) and enzyme concentration ${\left[\mathrm{E}\right]}_{\mathrm{0}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{9}}$ M. The assigned longitudinal relaxation rate constants were $R_{\mathrm{1}}^{\mathrm{S}}=R_{\mathrm{1}}^{\mathrm{ES}}=R_{\mathrm{1}}^{\mathrm{P}}=\mathrm{1}/\mathrm{60}$ s⁻¹. In the first instance, we set the longitudinal relaxation times of substrate, enzyme–substrate complex and product to be equal (this is discussed further below). The reaction rate constants were $k_{\mathrm{1}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{7}}$ M⁻¹ s⁻¹, $k_{-\mathrm{1}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{2}}$ s⁻¹, $k_{\mathrm{2}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ s⁻¹, and $k_{-\mathrm{2}}=\mathrm{0}$ M⁻¹ s⁻¹, such that $K_{\mathrm{M}}=\mathrm{5.1}\times {\mathrm{10}}^{-\mathrm{4}}$ M and $V_{\max }=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$ M s⁻¹. The time dependences of S^∗(t), S(t), P^∗(t) and P(t) are shown in Fig. 7b, left panel, subject to standard uni-directional Michaelis–Menten kinetics, and in Fig. 7c, left panel, the time dependence of ES^∗(t) and ES(t). The time dependence of ${\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)=\left[\mathrm{S}\left(t\right)\right]$ and ${\mathrm{P}}^{\ast }\left(t\right)+\mathrm{P}\left(t\right)=\left[\mathrm{P}\left(t\right)\right]$ are shown in Fig. 7b, right panel, and ${\mathrm{ES}}^{\ast }\left(t\right)+\mathrm{ES}\left(t\right)=\left[\mathrm{ES}\left(t\right)\right]$ and [E(t)] are shown in Fig. 7c, right panel, which recapture conventional chemical kinetics of [S(t)], [ES(t)], [P(t)] and [E(t)], as required for mathematical and physical consistency.

It is worth considering some of the consequences of Eq. (75) when studying enzyme-mediated reactions with hyperpolarized substrates. When the substrate concentration $\left[{\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)\right]$ is much greater than K_M, then the rate of product formation $\mathrm{d}\left(\left[{\mathrm{P}}^{\ast }\left(t\right)+\mathrm{P}\left(t\right)\right]\right)/\mathrm{d}t$ is given by $v=k_{\mathrm{2}}{\left[\mathrm{E}\right]}_{\mathrm{0}}=V_{\max }$, which is constant (i.e. it is effectively a zero-order reaction with respect to substrate concentration). The enzyme is said to be saturated; its rate is independent of substrate concentration, but V_max is proportional to the enzyme concentration [E]₀. When the substrate concentration $\left[{\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)\right]$ is much less than K_M, then the rate of product formation $\mathrm{d}\left(\left[{\mathrm{P}}^{\ast }\left(t\right)+\mathrm{P}\left(t\right)\right]\right)/\mathrm{d}t$ is given by $V=k_{\mathrm{2}}{\left[\mathrm{E}\right]}_{\mathrm{0}}\left[{\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)\right]/K_{\mathrm{M}}$ and the reaction is effectively first order with respect to substrate concentration. Nevertheless, the rate is still proportional to [E]₀ as well. The kinetics of enzyme systems, and indeed enzyme kinetics in general, are a composite of the two parameters K_M and V_max. The influences on one cannot be distinguished from the other on the basis of time course experiments alone; separate measurements are needed to estimate the total enzyme concentration.

Further simulations were performed to explore the influence of a much shorter value of $T_{\mathrm{1}}^{\mathrm{ES}}$ for the enzyme–substrate complex, while $T_{\mathrm{1}}^{\mathrm{S}}$ and $T_{\mathrm{1}}^{\mathrm{P}}$ were unchanged. Even if it were assumed to be very small viz. $T_{\mathrm{1}}^{\mathrm{ES}}=\mathrm{276.4}$ ms, the time evolution was indistinguishable from that presented in Fig. 7; the corresponding curves were superimposable. The signal that resided on the enzyme–substrate complex ES^∗ was 6 orders of magnitude lower than that of the substrate S^∗ and product P^∗. Therefore, the kinetics of signal evolution were dominated by $T_{\mathrm{1}}^{\mathrm{S}}$ and $T_{\mathrm{1}}^{\mathrm{P}}$, while changes in $T_{\mathrm{1}}^{\mathrm{ES}}$ could be ignored. An exception to this analysis might occur if the active site were next to a paramagnetic centre, such as is found in metalloproteins for which $T_{\mathrm{1}}^{\mathrm{ES}}$ could be very much shorter than predicted (see the relaxation theory section above).

5.3 Enzyme inhibition and hyperpolarized substrate kinetics

Our formalism can be readily extended to account for the influence of a ligand/solute to inhibit an enzyme. The simplest case is when a solute binds reversibly to the free enzyme E to form an enzyme–inhibitor complex EI; hence, the enzyme becomes unable to bind and react with its substrate S. To describe this scenario, Eq. (68) is modified to include an additional pathway for the loss of free enzyme:

$$\begin{array}{}\text{(76)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left[\mathrm{E}\left(t\right)\right]}{\mathrm{d}t}}& =-k_{\mathrm{1}}\left[\mathrm{E}\left(t\right)\right]\left[{\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)\right]\\ & +\left(k_{-\mathrm{1}}+k_{\mathrm{2}}\right)\left[{\mathrm{ES}}^{\ast }\left(t\right)+\mathrm{ES}\left(t\right)\right]\\ & -k_{-\mathrm{2}}\left[\mathrm{E}\left(t\right)\right]\left[{\mathrm{P}}^{\ast }\left(t\right)+\mathrm{P}\left(t\right)\right]\\ & -k_{\mathrm{3}}\left[\mathrm{E}\left(t\right)\right]\left[\mathrm{I}\left(t\right)\right]+k_{-\mathrm{3}}\left[\mathrm{EI}\left(t\right)\right].\end{array}\end{array}$$

The model is now extended to include differential equations describing the concentration of the inhibitor [I(t)] and the enzyme–inhibitor complex [EI(t)]:

$$\begin{array}{}\text{(77)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}\left[\mathrm{I}\left(t\right)\right]}{\mathrm{d}t}}=-k_{\mathrm{3}}\left[\mathrm{E}\left(t\right)\right]\left[\mathrm{I}\left(t\right)\right]+k_{-\mathrm{3}}\left[\mathrm{EI}\left(t\right)\right],\\ & {\displaystyle \frac{\mathrm{d}\left[\mathrm{EI}\left(t\right)\right]}{\mathrm{d}t}}=k_{\mathrm{3}}\left[\mathrm{E}\left(t\right)\right]\left[\mathrm{I}\left(t\right)\right]-k_{-\mathrm{3}}\left[\mathrm{EI}\left(t\right)\right].\end{array}\end{array}$$

Such equations can be incorporated into the Michaelis–Menten equations, and we develop this next.

5.3.1 Types of enzyme inhibition

There are three commonly encountered types of reversible enzyme inhibition (Kuchel, 2009): (i) a competitive inhibitor is structurally similar to the substrate and binds preferentially in the active site of the free enzyme, E, thus preventing the substrate from binding and reacting; (ii) an uncompetitive inhibitor binds only to the enzyme–substrate complex and therefore causes substrate-concentration-dependent inhibition; and (iii), a non-competitive inhibitor binds to both the free enzyme and to the enzyme–substrate complex; it causes a conformational change at the active site that inhibits (or even enhances) the reaction. Such an effect is referred to as allosteric inhibition (or activation).

Accounting for all three scenarios, the free enzyme concentration is given by

$$\begin{array}{}\text{(78)}& \begin{array}{rl}\left[\mathrm{E}\right(t\left)\right]& =[\mathrm{E}{]}_{\mathrm{0}}-[\mathrm{EI}\left(t\right)]-[{\mathrm{ES}}^{\ast }\left(t\right)+\mathrm{ES}\left(t\right)]\\ & -\left[{\mathrm{ESI}}^{\ast }\left(t\right)+\mathrm{ESI}\left(t\right)\right].\end{array}\end{array}$$

Substituting

$$\begin{array}{}\text{(79)}& \mathit{\alpha }=\mathrm{1}+{\displaystyle \frac{\left[\mathrm{I}\left(t\right)\right]}{K_{\mathrm{I}}}}\mathrm{and}{\mathit{\alpha }}^{\prime }=\mathrm{1}+{\displaystyle \frac{\left[\mathrm{I}\left(t\right)\right]}{K_{\mathrm{I}}^{\prime }}},\end{array}$$

where $K_{\mathrm{I}}=\left[\mathrm{E}\left(t\right)\right]\left[\mathrm{I}\left(t\right)\right]/\left[\mathrm{EI}\left(t\right)\right]$ and $K_{\mathrm{I}}^{\prime }=\left[\mathrm{ES}\left(t\right)\right]\left[\mathrm{I}\left(t\right)\right]/\left[\mathrm{ESI}\left(t\right)\right]$, yields

$$\begin{array}{}\text{(80)}& {\displaystyle \frac{\mathrm{d}\left(\left[{\mathrm{P}}^{\ast }\left(t\right)+\mathrm{P}\left(t\right)\right]\right)}{\mathrm{d}t}}={\displaystyle \frac{k_{\mathrm{2}}\left[\mathrm{E}{]}_{\mathrm{0}}\right[{\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)]}{\mathit{\alpha }{K}_{\mathrm{M}}+{\mathit{\alpha }}^{\prime }\left[{\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)\right]}}.\end{array}$$

The three types of enzyme inhibition can be distinguished by their influence on the kinetic parameters that are estimated in specially designed experiments performed on the enzyme over a range of substrate and inhibitor concentrations (Kuchel, 2009): (i) competitive inhibitors cause an increase in apparent K_M value, while V_max is unchanged; (ii) uncompetitive inhibitors cause a reduction in V_max, while the apparent K_M is unchanged; and (iii) non-competitive inhibitors cause both a reduction in V_max and an increase in apparent K_M.

An additional effect that can be considered is where either the substrate of the reaction [S(t)] or the product of the reaction, [P(t)], acts as the inhibitor, called unsurprisingly “substrate inhibition” and “product inhibition”, respectively. The relevant enzyme kinetic equations are composed by substituting $\left[\mathrm{I}\left(t\right)\right]=\left[{\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)\right]$ or $\left[\mathrm{I}\left(t\right)\right]=\left[{\mathrm{P}}^{\ast }\left(t\right)+\mathrm{P}\left(t\right)\right]$ in the above equations.

6 Cofactors and unlabelled pools – lactate dehydrogenase

We now consider a real system that is of contemporary interest for in vivo clinical studies using dDNP. It is lactate dehydrogenase (LDH; E.C. 1.1.1.27). Consider the LDH-catalysed reaction of a hyperpolarized substrate; it follows an ordered sequential reaction in which $\mathrm{E}+\mathrm{NADH}\leftrightarrow \mathrm{E}\cdot \mathrm{NADH}+{\mathrm{Pyr}}^{\ast }\leftrightarrow \mathrm{E}\cdot \mathrm{NAD}+{\mathrm{Lac}}^{\ast }\leftrightarrow \mathrm{E}+{\mathrm{NAD}}^{+}$. Again, we assume that relaxation of magnetization occurs through T₁ processes to form a pool of reactants $\mathrm{E}+\mathrm{NADH}\leftrightarrow \mathrm{E}\cdot \mathrm{NADH}+\mathrm{Pyr}\leftrightarrow \mathrm{E}\cdot \mathrm{NAD}+\mathrm{Lac}\leftrightarrow \mathrm{E}+{\mathrm{NAD}}^{+}$ as shown in Fig. 8a. The relevant differential equations used to describe the kinetics are (omitting the square brackets that denote concentration)

$$\begin{array}{}\text{(81)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{{\mathrm{dPyr}}^{\ast }\left(t\right)}{\mathrm{d}t}}& =-k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right){\mathrm{Pyr}}^{\ast }\left(t\right)\\ & +k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right){\mathrm{Lac}}^{\ast }\left(t\right)-R_{\mathrm{1}}^{\mathrm{P}}{\mathrm{Pyr}}^{\ast }\left(t\right),\end{array}\text{(82)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{{\mathrm{dLac}}^{\ast }\left(t\right)}{\mathrm{d}t}}& =k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right){\mathrm{Pyr}}^{\ast }\left(t\right)\\ & -k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right){\mathrm{Lac}}^{\ast }\left(t\right)-R_{\mathrm{1}}^{\mathrm{L}}{\mathrm{Lac}}^{\ast }\left(t\right),\end{array}\text{(83)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dPyr}\left(t\right)}{\mathrm{d}t}}& =-k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)\mathrm{Pyr}\left(t\right)\\ & +k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)\mathrm{Lac}\left(t\right)+R_{\mathrm{1}}^{\mathrm{P}}{\mathrm{Pyr}}^{\ast }\left(t\right),\end{array}\text{(84)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dLac}\left(t\right)}{\mathrm{d}t}}& =k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)\mathrm{Pyr}\left(t\right)\\ & -k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)\mathrm{Lac}\left(t\right)+R_{\mathrm{1}}^{\mathrm{L}}{\mathrm{Lac}}^{\ast }\left(t\right),\end{array}\text{(85)}& {\displaystyle \frac{\mathrm{dNADH}\left(t\right)}{\mathrm{d}t}}=-k_{\mathrm{1}}\mathrm{E}\left(t\right)\mathrm{NADH}\left(t\right)+k_{-\mathrm{1}}\mathrm{E}.\mathrm{NADH}\left(t\right),\text{(86)}& {\displaystyle \frac{\mathrm{dNAD}\left(t\right)}{\mathrm{d}t}}=k_{\mathrm{3}}\mathrm{E}.\mathrm{NAD}\left(t\right)-k_{-\mathrm{3}}\mathrm{E}\left(t\right)\mathrm{NAD}\left(t\right),\text{(87)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dE}.\mathrm{NADH}\left(t\right)}{\mathrm{d}t}}& =k_{\mathrm{1}}\mathrm{E}\left(t\right)\mathrm{NADH}\left(t\right)-k_{-\mathrm{1}}\mathrm{E}.\mathrm{NADH}\left(t\right)\\ & -k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)({\mathrm{Pyr}}^{\ast }\left(t\right)+\mathrm{Pyr}(t\left)\right)\\ & +k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)({\mathrm{Lac}}^{\ast }\left(t\right)+\mathrm{Lac}(t\left)\right),\end{array}\text{(88)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dE}.\mathrm{NAD}\left(t\right)}{\mathrm{d}t}}& =k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)({\mathrm{Pyr}}^{\ast }\left(t\right)+\mathrm{Pyr}(t\left)\right)\\ & -k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)({\mathrm{Lac}}^{\ast }\left(t\right)+\mathrm{Lac}(t\left)\right)\\ & -k_{\mathrm{3}}\mathrm{E}.\mathrm{NAD}\left(t\right)+k_{-\mathrm{3}}\mathrm{E}\left(t\right)\mathrm{NAD}\left(t\right),\end{array}\text{(89)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dE}\left(t\right)}{\mathrm{d}t}}& =-k_{\mathrm{1}}\mathrm{E}\left(t\right)\mathrm{NADH}\left(t\right)+k_{-\mathrm{1}}\mathrm{E}.\mathrm{NADH}\left(t\right)\\ & +k_{\mathrm{3}}\mathrm{E}.\mathrm{NAD}\left(t\right)-k_{-\mathrm{3}}\mathrm{E}\left(t\right)\mathrm{NAD}\left(t\right),\end{array}\end{array}$$

where E(t) is the concentration of free enzyme, NAD(t) and NADH(t) are the concentrations of the free cofactors, E.NAD(t) and E.NADH(t) are the concentrations of the enzyme–cofactor complexes and Pyr(t) and Lac(t) are the free substrates with relaxation rate constants $R_{\mathrm{1}}^{\mathrm{P}}$ and $R_{\mathrm{1}}^{\mathrm{L}}$, respectively.

Figure 8 Simulated kinetics of LDH for exchange of solutes, $\mathrm{E}+\mathrm{NADH}\leftrightarrow \mathrm{E}\cdot \mathrm{NADH}+{\mathrm{Pyr}}^{\ast }\leftrightarrow \mathrm{E}\cdot \mathrm{NAD}+{\mathrm{Lac}}^{\ast }\leftrightarrow \mathrm{E}+{\mathrm{NAD}}^{+}$, conforming to conservation of mass, assuming initial hyperpolarized magnetization of only solute ${\mathrm{Pyr}}^{\ast }\left(\mathrm{0}\right)=\mathrm{0.001}$ and ${\left[\mathrm{E}\right]}_{\mathrm{0}}=\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{9}}$ M. Longitudinal relaxation rate constants were $R_{\mathrm{1}}^{\mathrm{P}}=R_{\mathrm{1}}^{\mathrm{L}}=\mathrm{1}/\mathrm{60}$ s⁻¹. Rate constants were $k_{\mathrm{1}}=\mathrm{1.03}\times {\mathrm{10}}^{\mathrm{8}}$ M⁻¹ s⁻¹, $k_{-\mathrm{1}}=\mathrm{549}$ s⁻¹, $k_{\mathrm{2}}=\mathrm{6.72}\times {\mathrm{10}}^{\mathrm{6}}$ M⁻¹ s⁻¹, $k_{-\mathrm{2}}=\mathrm{3.44}\times {\mathrm{10}}^{\mathrm{4}}$ M⁻¹ s⁻¹, k₃=842 s⁻¹ and $k_{-\mathrm{3}}=\mathrm{9.12}\times {\mathrm{10}}^{\mathrm{5}}$ M⁻¹ s⁻¹. Initial cofactor concentrations were $\left[\mathrm{NADH}\left(\mathrm{0}\right)\right]=\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}$ M and $\left[\mathrm{NAD}\left(\mathrm{0}\right)\right]=\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{3}}$ M. (b) Simulated time dependence Pyr^∗(t), Pyr(t), Lac^∗(t) and Lac(t), left panel, [E(t)], [E.NAD(t)] and [E.NADH(t)], middle panel, and ${\mathrm{Pyr}}^{\ast }\left(t\right)+\mathrm{Pyr}\left(t\right)=\left[\mathrm{Pyr}\left(t\right)\right]$, ${\mathrm{Lac}}^{\ast }\left(t\right)+\mathrm{Lac}\left(t\right)=\left[\mathrm{Lac}\left(t\right)\right]$, [NAD(t)] and [NADH(t)], right panel. (c) Simulations of the time dependence of Lac^∗(t) under the conditions that [E]₀= (i) $\mathrm{0.6}\times {\mathrm{10}}^{-\mathrm{9}}$ M; (ii) $\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{9}}$ M; and (iii) $\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{9}}$ M, while all other parameters remained unchanged. (d) Simulations of the time dependence of Lac^∗(t) under the conditions that Lac(0)= (i) 0 mM; (ii) 20 mM; and (iii) 40 mM, while all other parameters remained unchanged.

Mass is conserved, as is confirmed by the fact that $\mathrm{d}\left({\mathrm{Pyr}}^{\ast }\left(t\right)+\mathrm{Pyr}\left(t\right)+{\mathrm{Lac}}^{\ast }\left(t\right)+\mathrm{Lac}\left(t\right)\right)/\mathrm{d}t=\mathrm{0}$. Enzyme concentration is conserved, as is confirmed by $\mathrm{d}\left(\mathrm{E}.\mathrm{NADH}\left(t\right)+\mathrm{E}.\mathrm{NAD}\left(t\right)+\mathrm{E}\left(t\right)\right)/\mathrm{d}t=\mathrm{0}$, and cofactor pools are conserved, as is confirmed by $\mathrm{d}\left(\mathrm{NADH}\left(t\right)+\mathrm{NAD}\left(t\right)+\mathrm{E}.\mathrm{NADH}\left(t\right)+\mathrm{E}.\mathrm{NAD}\left(t\right)\right)/\mathrm{d}t=\mathrm{0}$. Therefore, provided ${\mathrm{Pyr}}^{\ast }\left(\mathrm{0}\right)+\mathrm{Pyr}\left(\mathrm{0}\right)={\left[\mathrm{Pyr}\right]}_{\mathrm{0}}$ and ${\mathrm{Lac}}^{\ast }\left(\mathrm{0}\right)+\mathrm{Lac}\left(\mathrm{0}\right)={\left[\mathrm{Lac}\right]}_{\mathrm{0}}$, then ${\mathrm{Pyr}}^{\ast }\left(t\right)+\mathrm{Pyr}\left(t\right)=\left[\mathrm{Pyr}\left(t\right)\right]$ and ${\mathrm{Lac}}^{\ast }\left(t\right)+\mathrm{Lac}\left(t\right)=\left[\mathrm{Lac}\left(t\right)\right]$, respectively.

Equations (81)–(89) can be written in matrix vector form, given by Eq. (A90); see Appendix. We can apply a similarity transform, given by Eq. (A91) (see Appendix), to yield an equation of motion in the transformed basis vector given by Eq. (A92); see Appendix.

Figure 8b shows numerical simulations of the time evolution of the system that is described by Eqs. (81)–(89) with initial hyperpolarized signal/concentration (see above for a comment on this aspect) ${\mathrm{Pyr}}^{\ast }\left(t\right)=\mathrm{0.001}$ and longitudinal relaxation rate constants $R_{\mathrm{1}}^{\mathrm{P}}=R_{\mathrm{1}}^{\mathrm{L}}=\mathrm{1}/\mathrm{60}$ s⁻¹. The kinetic parameters used for lactate dehydrogenase were as previously published (Zewe and Fromm, 1962; Witney et al., 2011) for the rabbit muscle enzyme. Enzyme concentration was ${\left[\mathrm{E}\right]}_{\mathrm{0}}=\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{9}}$ M and rate constants were $k_{\mathrm{1}}=\mathrm{1.03}\times {\mathrm{10}}^{\mathrm{8}}$ M⁻¹ s⁻¹, $k_{-\mathrm{1}}=\mathrm{549}$ s⁻¹, $k_{\mathrm{2}}=\mathrm{6.72}\times {\mathrm{10}}^{\mathrm{6}}$ M⁻¹ s⁻¹, $k_{-\mathrm{2}}=\mathrm{3.44}\times {\mathrm{10}}^{\mathrm{4}}$ M⁻¹ s⁻¹, k₃=842 s⁻¹, and $k_{-\mathrm{3}}=\mathrm{9.12}\times {\mathrm{10}}^{\mathrm{5}}$ M⁻¹ s⁻¹.

The computed time dependence of polarized and unpolarized pools Pyr^∗(t), Pyr(t), Lac^∗(t) and Lac(t) are shown in Fig. 8b, left column of plots. The time dependences of [E(t)], [E.NAD(t)] and [E.NADH(t)] are shown in Fig. 8(b), middle panel. The time dependences of ${\mathrm{Pyr}}^{\ast }\left(t\right)+\mathrm{Pyr}\left(t\right)=\left[\mathrm{Pyr}\left(t\right)\right]$, ${\mathrm{Lac}}^{\ast }\left(t\right)+\mathrm{Lac}\left(t\right)=\left[\mathrm{Lac}\left(t\right)\right]$, [NAD(t)] and [NADH(t)] are shown in Fig. 8b, right panels. Several interesting features are evident. First, the model predicted the expected time dependences of both hyperpolarized pyruvate Pyr^∗(t) and its conversion to Lac^∗(t). Under the conditions of the simulation, the free enzyme [E(t)] was rapidly depleted to form an equilibrium of [E.NAD(t)] and [E.NADH(t)]. During the reaction with Pyr^∗(t), the equilibrium position of the reaction was altered to give a final equilibrium position that could then be appreciated from the total pools of ${\mathrm{Pyr}}^{\ast }\left(t\right)+\mathrm{Pyr}\left(t\right)=\left[\mathrm{Pyr}\left(t\right)\right]$ and ${\mathrm{Lac}}^{\ast }\left(t\right)+\mathrm{Lac}\left(t\right)=\left[\mathrm{Lac}\left(t\right)\right]$, which predicts a net conversion of [Pyr(t)] to [Lac(t)] of ∼10 %. Also note, from this simulation, that the activity of LDH switches off at t=200 s since the concentration of [NADH(t)] is limiting in this simulation; i.e. it becomes depleted. This does not happen if [NADH(t)] is increased. In a normal cellular context NADH would be regenerated by glyceraldehyde-3-phosphate dehydrogenase during glycolysis.

Finally, we consider real case scenarios that are reported in the literature, i.e. measurement of hyperpolarized [1-¹³C] pyruvate kinetics in living cells (Andersson et al., 2007; Day et al., 2007; Karlsson et al., 2007; Hill et al., 2013a, b; Lin et al., 2014; Pagès et al., 2014; Beloueche-Babari et al., 2017). Figure 8c shows the situation where the LDH expression level is altered, e.g. by the progression of disease (LDH expression is known to be upregulated in more aggressive cancer phenotypes – Albers et al., 2008) or downregulation during therapy (Ward et al., 2010), which can be explored through the value of [E]₀. Figure 8c shows simulations of the Lac^∗(t) signal under the conditions that [E]₀= (i) $\mathrm{0.6}\times {\mathrm{10}}^{-\mathrm{9}}$ M, (ii) $\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{9}}$ M, and (iii) $\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{9}}$ M, while all other parameters remained unchanged, relative to those used for Fig. 8b. It is apparent that increased enzyme expression leads to an increase in the apparent rate of conversion of Pyr^∗(t) to Lac^∗(t) even in the absence of a change in enzyme activity, as seen in real experiments. Another situation that is frequently encountered is the change in the pool size of endogenous lactate, e.g. in response to hypoxia, which can be explored through the parameter Lac(0). Figure 8d shows simulations of the Lac^∗(t) signal under the conditions that Lac(0)= (i) 0 mM, (ii) 20 mM, and (iii) 40 mM, while all other parameters remained unchanged relative to those used to generate Fig. 8b. The model therefore predicts that an increased pool of endogenous unpolarized lactate leads to an increase in the rate of conversion of Pyr^∗(t) to Lac^∗(t), as reported in the literature (Day et al., 2007).

7 Conclusions

We have described an approach to formulating the kinetic master equations that describe the time evolution of hyperpolarized ¹³C NMR signals in reacting (bio)chemical systems, including enzymes with two or more substrates, and various enzyme reaction mechanisms as classified by Cleland. The modelling can be the basis for simulating many pertinent features that are seen in dDNP experiments. Derivation of the Michaelis–Menten equation in the context of dDNP experiments illustrates why formation of a hyperpolarized enzyme–substrate complex does not (typically) cause an appreciable loss of the signal from the substrate or product. It was also able to answer why the concentration of an unlabelled pool of substrate, for example ¹²C lactate, causes an increase in the rate of exchange of the ¹³C-labelled pool and to what extent the equilibrium position of an enzyme-catalysed reaction, for example LDH, is altered upon adding hyperpolarized substrate. The formalism described here should contribute to a fuller mechanistic understanding of the time courses derived from dDNP experiments and will be relevant to ongoing clinical applications using dDNP.

Appendix A

$$\begin{array}{}\text{(A69)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}{\mathrm{S}}^{\ast }\left(t\right)\\ {\mathrm{ES}}^{\ast }\left(t\right)\\ {\mathrm{P}}^{\ast }\left(t\right)\\ \mathrm{S}\left(t\right)\\ \mathrm{ES}\left(t\right)\\ \mathrm{P}\left(t\right)\\ \mathrm{E}\left(t\right)\end{array}\right]=\left[\begin{array}{ccccccc}-k_{\mathrm{1}}\mathrm{E}\left(t\right)-R_{\mathrm{1}}^{\mathrm{S}}& k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ k_{\mathrm{1}}\mathrm{E}\left(t\right)& -k_{-\mathrm{1}}-k_{\mathrm{2}}-R_{\mathrm{1}}^{\mathrm{ES}}& k_{-\mathrm{2}}\mathrm{E}\left(t\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& k_{\mathrm{2}}& -k_{-\mathrm{2}}\mathrm{E}\left(t\right)-R_{\mathrm{1}}^{\mathrm{P}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ R_{\mathrm{1}}^{\mathrm{S}}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}\mathrm{E}\left(t\right)& k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& R_{\mathrm{1}}^{\mathrm{ES}}& \mathrm{0}& k_{\mathrm{1}}\mathrm{E}\left(t\right)& -k_{-\mathrm{1}}-k_{\mathrm{2}}& k_{-\mathrm{2}}\mathrm{E}\left(t\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& R_{\mathrm{1}}^{\mathrm{P}}& \mathrm{0}& k_{\mathrm{2}}& -k_{-\mathrm{2}}\mathrm{E}\left(t\right)& \mathrm{0}\\ -k_{\mathrm{1}}\mathrm{E}\left(t\right)& k_{-\mathrm{1}}+k_{\mathrm{2}}& -k_{-\mathrm{2}}\mathrm{E}\left(t\right)& -k_{\mathrm{1}}\mathrm{E}\left(t\right)& k_{-\mathrm{1}}+k_{\mathrm{2}}& -k_{-\mathrm{2}}\mathrm{E}\left(t\right)& \mathrm{0}\end{array}\right]\\ & \left[\begin{array}{c}{\mathrm{S}}^{\ast }\left(t\right)\\ {\mathrm{ES}}^{\ast }\left(t\right)\\ {\mathrm{P}}^{\ast }\left(t\right)\\ \mathrm{S}\left(t\right)\\ \mathrm{ES}\left(t\right)\\ \mathrm{P}\left(t\right)\\ \mathrm{E}\left(t\right)\end{array}\right]\end{array}\end{array}$$

$$\begin{array}{}\text{(A70)}& U=\left[\begin{array}{ccccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\end{array}$$

$$\begin{array}{}\text{(A71)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}{\mathrm{S}}^{\ast }\left(t\right)\\ {\mathrm{ES}}^{\ast }\left(t\right)\\ {\mathrm{P}}^{\ast }\left(t\right)\\ {\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)\\ {\mathrm{ES}}^{\ast }\left(t\right)+\mathrm{ES}\left(t\right)\\ {\mathrm{P}}^{\ast }\left(t\right)+\mathrm{P}\left(t\right)\\ \mathrm{E}\left(t\right)\end{array}\right]=\left[\begin{array}{ccccccc}-k_{\mathrm{1}}\mathrm{E}\left(t\right)-R_{\mathrm{1}}^{\mathrm{S}}& k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ k_{\mathrm{1}}\mathrm{E}\left(t\right)& -k_{-\mathrm{1}}-k_{\mathrm{2}}-R_{\mathrm{1}}^{\mathrm{ES}}& k_{-\mathrm{2}}\mathrm{E}\left(t\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& k_{\mathrm{2}}& -k_{-\mathrm{2}}\mathrm{E}\left(t\right)-R_{\mathrm{1}}^{\mathrm{P}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}\mathrm{E}\left(t\right)& k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{1}}\mathrm{E}\left(t\right)& -k_{-\mathrm{1}}-k_{\mathrm{2}}& k_{-\mathrm{2}}\mathrm{E}\left(t\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{2}}& -k_{-\mathrm{2}}\mathrm{E}\left(t\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}\mathrm{E}\left(t\right)& k_{-\mathrm{1}}+k_{\mathrm{2}}& -k_{-\mathrm{2}}\mathrm{E}\left(t\right)& \mathrm{0}\end{array}\right]\\ & \left[\begin{array}{c}{\mathrm{S}}^{\ast }\left(t\right)\\ {\mathrm{ES}}^{\ast }\left(t\right)\\ {\mathrm{P}}^{\ast }\left(t\right)\\ {\mathrm{S}}^{\ast }\left(t\right)+\mathrm{S}\left(t\right)\\ {\mathrm{ES}}^{\ast }\left(t\right)+\mathrm{ES}\left(t\right)\\ {\mathrm{P}}^{\ast }\left(t\right)+\mathrm{P}\left(t\right)\\ \mathrm{E}\left(t\right)\end{array}\right]\end{array}\end{array}$$

$$\begin{array}{}\text{(A90)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}{\mathrm{Pyr}}^{\ast }\left(t\right)\\ {\mathrm{Lac}}^{\ast }\left(t\right)\\ \mathrm{Pyr}\left(t\right)\\ \mathrm{Lac}\left(t\right)\\ \mathrm{NADH}\left(t\right)\\ \mathrm{NAD}\left(t\right)\\ \mathrm{E}.\mathrm{NADH}\left(t\right)\\ \mathrm{E}.\mathrm{NAD}\left(t\right)\\ \mathrm{E}\left(t\right)\end{array}\right]=\\ & \left[\begin{array}{ccccccccc}-k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)-R_{\mathrm{1}}^{\mathrm{P}}& k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)& -k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)-R_{\mathrm{1}}^{\mathrm{L}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ R_{\mathrm{1}}^{\mathrm{P}}& \mathrm{0}& -k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)& k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& R_{\mathrm{1}}^{\mathrm{L}}& k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)& -k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}\mathrm{E}\left(t\right)& \mathrm{0}& k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{-\mathrm{3}}\mathrm{E}\left(t\right)& \mathrm{0}& k_{\mathrm{3}}& \mathrm{0}\\ -k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)& k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)& -k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)& -k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)& k_{\mathrm{1}}\mathrm{E}\left(t\right)& k_{-\mathrm{3}}\mathrm{E}\left(t\right)& -k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)& -k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)& k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)& -k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{3}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{1}}\mathrm{E}\left(t\right)& -k_{-\mathrm{3}}\mathrm{E}\left(t\right)& k_{-\mathrm{1}}& k_{\mathrm{3}}& \mathrm{0}\end{array}\right]\\ & \left[\begin{array}{c}{\mathrm{Pyr}}^{\ast }\left(t\right)\\ {\mathrm{Lac}}^{\ast }\left(t\right)\\ \mathrm{Pyr}\left(t\right)\\ \mathrm{Lac}\left(t\right)\\ \mathrm{NADH}\left(t\right)\\ \mathrm{NAD}\left(t\right)\\ \mathrm{E}.\mathrm{NADH}\left(t\right)\\ \mathrm{E}.\mathrm{NAD}\left(t\right)\\ \mathrm{E}\left(t\right)\end{array}\right].\end{array}\end{array}$$

$$\begin{array}{}\text{(A91)}& U=\left[\begin{array}{ccccccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],\end{array}$$

$$\begin{array}{}\text{(A92)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\begin{array}{c}{\mathrm{Pyr}}^{\ast }\left(t\right)\\ {\mathrm{Lac}}^{\ast }\left(t\right)\\ {\mathrm{Pyr}}^{\ast }\left(t\right)+\mathrm{Pyr}\left(t\right)\\ {\mathrm{Lac}}^{\ast }\left(t\right)+\mathrm{Lac}\left(t\right)\\ \mathrm{NADH}\left(t\right)\\ \mathrm{NAD}\left(t\right)\\ \mathrm{E}.\mathrm{NADH}\left(t\right)\\ \mathrm{E}.\mathrm{NAD}\left(t\right)\\ \mathrm{E}\left(t\right)\end{array}\right]=\\ & \left[\begin{array}{ccccccccc}-k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)-R_{\mathrm{1}}^{\mathrm{P}}& k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)& -k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)-R_{\mathrm{1}}^{\mathrm{L}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)& k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)& -k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{1}}\mathrm{E}\left(t\right)& \mathrm{0}& k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{-\mathrm{3}}\mathrm{E}\left(t\right)& \mathrm{0}& k_{\mathrm{3}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)& k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)& k_{\mathrm{1}}\mathrm{E}\left(t\right)& k_{-\mathrm{3}}\mathrm{E}\left(t\right)& -k_{-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& k_{\mathrm{2}}\mathrm{E}.\mathrm{NADH}\left(t\right)& -k_{-\mathrm{2}}\mathrm{E}.\mathrm{NAD}\left(t\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& -k_{\mathrm{3}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& k_{\mathrm{1}}\mathrm{E}\left(t\right)& -k_{-\mathrm{3}}\mathrm{E}\left(t\right)& k_{-\mathrm{1}}& k_{\mathrm{3}}& \mathrm{0}\end{array}\right]\\ & \left[\begin{array}{c}{\mathrm{Pyr}}^{\ast }\left(t\right)\\ {\mathrm{Lac}}^{\ast }\left(t\right)\\ {\mathrm{Pyr}}^{\ast }\left(t\right)+\mathrm{Pyr}\left(t\right)\\ {\mathrm{Lac}}^{\ast }\left(t\right)+\mathrm{Lac}\left(t\right)\\ \mathrm{NADH}\left(t\right)\\ \mathrm{NAD}\left(t\right)\\ \mathrm{E}.\mathrm{NADH}\left(t\right)\\ \mathrm{E}.\mathrm{NAD}\left(t\right)\\ \mathrm{E}\left(t\right)\end{array}\right].\end{array}\end{array}$$