Continuous-time Markov-chain models for reaction systems: fast and slow processes

It is sometimes of interest to identify the fast and slow processes in a reaction system. We present here an approach to this problem which is based on a simple stochastic model, a continuous-time Markov chain on a small number of states. We show how it is possible to use such a stochastic model to find and plot the time-courses of concentrations, and to find simple short-time and long-time approximations to these time-courses; that is, to separate the fast and the slow processes. The most significant computation involved is the exponentiation of many small matrices, which is easily accomplished in the computing environment R.


Introduction
In the modelling of chemical and other reaction systems, there are essentially two alternatives: deterministic models and stochastic (i.e. probabilistic) ones. The former type enables one to draw conclusions about (e.g.) concentrations from a system of ordinary differential equations, the rate equations, the latter typically seeks to answer similar questions by finding transition probabilities in a continuous-time Markov chain representing the reactions.
To many chemists and other physical scientists, deterministic models and their associated differential equations are probably the more familiar approach, but that is not necessarily the simpler route to a useful conclusion, merely the more familiar of 1 3 the two. Similarly, those of a statistical background would tend to choose stochastic models. The books of Érdi and Tóth [4] and Érdi and Lente [3] cover (respectively) both deterministic and stochastic models, and stochastic models only.
One sometimes sees authors writing as if the system of rate equations IS the system of chemical reactions, as opposed to one very useful representation thereof; we prefer to distinguish the two systems. For example, there are literature references to Robertson [12] which cite that work as the source of the 'Robertson reaction'. There are no reactions at all in [12], but there is discussed in that work a system of differential equations (DEs) that could very well arise as the rate equations implied by a system of chemical reactions. Furthermore, it is in our view useful to have at one's disposal both deterministic and stochastic modelling. There may be insights readily available from one of the techniques that are not as readily available from the other. See, for instance, MacDonald [10, pp. 3,4], in which continuous-time Markov chains are used (instead of the rate equations) as a simple route to clear-cut conclusions about a unimolecular triangle reaction system.
Alvarez-Ramirez et al. [1] use the rate equations, the singular value decomposition (SVD) of a certain matrix arising in those equations, and singular perturbation theory in order to separate the fast and slow processes in a reaction system. We present here an alternative formulation of their worked examples which demonstrates how continuous-time Markov chains (CTMCs) can instead be used to study the time-courses of the concentrations of the species involved, and thereby draw conclusions about the fast and slow processes in the system. It is to some extent possible to derive implied fast or slow manifolds, but examination of these time-courses is also informative, and apparently easier to carry out.
Our approach is essentially that described by MacDonald [10]; states are defined, the generator matrix of the (time-homogeneous) CTMC representing the system is identified, and the generator is then exponentiated in order to find the transition probabilities, which can deduce the vector of concentrations at a given time from the initial concentrations. The time-courses of the concentrations can thereby be computed and plotted.

The reaction system
The first example of Alvarez-Ramirez et al. [1] relates to the kinetics of three-lump hydrocracking (HDC). The system, consisting of three monomolecular reactions, is as follows: and can be depicted as in Fig. 1. The symbol R denotes heavy oil residue, and L and G are the liquid and the gas products. (In passing: the use of k 3 , not k 2 , in the second reaction above is intentional.) As noted by Alvarez-Ramirez et al., the system of ordinary differential equations (ODEs) representing the reaction system is: The linearity and triangular structure of the Eq. 1 make it easy to solve them explicitly, 1 but we take a slightly different route.
In view of the explicit solution, the approximations based on the SVD and singular perturbation may seem redundant, but the implications of such approximations can if necessary be compared to those of the exact solution. In the example presented by Alvarez-Ramirez et al., the values of the rate constants (per hour) are: k 1 = 0.2787 , k 2 = 0.0361 , k 3 = 0.0101 ; these values arise in the work of Soto-Azuara et al. [13].

The Markov-chain model
With the states being R, L and G, the generator of the associated CTMC is (We use the convention that the row sums of such a matrix are zero.) The last row reflects the fact that G is an absorbing state. Since Q is upper-triangular, its eigenvalues are −(k 1 + k 2 ) , −k 3 and 0. Apart from certain obvious special cases (which we exclude here), these eigenvalues are distinct and the generator therefore diagonalizable. A matrix of right (i.e. column) eigenvectors, with the same ordering, is the inverse of which can be found explicitly: One can therefore find the required transition probability matrix (TPM) P(t) = exp(tQ) explicitly as where D t is the diagonal matrix with diagonal (e −(k 1 +k 2 )t , e −k 3 t , 1) . In any case, it follows that each transition probability p ij (t) is of the form where a, b and c are independent of t. For instance, the second of which peaks at time t = The (relative) concentrations at time t, adding to 1 and corresponding to the three states, are supplied by the row vector where P(t) is given either by the explicit expression derived above, or numerically by P(t) = exp(tQ), the matrix exponential being evaluated by (e.g.) the excellent R routine expm [5]. (If a matrix is diagonalizable, it is easily exponentiated; expm does not assume diagonalizability.) For ease of exposition we shall assume here that u(0) = (1, 0, 0) , but that assumption is inessential. The concentrations at time t are then with the transition probabilities p 1j (t) as given by (2). Differences or sums of concentrations might well be of interest, and are of course easily deduced.

Fast and slow processes
Since we have simple explicit expressions for the concentrations, we can directly investigate their approximate behaviour for small t and for large t, and thereby possibly find fast and slow 'manifolds'. Firstly, note that, if (k 1 + k 2 )t and k 3 t are sufficiently small, and Hence This identifies a fast manifold which is almost exactly in agreement with that given by Eq. 21 of Alvarez-Ramirez et al.: It is obviously less in agreement with their Eq. 33: If one makes the assumption (only) that k 3 = 0 , then it follows that (exactly) so the fast manifold discussed above is an exact consequence of the assumption that k 3 = 0.
Secondly, note that, if (k 1 + k 2 )t is sufficiently large, u 1 (t) ≈ 0 and so u 2 (t) + u 3 (t) ≈ 1 , which identifies a slow manifold. But the search for fast and slow manifolds does not seem to add much information if it is easily possible to plot u(t) for a given generator Q and initial distribution u(0).
If (as in this example) k 3 is less than k 2 and very much less than k 1 , it is clear that conversion to gas is slow, as seen in Fig. 2, which can be compared with Fig. 1 of Alvarez-Ramirez et al.

Generator lacking distinct eigenvalues
It has been suggested that we should not confine ourselves here to the case of the generator Q having distinct eigenvalues. As pointed out by Hirsch et al. [6, p. 100], 'most' matrices have distinct eigenvalues. Hence 'most' (square) matrices are diagonalizable. More precisely, the set of real n × n matrices that have n distinct eigenvalues is (open and) dense in the set of all real n × n matrices [6, p. 101]. The case of Q being deficient in that respect is therefore of theoretical rather than practical significance. To proceed with a theoretical analysis, one would need to write the generator in Jordan form or Schur triangular form; for the latter, see Meyer [11, p. 524]. But there would seem to be little to be gained by doing so. If such a deficient generator Q were to turn out to be of interest, it could of course be nondiagonalizable, but that would not prevent one from exponentiating tQ with appropriate software, e.g. expm in R.

Kinetics of four-lump hydrocracking
In addition to three-lump hydrocracking, Alvarez-Ramirez et al. also discuss more briefly the four-lump HDC model of Soto-Azuara et al [13]. Here the four states are (in order) R, HLP (heavy liquid product), LLP (light liquid product) and G, and the generator of our CTMC model is Its eigenvalues are −(k 1 + k 2 + k 3 ) , −(k 4 + k 5 ) , −k 6 and 0. Again, apart from certain special cases the eigenvalues are distinct and the generator therefore diagonalizable. But in the example considered, one of those special cases does indeed arise, since k 4 = k 5 = 0 . Zero is then a repeated eigenvalue, and states 2 (HLP) and 4 (G) are both absorbing states. Alvarez-Ramirez et al. describe the resulting model as 'physically meaningless', but it seems worthwhile to explore the consequences, for the CTMC, of the assumption that k 4 = k 5 = 0.
The full set of rate constants considered by Soto-Azuara et al. and Alvarez-Ramirez et al. is: k = (0.2323, 0.0487, 0.0336, 0, 0, 0.1244) . The two absorbing states give rise to an infinity of equilibrium distributions (all convex linear combinations of the two row vectors (0, 1, 0, 0) and (0, 0, 0, 1)). But from the structure of the resulting network it seems plausible that, for a given initial distribution u(0) , there is a unique limiting distribution.
That conjecture can be proved as follows. In spite of the repeated eigenvalue (zero) the generator Q can be diagonalized: two linearly independent right eigenvectors, both corresponding to the eigenvalue zero, are (k 1 ∕(k 1 + k 2 + k 3 ), 1, 0, 0) T and (1, 1, 1, 1) T . More simply, perhaps: two linearly independent left eigenvectors, corresponding to the eigenvalue zero, are (0, 1, 0, 0) and (0, 0, 0, 1). Hence, for some nonsingular matrix U, provided that k 1 + k 2 + k 3 and k 6 are strictly positive. It then follows that The model may indeed be physically meaningless, in that it precludes conversion of HLP to LLP or G, but its properties are not difficult to establish. Fig. 3 displays the time-courses of the concentrations of the four species corresponding to k as given above and initial concentrations u(0) = (1, 0, 0, 0).

The reaction system
The other example of Alvarez-Ramirez et al. is the following system: which can also be displayed as follows: Alvarez-Ramirez et al. give the four ODEs for the four concentrations In the second one (their Eq. 44b) there appears to be a typographical error; the term −k 2 C 3 should apparently be +k 2 C 3 . Note in particular that three of these four rate equations (Eqs. 44a-44c) involve a nonlinear term; the system of rate equations is nonlinear.
As Alvarez-Ramirez et al. note, the traditional Michaelis-Menten enzymatic reaction scheme is an example of such a system, with the species A 1 -A 4 being (respectively) the substrate, enzyme, intermediate complex ('adduct') and product. (But advances in single-molecule experimental technology suggest that the Michaelis-Menten scheme may, at least in some cases, be an oversimplification of such reactions: see Kou et al. [8] and Kou [7]. Kou [7] describes a very detailed CTMC model for the case of the enzyme -galactosidase.)

The Markov-chain model
In the case of this system it is less obvious than in the preceding examples how one should define the states of a CTMC. We explore here a simple three-state CTMC involving the four species, and investigate what useful conclusions can be drawn from it, including conclusions relating to fast and slow processes. We define the following three states: -state 1, (1,1,0,0), that is, one molecule of A 1 and one of A 2 ; -state 2, (0,0,1,0), one molecule of A 3 ; -state 3, (0,1,0,1), one molecule each of A 2 and A 4 .
Unless one is considering only first-order reactions, one cannot assume that the transition rates of the CTMC are in the same units as are the rate constants of the deterministic model, but for notational convenience we shall also use the symbols k i for the three positive transition rates in the CTMC. The generator is then We assume (unless it is otherwise indicated) that all the k i s are strictly positive. The model is a random process that is in exactly one of the three states at any given time, although once it reaches state 3 it remains there. It is perhaps worth stressing that the three probabilities u j (t) add to 1, unlike the four concentrations C i (t) of Alvarez-Ramirez et al.: note (e.g.) the illustrative values C 1,0 = 3 and C 2,0 = 2 that they use. In our model, the relation between the three probabilities u j (t) and the four concentrations [A i ] is as follows: The generator used here is not triangular; but it is tridiagonal, which implies detailed balance. Its nonzero eigenvalues are which are real (and negative). Denote these eigenvalues by 1 and 2 , with 1 < 2 .
The following is a matrix of right eigenvectors corresponding to (respectively) 1 , 2 and zero: its inverse is: (The last row is the equilibrium distribution (0,0,1), a left eigenvector corresponding to eigenvalue zero.) The TPM P(t) can then be found explicitly by diagonalization. For instance: (Notice that p 11 (t) + p 12 (t) = ( 2 e 1 t − 1 e 2 t )∕( 2 − 1 ) .) If we assume that u(0) = (1, 0, 0) , we therefore have explicit expressions for u j (t) = p 1j (t) (j=1, 2, 3). By one means or another, therefore, it is possible to compute the probabilities u j (t) corresponding to states 1, 2 and 3, for a given time t and given initial distribution For case II ( k 3 dominant), Fig. 5 displays the corresponding plot. There is a marked difference between the two plots as regards [A 3 ] , the concentration of the intermediate complex (the dashed red line). In case I, the concentration of A 3 peaks at height about 0.6 before descending to its limiting value, zero. In case II, A 3 is (as one would expect) very quickly converted to A 2 + A 4 , and its concentration never rises above very low values.
By plotting u 2 (t) (the concentration of A 3 , the intermediate complex) against u 1 (t) (the concentration of A 1 , the substrate) it is also possible to produce plots similar to the phase plots of Alvarez-Ramirez et al., their Figs. 3 and 4. The initial conditions differ, however; our initial condition is u(0) = (1, 0, 0) , i.e. (with probability 1) one molecule each of A 1 (substrate) and of A 2 (enzyme). Our

Fast and slow processes
As regards identifying fast and slow processes, we explore first the consequences of the simple assumption that k 3 = 0 (in order to identify the fast process). If k 3 = 0 is assumed, it follows that 1 = −(k 1 + k 2 ) , 2 = 0 and, rather obviously, that p 11 (t) + p 12 (t) = 1 . With the initial condition u(0) = (1, 0, 0) , this produces the relation u 1 (t) + u 2 (t) = 1 . If k 3 is 'small' (as in case I), this assumption should produce a reasonable approximation for small t, i.e. the fast dynamics. An alternative route to the same conclusion proceeds as follows. Equation 4 implies that hence that p 11 (t) + p 12 (t) ≈ 1 for small t. This approximation could be improved by not ignoring the (rather simple) term of order t 2 .
In order to approximate the slow dynamics, consider the ratio p 11 (t)∕p 12 (t) , and the limit thereof as t tends to infinity. From Eq. 4 we have (Recall that 1 < 2 < 0 . It can also be shown that 1 < −k 1 , so this limit is indeed positive, as expected.) With the usual assumption that u(0) = (1, 0, 0) , this implies that u 2 (t)∕u 1 (t) tends to k 1 −(k 1 + 1 ) , hence the approximation, for large t, that u 2 (t)∕u 1 (t) ≈ k 1 −(k 1 + 1 ) . This approximation was used in Figs. 6 and 7, and appears there as a dashed and dotted red line.

A variation
A variation that may be of interest is to assume that the concentration of A 1 (in the Michaelis-Menten scheme, the substrate) is constant. That is the context considered by Kou [7]. The system under discussion here then becomes: and the resulting three ODEs can be written down. The corresponding generator is The TPM P(t) , and anything based on it, e.g. probabilities of the kind given by Eq. 4, then follow.

Nonlinearity
It might seem that the use of CTMCs as models merely provides a different route, in different notation and terminology, to conclusions that are also available from the theory of linear dynamical systems. In our opinion that would not be correct. In fact, one of the examples already considered (the catalytic reaction system) involves significant nonlinearity. In that case, three of the four rate equations include a nonlinear term. The formulation of the system as a CTMC enabled us to proceed with an explicit analysis in spite of this nonlinearity.
We provide here two further examples of cases in which the rate equations are nonlinear, yet the technique of using CTMCs applies without any essential alteration.
Consider the very simple reaction A 1 + A 2 ⟶ A 3 , with rate proportional to the product of the concentrations [A 1 ] and [A 2 ] . That is, the reaction is of first order in both A 1 and A 2 . The resulting rate equations are all nonlinear. With C i denoting [A i ] , and the prime symbol denoting differentiation with respect to time, these are, for some k > 0 : A rather more complex example would be the 'Robertson reaction': (This system of reactions is consistent with the system of three nonlinear ODEs discussed by Robertson [12], and does appear in some of the literature, e.g. Amat et al. [2].) If one adds the nonreacting species A 3 to the first reaction, the system can be written as Now take A 1 + A 3 , A 2 + A 3 and 2A 2 (in that order) to be the states of a CTMC with generator of the form That is, state 1 is one molecule of A 1 plus one of A 3 , state 2 is one each of A 2 and A 3 , and state 3 is two of A 2 . The generator is tridiagonal, which implies detailed balance. The probability distribution over the three states at time t is given, as usual, by the row vector u(t) = u(0) exp(tQ) . The concentrations C j ≡ [A j ] are given in terms of u(t) by the following (invertible) relationship:

Discussion and conclusion
Both the three-lump HDC model and (as formulated here) the catalytic reaction are special cases of the general three-species monomolecular reaction system in which some of the rates q ij are zero; see Fig. 8. But regardless of whether some rates are

(t) and then use u(t) = u(0)P(t).
A phenomenon that can cause difficulty, in both stochastic and deterministic models for reaction systems, is the coexistence of very fast and very slow reactions. In the deterministic case, the system of differential equations is then 'stiff', and special-purpose numerical software designed for stiff systems can solve the problem. In the stochastic case, one can use the technique known as uniformization (see e.g. van Dijk et al. [14]). The basis of that technique is that the transition probability matrix P(t) = exp(tQ) (surprisingly?) is exactly if ≥ max i |q ii | and the (discrete-time) transition probability matrix R is defined by R = I + Q∕ . For this result see Lange [9, p. 207]. The TPM P(t) is then approximated by taking sufficient terms in the above infinite series.
It is interesting to compare the structures of the two models for the catalytic reaction. One consists of four ODEs (three being nonlinear) in the four concentrations, the other is a (time-homogeneous) CTMC on three states. Both seek to represent the evolution of the process over time, and in this endeavour they use the same information: the structure of the reaction system, and the three quantities k i , although these will not all be in the same units for the two models. The CTMC is more restricted in the initial conditions that it can represent.
Some of what we have described here can certainly be done by using the standard 'dynamical system' approach to the rate equations, although (unfortunately) there are differences in commonly used notation which lead to the conclusions of one approach being a transposed version of those of the other. The use of CTMCs as models for systems of reactions is of course not new, and it may be a matter of taste and background which approach one prefers. But it would seem that the two approaches can answer similar questions, and the use of a CTMC seems simple and direct, at least as regards the derivation of fast-and slow-process approximations; singular perturbation theory is not required. There would be no great difficulty in using a CTMC in a case that required more states than the three or four needed in the examples discussed here. Furthermore, the CTMC technique can avoid the complications which result from nonlinearity in the rate equations.
Acknowledgements The editor-in-chief and the reviewers are thanked for their helpful comments and suggestions.
Funding Open access funding provided by University of Cape Town.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission