Theoretical Analysis of Discreteness-Induced Transition in Autocatalytic Reaction Dynamics

Transitions in the qualitative behavior of chemical reaction dynamics with a decrease in molecule number have attracted much attention. Here, a method based on a Markov process with a tridiagonal transition matrix is applied to the analysis of this transition in reaction dynamics. The transition to bistability due to the small-number effect and the mean switching time between the bistable states are analytically calculated in agreement with numerical simulations. In addition, a novel transition involving the reversal of the chemical reaction flow is found in the model under an external flow, and also in a three-component model. The generality of this transition and its correspondence to biological phenomena are also discussed.

As the total number of molecules in a system decreases, transitions in the qualitative behavior of autocatalytic chemical reaction dynamics may appear. By adapting a method based on the discrete-time Markov process, we provide here a general analytic tool for the reaction dynamics of systems with a small number of molecules. Bistability induced by the small-number effect in a two-component model is analyzed and excellent agreement with the simulated mean switching time is demonstrated. A novel transition involving the reversal of the chemical reaction flow is also found in a three-component model and quantified analytically using the proposed method. Temporal changes in chemical concentrations are often analyzed by using the rate equation of the reaction kinetics in which a set of deterministic ordinary differential equations is adopted. Fluctuation around the average change in the concentration is neglected by assuming that the total number of molecules is sufficiently large; however, stochasticity in the reaction due to fluctuation in the number of reactants does exist, and it is non-negligible, especially when the number of molecules is small. Such stochasticity can introduce a qualitative change in the behavior of the reaction dynamics. The fluctuation around the average behavior is typically analyzed using the linear noise approximation (LNA) [1], which is represented by the Langevin equation with additive noise or the corresponding Fokker-Planck equation derived from van Kampen's system-size expansion [1,2]. Recently, several theoretical studies have reported "small-number effects" or "discretenessinduced transitions" that lead to a qualitative deviation from the behavior expected by the LNA, due to the small number of components; these effects have been reported in catalytic reaction dynamics [3][4][5][6][7][8][9][10], reaction-diffusion systems [11][12][13], gene regulatory circuits [14], and ecology [15,16].
One remarkable example of a discreteness-induced transition, the emergence of multi-stability, was reported by Togashi and Kaneko [3]. In their investigation of a catalytic chemical reaction system composed of four chemical species, it was revealed that when the total number of molecules is small, the system exhibits temporal switching between quasi-steady states in which only two chemical species are abundant and the others are extinct. The model was then simplified into two components and has been analyzed [4,5,16], but the small-number regime requires further mathematical analysis.
In this Letter, we propose a method based on the discrete-time Markov process for analyzing the smallnumber effect in chemical reactions. After demonstrating the validity of the method in an analysis of the two-component Togashi-Kaneko (2TK) model, we introduce a three-component model that exhibits a novel discreteness-induced transition: the current of the chemical reaction reverses when the total number of molecules becomes small. By applying both the method proposed here and the Fokker-Planck equation, the transition caused by a decrease in the number of molecules is explained with quantitative agreement with numerical simulations.
The 2TK model [4,16] consists of the following four chemical reactions involving the two chemicals A and B: where k, u, and v are the rate constants of the reactions. Note that the total number of molecules, denoted by N , is conserved in this model. Without loss of generality, we can set k = 1 by rescaling the time scale and set u/k → u and v/k → v. By denoting the number of molecules of A and B as i and N − i, respectively, the model can then be defined in the state space i = 0, 1, . . . N . The transition rates from state i to i + 1 and from i to i − 1 are given by where a volume of the system is set to be N . Examples of the time series of i/N for u = v = 0.01, obtained numerically using the Gillespie algorithm [17], are shown in Fig. 1(a). Although the result for large N (N = 10000), shown by the black line, converges to the fixed-point con- centration i/N = v/(u + v) predicted by the rate equation, switching behavior between i/N = 0 and i/N = 1 emerges for small N (N = 50), as shown by the blue line. The emergence of this switching behavior is an example of a discreteness-induced transition [3] or noise-induced bistability [16]. Previous studies [5,16] give the critical value of N for the appearance of switching behavior to be N c ≃ uN for the case v = u. It is especially noteworthy that this model is similar to the two-allele Moran model with mutation in population genetics [18].
Although an analysis based on the Fokker-Planck equation is often adopted for chemical reaction systems, in general, it is not applicable when N is small, which is the case in which we are interested. Thus, we apply here a method based on a discrete-time Markov process without assuming a large N . To illustrate this method, we estimate the mean switching time from state i = 0 to state i = N , since this represents the characteristic time scale of the reaction dynamics, namely, the time required for switching. Similar treatments have been adopted in population genetics [19] (see also [15]). Assuming that no more than one of the reactions in Eq. (1) occurs within a time unit [20], the transition probabilities within the time unit from state i to i + 1, λ i , and from i to i − 1, µ i , are given as To estimate the mean switching time from state i = 0 to state i = N , we set i = 0 as the initial condition and i = N as the final absorbing state. Thus, the chemical reaction system in Eq. (1) is now described by a discretetime Markov process with the transition probabilities λ i and µ i . The time spent in state j starting from state i, denoted by t ij , satisfies where δ ij is the Kronecker delta. Noting that the boundary conditions are µ N = 0 and t N,j = 0 for j = 0, . . . N , the time t 0,j is calculated as where j k=j+1 (µ k /λ k ) = 1. This can be obtained approximately as shows the estimated T 0 , scaled by uT 0 , which agrees well with that obtained from numerical simulations. Note that Biancalani et al. [16] calculated the mean switching time using the Fokker-Planck equation for large N and estimated that for small N in a heuristic manner. Our treatment, in contrast, gives a single expression, Eq. (6), that shows remarkable agreement for all N . It should be emphasized here that the expressions in Eqs. (4)- (6) are not limited to this specific model but are generally applicable to any λ i and µ i . The analogy between the Moran and 2TK models suggests that for small N the reaction dynamics involving an autocatalytic reaction such as A + B → 2A/2B generally give rise to "fixation" to a state in which only a single species dominates while the others are extinct. The outcome of the small-number effect is thus not limited to the emergence of bistability. As an example, we introduce here a three-component model, with autocatalytic reaction motifs, that exhibits a novel transition, i.e., the reversal of the net current of the chemical reactions, induced by the small-number effect.
The three-component model, illustrated in Fig. 2(a), involves the following reactions of the three chemical species A, B and C in a container with a volume N : model is given bẏ where x, y, and z are the concentrations of A, B, and C, respectively. Since the detailed balance condition is not satisfied, a cyclic molecular current can exist, and we define the net reaction current J as J = J a→b + J b→c + J c→a . Here, J a→b is given by Although the rate equation in Eq. (9) predicts a stationary clockwise current J > 0 for k c > 3u, the numerical simulation gives a counterclockwise current J < 0 when the total number of molecules N is small. As shown in Fig. 2(b), (c), and (d), the net current J decreases with decreasing N and becomes negative. Hence, the smallnumber effect leads to a reversal of the net current of the chemical reaction. We analyze this transition by using both the Fokker-Planck equation derived from the master equation for large N and the proposed method for small N .
The master equation of this model is expressed aṡ where T (xy|x ′ y ′ ) is obtained from Eq. (8) by letting x = a/N , y = b/N , and z = c/N = 1 − x − y. Using a Kramers-Moyal expansion to second order, the Fokker-Planck equation can be obtained from the master equation aṡ where Note that · represents the average over the transition rates. Here, the stationary solution of Eq. (12) cannot be obtained analytically, since the transition rates do not satisfy the detailed balance condition. As an analytical estimate, we approximate the stationary solution of Eq. (12) by a Gaussian distribution of mean (m x , m y ) and variance-covariance matrix Σ as P 0 (x, y) ∝ exp −vΣ −1 v/2 , where v = (x − m x , y − m y ). Under this approximation, m x , m y , the variances V x and V y , and the covariance V xy are obtained as m x = m y = 1/3, .
Note that this Gaussian approximation is invalid when V x or V y are large because the estimated probability distribution P st (x, y) then extends beyond the domain x is a necessary condition for the validity of the Gaussian approximation.
The net current estimated from the Gaussian approximation is expressed as The second term in this equation is a result of replacing the probability mass of P 0 (x, y) outside Ω with a probability of states with (x, y) = (1, 0), (0, 1), or (0, 0). As shown in Fig. 2(b), (c), and (d), the current J calculated using Eq. (14) agrees rather well with the numerical results, as long as V x ≤ m 2 x holds. For small N , the variance is much larger, and so the Fokker-Planck equation is no longer valid. The proposed method, however, is applicable. In contrast to the case with the 2TK model, the three-component model has two degrees of freedom, and thus the relationship in Eq. (4) cannot be applied straightforwardly. Nevertheless, we can apply Eq. (4) by focusing only on a relatively short time scale within which the process can be approximated as a single-variable Markov process. In the present model, autocatalytic reactions (e.g., A + C → 2A/2C) tend to push the system into a state in which only a single type of molecule dominates, while the other reactions promote the coexistence of the three molecular species. Therefore, the model can be regarded as a molecule (e.g., A or C) extinction process if the autocatalytic reactions dominate.
Let us assume that only A molecules exist in the initial condition and that the system transits from (a, b, c) = (N, 0, 0) to (a, b, c) = (N − 1, 0, 1) by the reaction A → C at t = 0. If u is sufficiently small, one can assume that one-body reactions (e.g., A → C) do not occur within the time interval 1/N u but that a molecule of A or C becomes extinct within the interval. These assumptions allow us to approximate the behavior of the model by a single-variable Markov process with the following transi-tion probabilities λ i and µ i : where i now represents the number of A molecules. In these equations, both i = 0 and i = N are absorbing states. During this transition process, any B molecules that are generated are counted as C molecules. The time t N −1,j spent in state j before the absorption starting from the initial state i = N −1 is calculated using Eq. (4) with the boundary condition t 0,j = t N,j = 0; for k c = 0, where α = 1 + k c and t N −1,j = (N − j) −1 for k c = 0. The total numbers of molecules that flow clockwise and counterclockwise within the interval 1/N u are then calculated as This estimate is valid only when the time until absorption T N −1 = N −1 j=1 t N −1,j is smaller than 1/N u. The current J, shown by the blue lines in Fig. 2(b), (c), and (d), agrees rather well with the numerical results when T N −1 < 1/N u holds. Figure 3 shows phase diagrams of the regime for clockwise or counterclockwise currents. The phase boundary is calculated by the criterion J = 0. For small N , it is estimated from the solution of the equation: 2α 1−N + k c N − 2α = 0 [21]. For large N , the boundary is estimated from Eq. (14). Both the estimated boundaries agree well with the simulation. Figure 3 (c) and (d) clearly show that the phase boundary is independent of u as long as T N −1 < 1/N u holds.
In this Letter, we have presented a general method based on a discrete-time Markov process for analyzing small-number effects in chemical reactions. The method has been applied to analyze the emergence of bistability in the 2TK model, and also to analyze the reversal of the net reaction current in the three-component model. The three-component model here consists of homogeneous cyclic reactions, but extensions to chain-like reactions and systems with inhomogeneous reaction coefficients [6] are fairly straightforward. In the case of chain-like reactions, the dominant molecule in the steady state may change depending on N because the direction of the reaction current may change. Both our analytical methodology and the insight provided by these findings will assist in understanding several of the small-number effects reported recently in real cellular systems [14,22]. This work was partially supported by a Grant-in-Aid for Scientific Research (No. 21120004) on Innovative Areas "Neural creativity for communication" (No. 4103) and the Platform for Dynamic Approaches to Living System from MEXT, Japan. We acknowledge helpful discussions with Shuji Ishihara.