An Analytical Framework for Studying Small-Number Effects in Catalytic Reaction Networks: A Probability Generating Function Approach to Chemical Master Equations

Cell activities primarily depend on chemical reactions, especially those mediated by enzymes, and this has led to these activities being modeled as catalytic reaction networks. Although deterministic ordinary differential equations of concentrations (rate equations) have been widely used for modeling purposes in the field of systems biology, it has been pointed out that these catalytic reaction networks may behave in a way that is qualitatively different from such deterministic representation when the number of molecules for certain chemical species in the system is small. Apart from this, representing these phenomena by simple binary (on/off) systems that omit the quantities would also not be feasible. As recent experiments have revealed the existence of rare chemical species in cells, the importance of being able to model potential small-number phenomena is being recognized. However, most preceding studies were based on numerical simulations, and theoretical frameworks to analyze these phenomena have not been sufficiently developed. Motivated by the small-number issue, this work aimed to develop an analytical framework for the chemical master equation describing the distributional behavior of catalytic reaction networks. For simplicity, we considered networks consisting of two-body catalytic reactions. We used the probability generating function method to obtain the steady-state solutions of the chemical master equation without specifying the parameters. We obtained the time evolution equations of the first- and second-order moments of concentrations, and the steady-state analytical solution of the chemical master equation under certain conditions. These results led to the rank conservation law, the connecting state to the winner-takes-all state, and analysis of 2-molecules M-species systems. A possible interpretation of the theoretical conclusion for actual biochemical pathways is also discussed.


INTRODUCTION
Biochemical systems consist of a variety of chemicals, including proteins, nucleic acids, and also small metabolites. Enzymatic reactions, which play an important role in catalyzing many biological reactions, are particularly important to maintain the structure and activity of these systems. Hence, biochemical systems are often modeled as catalytic reaction networks.
These networks are typically analyzed by using deterministic ordinary differential equations with respect to concentrations of chemical species, so-called reaction rate equations [or partial differential equations (reaction-diffusion equations) for spatially distributed, non-uniform cases]; i.e., the concentrations are represented by continuous variables. However, because each chemical actually consists of molecules, the concentration of each species should be a discrete variable. The effects of such discreteness as well as finite-size fluctuations in stochastic reactions become non-negligible if the number of molecules in the system is small. In theory, situations such as these can result in phenomena that cannot be described by rate equations (as well as those equations with additional noise; Kaneko, 2001, 2007;Kaneko, 2007, 2009).
In contrast, gene regulations are often modeled as a combination of binary (on/off) switches, typically represented by Boolean network models (Kauffman, 1969). However, this approach does not consider the quantities of chemicals such as DNA, mRNA, and proteins. Even for the seemingly digital expression of genes, a stoichiometry involving DNA cannot be ignored, as seen in X-chromosome inactivation (2 to 1) and the trisomy syndrome (2 to 3). Thus, the use of a binary (on/off) representation would also be inappropriate. We therefore need to consider the number of molecules.
Recent experiments have shown the existence in the cell of proteins that only consist of a few molecules each (Taniguchi et al., 2010), and these potential biological phenomena (sometimes referred to as the small-number issue) have been gradually recognized by biologists. General theories that would enable predictions in small number situations would be helpful to every biological scientist seeking to understand biochemical processes on any level. Of course, the distributional behavior of such discrete stochastic systems is described by the chemical master equation. However, it is generally difficult to obtain its solution; hence, most efforts have been devoted to the development of approximation and simulation methods and their application (Gillespie, 1976(Gillespie, , 1977Munsky and Khammash, 2006;Lee et al., 2009;Kim and Lee, 2012).
The effects of small numbers in particular systems, e.g., small autocatalytic systems have been mathematically analyzed (Ohkubo et al., 2008;Biancalani et al., 2012Biancalani et al., , 2014Houchmandzadeh and Vallade, 2015;Saito and Kaneko, 2015). In addition, a parameter representing the degree of discreteness in the number of molecules has been introduced (Haruna, 2010(Haruna, , 2015. In this work, we pursue analytical frameworks for studying the effects of small numbers in catalytic reaction networks by following an approach that is as general as possible. Our aim is to obtain the steady-state analytical solution for the chemical master equation without specifying parameters, rather than developing numerical schemes to solve it as was previously done (Kim and Lee, 2012). Furthermore, we try to theoretically describe the long-term behavior of the system by only using information about relationships between elements, thereby implying that we aim to produce results that could be applied to studies in any field.
Our framework provides good operability because our formulas have a specific and satisfyingly simple form, and enables us to obtain the steady state for a wide class of catalytic reaction networks because our framework never specifies any parameters for the networks. We use a probability generating function approach. The probability generating function approach to stochastic chemical kinetics itself has been proposed a long time ago (e.g., Krieger and Gans, 1960;McQuarrie, 1967), after which it spread to biological stochastic kinetics studies (e.g., stochastic gene expression Thattai and van Oudenaarden, 2001;Shahrezaei and Swain, 2008); thus, physicists as well as chemists and mathematical biologists are familiar with the approach. Our main contributions in this paper relate to the efficient usage of the probability generating function (as in Gadgil et al., 2005 regarding first-order reaction networks). Therefore, our approach for obtaining the steady state is understandably easier than that followed in a previous study (Anderson et al., 2010) while, as a consequence, our results are consistent with that study (see Theorem 4.1 and 4.2 in Anderson et al., 2010). Furthermore, since our method uses a procedure based on analytical calculations, it can be easily converted to a computer algorithm.
The present paper is organized as follows. In Section 2.1, we define the catalytic reaction network considered in this paper. The chemical master equation (CME) is provided in Section 2.2. We introduce the probability generating function (PGF) and derive the generating function equation (GFE) in Section 2.3. In subsequent Sections (2.3.1-2.3.3), we show that the GFE introduces the time evolution equation of the first-order and second-order moments of concentrations (we refer to the firstorder moment time evolution equation as the pre-rate equation, PRE), and the second-order moment expression of time-averaged concentrations (SME). Section 2.4 is devoted to obtaining the steady-state solutions of the GFE. To simplify the GFE, we neglect the non-catalytic reactions considered as perturbation for the catalytic reaction network if the system is "entirely ergodic." In Section 2.4.2.3 as the main result, we obtain the probability generating function without winner-takes-all states (PGFwoWTAS) including the solutions of the corresponding rate equation. In Section 3, we describe applications of these results: the rank conservation law, the connecting state to the winner-takes-all state, analysis of 2-molecules M-species systems, and non-autocatalyzation of autocatalytic reaction networks. The prospects of our theory in terms of the small-number issue are briefly discussed in Section 4.

Catalytic Reaction Networks
Consider an abstract catalytic reaction network consisting of M chemical species and N molecules in a well-stirred reactor of volume V, as in Awazu and Kaneko (2007). Each species is labeled by each of the integers between 1 and M. The total number of molecules N is always conserved in reaction processes. This chemical reaction system involves both catalytic reactions and non-catalytic reactions, to prevent catalytic reactions from stopping: (i) Catalytic reactions (two-body catalysis): where the species i, j, k ∈ [1, M] represent a substrate, a catalyst, and a product. The reaction rate constant is represented by R ijk > 0. If this catalytic reaction does not exist, we specify R ijk = 0. Therefore, the catalytic reaction networks are determined by R ijk . In this paper, we impose the following conditions for the catalytic reaction networks; One product against a substrate and a catalyst.
(ii) Non-catalytic reactions (one-body reactions): This reaction exists for all combinations between each species, but a product j ∈ [1, M] is uniformly-randomly selected from all species 1 to M. The rate constant is set ε > 0 in common, where ε is very small, i.e., ε ≪ min{R ijk > 0}.
The state of this catalytic reaction network is specified by the combination of M natural numbers n = (n 1 , n 2 , · · · , n M ), where n i ∈ [0, N] is the number of molecules of the ith-species.
For later convenience, we introduce the following notations: the state space of the catalytic reaction network W M,N (abbr. W) is represented by which consists of M N : = (N+M−1)! (M−1)!N! points; a collection consisting of (one species) winner-takes-all states I M,N (abbr. I) is represented by which consists of M points. Of course, the winners-take-all states of more than one species can be considered, but we focus on the winner-takes-all states of one species by supposing the system satisfies a certain condition, entire ergodicity (see Section 2.4.2.4).
In the present paper, we are interested in the N-dependence of the concentration of each species x i = n i /V. We basically consider a situation in which the total density of molecules ρ = N/V is conserved, even if the total number of molecules N is changed.

The Chemical Master Equation
The rate constant R ijk in the catalytic reaction is defined as the number of reactions per unit volume, unit concentration, and unit time. Therefore, the number of reactions per unit time in the catalytic reaction i j −→ k, such that the concentrations are x i and x j , is where ρ = N/V is the total density of molecules in the vessel. On the other hand, the number of reactions per unit time in the non-catalytic reaction i −→ j with probability 1/M is Then, the time-evolution of the probability P(n, t), with which the system is in the state n at time t, obeys the chemical master equation (CME): where E ±m i are step operators, i.e., E ±m i f (n 1 , · · · , n i , · · · , n M ): = f (n 1 , · · · , n i ± m, · · · , n M ).
Of course, the state space on which the probability P(·, t) is supported, is the previously described W M,N .

Probability Generating Function Method
The probability generating function (PGF) is useful to analyze the CME: Note that the following expressions are translated to differential forms of the PGF: Frontiers in Physiology | www.frontiersin.org Therefore, rewriting the CME to the equation governing the PGF enables the generating function equation (GFE) to be obtained: The GFE consists of continuous variables, unlike the CME, which consists of discrete variables.
Once we obtain the PGF φ as a solution to the GFE, we can derive all statistics of the catalytic reaction network; for example, the ensemble averages (first-order moments) and second-order moments become and the marginal distributions become where p i (n, t) is defined by 2.3.1. The Pre-rate Equation Rate equations are differential equations for the concentrations x i of chemical species i = 1, · · · , M when N → ∞. We use the GFE to derive a formula, which is reminiscent of the rate equation. Differentiating both sides of Equation (11) by z i , substituting z = 1, and writing x i = n i /V, leads to the following pre-rate equation (PRE); Note that, as we are considering a system of which the total number of molecules N is conserved, the concentration conservation law (CCL) must be satisfied; where the PRE (Equation 15) is consistent with the CCL (Equation 16). If the independence x i x j = x i x j (i = j) holds, the PRE (Equation 15) becomes the same expression as the rate equation. On the other hand, the PRE does not explicitly include extensive variables, such as the total number of molecules N. Therefore, the formula (Equation 15) always holds for an arbitrary molecular number N unless the total molecular density ρ is changed.

Second-Order Moment Expression for Time-Averaged Concentrations
We suppose the following ergodicity to replace ensembleaverages with time-averages; n i = n i * i.e., where P * (n) is the steady-state solution of the CME corresponding to almost all initial conditions. Considering the steady state in the PRE [that is, taking t → ∞ on both sides of Equation (15)], one obtains the following second-order moment expression (SME) for time-averaged concentrations: The second term on the right hand side of Equation (18) represents the difference from the uniform concentration ρ/M. If the independence x i x j = x i x j holds, Equation (18) becomes an equation that only includes the concentrations x i . Therefore, in combination with the time-averaged version of the CCL (Equation 16), the concentrations x i can be determined without including unknown quantities. However, the actual concentrations x i depend on their second-order moments x i x j .

Time-Evolution of Second-Order Moments
Determining the concentrations x i from the SME (Equation 18) requires us to know their second-order moments x i x j .
Differentiating both sides of the GFE (Equation 11) by z l and z m (l = m), substituting z = 1, and writing x i = n i /V, then, after simplification, the following time-evolution equation of second-order moments (TESM) is obtained: On the other hand, twice differentiating both sides of the GFE (Equation 11) by z l , substituting z = 1, and writing x i = n i /V, then after some simplification, the following alternative TESM is obtained: The expression of TESMs (Equation 19) consists of M(M + 1)/2 equations containing the third-order moments x i x j x k (t). Therefore, the TESMs are not effective unless the systems are restricted such that the third-order moments vanish (e.g., 2-molecules systems). If N = 2, because ∂ z i ∂ z j ∂ z k φ = 0 (φ is a second-order polynomial), one can calculate the 2-molecules version of the TESMs (2mTESMs); Treatment of the 2mTESMs (Equation 20) to demonstrate their effectiveness appears in a later section.

Steady-State Solutions of GFE
If the GFE (Equation 11) can be solved, this would enable us to obtain all the statistics of the catalytic reaction networks. However, it is generally difficult to solve. Here, we focus on the steady-state solutions of the GFE and consider the case that the ε-term in the GFE can be ignored. Through the following discussion, we see that the approximation is effective only if the system is ergodic.

The Case of Non-catalytic Reactions Only
First, we consider the steady-state solutions of noncatalytic reactions only as an introduction. The PGF φ nc * (z), corresponding to the steady state in the case of non-catalytic reactions only, should satisfy the following equation: By exchanging the subscripts i and j in the second term, one obtains Because the coefficients of variables z i must be zero, φ nc * must satisfy the following equations: equivalently, Equation (24) implies that all of ∂ z i φ nc * (i = 1, · · · , M) are equal to each other, i.e., Considering that the PGF φ nc * (z) is an Nth order polynomial of z i and must satisfy the condition φ nc * (1) = 1 by definition, the following solution of Equation (25) can be found: Therefore, we obtain the following stationary distribution in the case of non-catalytic reactions only: where N n 1 ,n 2 ,··· ,n M : = N! n 1 !n 2 !···n M ! are multinomial coefficients. Furthermore, Equation (13) can also be used to derive the marginal distribution of the i-th species: which is the binomial distribution with parameters N and 1/M. If we suppose the ergodicity n i = n i , the following statistics can be calculated from Equation (12): where

The Case of Catalytic Reactions Only
Next, we consider the steady-state solutions of catalytic reactions only, assuming that the ε-term in the GFE (Equation 11) can be ignored. The steady-state solutions are assumed to have a form similar to Equation (26), including undetermined coefficients (λ i ) deriving from the network structure (R ijk ).

A condition for finding the steady state: λ-condition
The PGF φ c * (z), corresponding to the steady state in the case of catalytic reactions only, should satisfy the following equation; Construction of the steady-state solution requires us to find a particular solution of Equation (30) as bases of linear space consisting of N-th order polynomials. Accordingly, let us assume the following extended form of Equation (26) by introducing parameters {λ i } M i=1 , which eventually correspond to the concentrations (per total density ρ) of each species in the continuous limit N → ∞ as shown in Section 2.4.2.3: because φ * (1) = 1. Substituting Equation (31) into Equation (30) and setting the coefficients of variables z i z j as zero, gives the following condition for {λ i λ j }: The λ-condition (Equation 33) represents M(M − 1)/2 homogeneous equations for λ i λ j ; therefore, λ i can be calculated by combining the λ-condition (Equation 33) with Equation (32). The λ-condition has trivial solutions: These solutions represent the states for which the i-th species take all molecules (the winner-takes-all state). On the other hand, there is also a non-trivial solution. In the following paragraph, we treat three-species systems (M = 3) to demonstrate that nontrivial solutions do exist. The following procedures can easily be extended to those for arbitrary M species systems.
Therefore, the non-trivial solution of Equation (35), i.e., the eigenvector corresponding to the eigenvalue 0 of the matrix A, certainly exists and it has the following expression: The proportionality constant (> 0) can be determined by the condition (Equation 32), and thus the desired non-trivial solution of the λ-condition (Equation 33) for M = 3 is obtained: Note that it is not guaranteed that the above solution always represents a non-trivial solution of the λ-condition (Equation 33). For example, the case of 1 = 0 (where 2 and 3 are not zero) implies the existence of a trivial solution (λ 1 , λ 2 , λ 3 ) = (1, 0, 0). In another example, the case of 1 = 2 = 0 implies that the denominator becomes zero; thus, the expression becomes indefinite.

PGF without winner-takes-all states
We are interested in those states in which any species does not take all molecules, because the actual simulations are performed by using the initial states excluding the winner-takes-all states. The PGF without the winner-takes-all states is represented by a linear summation of winner-takes-all states z N i and Equation (31) if the λ-condition has a non-trivial solution like Equation (39); where b 1 + · · · + b M+1 = 1. The take-all exclusion conditions are that the coefficients of z N i are zero; Therefore, we obtain the desired PGF without the winner-takesall states (PGFwoWTAS); which immediately implies the following stationary distribution in the case of catalytic reactions only; where N n 1 ,n 2 ,··· ,n M : = N! n 1 !n 2 !···n M ! are multinomial coefficients. Furthermore, we can calculate the marginal distribution of the i-th species from Equation (13); If we suppose the ergodicity n i = n i , the time-averaged concentrations can be derived from Equation (12); The above equations indicate that λ i means the concentration per total density ρ in the continuous limit N → ∞, that is, λ i should be the solution of the classical rate equation. We can also calculate the second-order moments x i x j and the variances (12); In the continuous limit N → ∞, the above Equations (45b) and (45c) become λ i λ j and 0, respectively, that is, the concentrations become mutually independent without fluctuating variables. We compare our formulas with simulation results that are obtained by applying the Gillespie algorithm (Gillespie, 1977) to the following three-species system ( Figure 1A) as an example: which implies the following by Equation (39): The marginal distributions are shown in Figure 2. Our formula Equation (44) is entirely in agreement with the simulation results. When the total number of molecules N is large, the marginal distributions can be approximated by normal distributions. Similarly, the formulas of the concentrations and the variances, Equations (45a) and (45c), are completely in agreement with the simulation results as shown in Figure 3. One can see that the variances monotonically decrease as N becomes larger.
Note that if the λ-condition (Equation 33) does not have a non-trivial solution like Equation (39), the expression for the PGF (Equation 42) cannot be applied. Such special cases are treated in the following section.

Ergodicity as a sufficient condition for applying our PGF
The PGFwoWTAS (Equation 42) would be applicable if the catalytic reaction network was "entirely ergodic, " which means the following in this paper (it is reminiscent of the ω-limit set):  Cross symbols (red, green, blue) represent the simulation results obtained numerically using the Gillespie algorithm (Gillespie, 1977), which is performed under the following condition; the total number of reactions: 10 8 , the number of reactions for transient exclusion: 10 7 , and the initial value (n 1 (0), n 2 (0), n 3 (0)) is randomly selected from W \ I such that the average per one-species is N/M. Empty symbols (circle, triangle, square) represent the theoretical expression Np c i * (Nx) [Equation (44)] for λ 1 = 2/11, λ 2 = 3/11, and λ 3 = 6/11.  Figure 1A. Empty symbols are numerically obtained using the Gillespie algorithm, where the time averages are performed over a long time series n i (t) such that the total number of reactions reaches 10 8 . The initial values are randomly selected from W \ I such that the average per one-species is N/M, and the number of reactions for transient exclusion is 10 7 . Lines in each figure represent the theoretical expressions, Equations (45a) and (45c), for λ 1 = 2/11, λ 2 = 3/11, and λ 3 = 6/11. One can see that the rank of concentrations is conserved but the rank of variances is exchanged between N = 2 and 4. where and ξ represents one of the trials of the stochastic process n(t) according to the catalytic reaction network. We use a specific three-species system to intuitively illustrate what Equation (48) means. In the case of the three-species system of Figure 1A, the ergodic condition (Equation 48) can be rewritten in terms of simplified conditions. As previously described, the state space of catalytic reaction networks consists of M N = (M+N−1)! (M−1)!N! points. Therefore, in the case of M = 3 for example, there are (N + 2)!/(2N!) points. The state space of the three-species system forms a regular triangle in three-dimensional Euclidean space (see Figure 4A). The possible motion from each state point is shown in Figure 4B, where FIGURE 4 | (A) State space of the three-species system of Figure 1A in the case of N = 10. Circles (red) represent the states that a single trajectory can visit, which means that the system is "entirely ergodic." (B) Motion of the state point when one reaction occurs in the three-species system of Figure 1A. There are 6 possible directions in which to move from one state point. Each direction is randomly selected in proportion to the specified probability. In the case of Figure 1A, the state point cannot move in the direction of R 321 since R 321 = 0. Note that the state point on the boundary (e.g., n 1 = 0) cannot move parallel to the boundary (in this case, the directions of R 213 and R 312 ).
there potentially exist six possible directions, but the network structure forbids the direction of R 321 in the case of the network of Figure 1A. Note that the state point on the boundary cannot move in a direction parallel to the boundary. In this case obviously, a trajectory starting from an initial point n 0 ∈ W \ I visits every point in W \ I as shown by the circles in Figure 4A (the set I is equivalent to the vertexes of a triangle). Therefore, as seen in Figure 4B, the necessary and sufficient condition to hold the ergodic condition (Equation 48) for M = 3 is: where the first condition represents that at least one direction for moving away from the boundary n k = 0 is allowed, and the second condition represents that at least one direction for approaching the boundary n k = 0 is allowed. Note that the condition (Equation 50) is no longer a sufficient condition for entire ergodicity in the case of four-species systems. In fact, in the following four-species system (Figure 5A), the condition (Equation 50) does not imply entire ergodicity; This four-species system obviously satisfies the condition (Equation 50), but most initial states [in particular, n k (0) > 0 (∀k)] eventually fall into any of 2-species winners-take-all states, i.e., {n 1 = n 2 = 0, n 3 + n 4 = N}.
In the case of M-species systems, a necessary and sufficient condition for entire ergodicity is difficult to derive, although we may consider a sufficient condition instead. We discuss the sufficient condition for entire ergodicity by introducing the collection of all l-species winners-take-all states (l = 1, 2, · · · , M): where We omit the subscripts M and N provided there is no concern for confusion. The above collections are not empty sets if l ≤ N. Note that I = I (1) and W = M l=1 I (l) . One can see that the following diagram holds if the catalytic network is fully connected [i.e., R ijk > 0 for all i, j, k (i = j = k = i)]; where each arrow represents the direction in which a state point can move by one reaction. Each I (l) (i 1 , · · · , i l ) is a (l − 1)-dimensional object in the whole (M − 1)-dimensional state space W. For example, the state space of four-species systems forms a regular tetrahedron W M = 4,N (see Figure 5B); I (4) (1, 2, 3, 4) is the interior of the regular tetrahedron; I (3) (1, 2, 3), I (3) (1, 2, 4), I (3) (1, 3, 4), and I (3) (2, 3, 4) are regular triangles that form the boundaries of I (4) ; and I (2) (1, 2), I (2) (1, 3), I (2) (1, 4), I (2) (2, 3), I (2) (2, 4), and I (2) (3, 4) are line segments that form the boundaries of I (3) . Note that I (2) (i 1 , i 2 ) → I (2) (j 1 , j 2 ) for (i 1 , i 2 ) = (j 1 , j 2 ) is possible but I (2) (i 1 , i 2 ) → I (2) (i 1 , i 2 ) is always impossible as seen from Figure 5C. By considering the diagram (Equation 54) and keeping in mind Figures 5B,C, we can expect the following fact: if entire ergodicity holds in the restricted system on I (l) for l ≥ 3, then entire ergodicity holds in the whole state space W. In other words, in systems consisting of M-species, the sufficient condition under which the ergodic condition (Equation 48) holds can be geometrically described as follows: on each I (l) (i 1 , · · · , i l ) for l ≥ 3, at least one direction for moving away from each boundary I (l−1) (j 1 , · · · , j l−1 ) should be allowed, and simultaneously at least one direction for approaching each boundary I (l−1) (j 1 , · · · , j l−1 ) should also be allowed, where {j 1 , · · · , j l−1 } ⊂ {i 1 , · · · , i l }. . Red and blue arrows represent upward and downward arrows, respectively. There are 12 possible directions in which to move from one state point. Each direction is randomly selected in proportion to the specified probability. Note that the state point on each I (2) (i 1 , i 2 ) cannot move parallel to itself (e.g., on I (2) (1, 2), the directions of R 231 , R 241 as well as R 132 , R 142 are not allowed).

The Case of Catalytic-Noncatalytic Mixed Reactions
We could not obtain the steady-state solution of the general GFE (Equation 11) in the case of catalytic-noncatalytic mixed reactions (ε > 0), but we expect our PGF (Equation 42) to be a good approximation for mixed-reaction systems if ε is sufficiently small (ε ≪ min{R ijk > 0}). More specifically, we expect the PGF (Equation 42) to be robust against non-catalytic reactions if the catalytic reaction system constituting the mixed reaction system has an ergodic component spread across the entire state space (entire ergodicity). Otherwise, if the catalytic reaction system has several ergodic components in the entire state space, the non-catalytic reactions may imply that a certain highly stable ergodic component attracts all possible trajectories, which obviously means that our PGF is not applicable to such a mixed reaction system. Although we were unable to prove this mathematically, we used numerical simulations to determine whether our PGF is a good approximation in the case of "entirely ergodic." Figure 6 represents the non-catalytic reaction rate constant ε dependencies of time-averaged concentrations x i for the threespecies system of Figures 1A,B. (As shown later, the system shown in Figure 1B is not entirely ergodic, which is the reason why the systems shown in Figures 1A,B are compared here.) It can be seen that the differences of x i in Figure 6A between ε = 0 and ε = 0.1 are smaller than those in Figure 6B. We consider the robustness in the case of Figure 6A to originate from the entire ergodicity of Figure 1A, which means that the ergodic component spreads on the entire state space except for winner-takes-all states (see Figure 4A).

Rank Conservation Law for Concentrations
We show that the rank of concentrations is conserved even if the total number of molecules changes in catalytic reaction networks (excluding non-catalytic and auto-catalytic reactions).
Suppose the concentration of the ith-species is expressed by Equation (45a) in the state corresponding to the PGFwoWTAS. When determining the rank conservation, it suffices to confirm that the relation between the amount of two arbitrary species is unchanged if the total number of molecules is changed. Let λ 1 , λ 2 be the concentrations per total density in the continuous limit such that λ 1 < λ 2 . Because M i=1 λ i = 1, the inequality λ 1 +λ 2 < FIGURE 6 | ε-dependence of the time-averaged concentrations x i (N = 10, ρ = 1) in the case of (A) the three-species system of Figure 1A with non-catalytic reactions and (B) that of Figure 1B with non-catalytic reactions. Points in the figures are numerically obtained using the Gillespie algorithm in the same way as Figure 3. The catalytic reaction network of Figure 1A is entirely ergodic; thus, the non-catalytic reactions may be treated as perturbation. On the other hand, that of Figure 1B is non-ergodic; thus, the non-catalytic reactions introduce relatively large differences compared with the case of ε = 0.  Figure 1C. Empty symbols are numerically obtained using the Gillespie algorithm in the same way as Figure 3. Lines in each figure represent the theoretical expressions, Equations (45a) and (45c), for λ 1 = 1/1000, λ 2 = 1/3, and λ 3 = 1997/3000. It is clear that the rank of time-averaged concentrations is conserved.
1 must be satisfied. Therefore, the following evaluation holds: This is the rank conservation law of concentrations. Note that the rank of the variances of concentrations is generally not conserved when the total number of molecules changes. For example, let us consider the following three-species system ( Figure 1C); R 123 = 1997/3, R 132 = 1000/3, R 231 = R 321 = 1, As shown in Figure 7, this system depends on the total number of molecules for the time-averaged concentrations x i and the variances of concentrations Var[x i ], represented by Equations (45a) and (45c) with λ 1 = 1/1000, λ 2 = 1/3, and λ 3 = 1997/3000. The rank of time-averaged concentrations is always conserved, but the rank of variances is exchanged at certain N (see Figure 8).

A Special Case: The Connecting State to the Winner-Takes-All State
There exists an N-molecular state, which connects to the winnertakes-all state in the continuous limit N → ∞. As is later shown, such a special case is not the "weakly reversible" case (Anderson et al., 2010). We show this by taking the following limit in the PGFwoWTAS (Equation 42):  Figure 1C.
The ranks (lines) in the figures are depicted using the formulas (Equations 45a and 45c). Clearly, the rank of concentrations is conserved, although the rank of variances is exchanged between N = 2 and 3, and also N = 6 and 7.
where we suppose the following constraints are always satisfied: where {κ i } are positive constants, which are determined from the network structure {R ijk } (it is explained later). Evidently, the following holds; The PGFwoWTAS (Equation 42) has the following limiting expression: The state corresponding to the PGF (Equation 60) is the connecting state to the winner-takes-all state (CStoWTAS). The stationary distribution corresponding to the CStoWTAS is immediately obtained: where J M,N (abbr. J) is defined as Furthermore, the marginal distributions of the i-th species can be derived by Equation (13): p cs i * (n) = δ n,N−1 for i = 1, (1 − κ i )δ n,0 + κ i δ n,1 for i = 2, · · · , M, The time-averaged concentrations calculated by Equation (12) are and the variances of the time-averaged concentrations calculated by Equation (12) are Var Next, we derive the relation between the positive constants {κ i } and the network structure {R ijk }. The λ-condition (Equation 33) can be converged to conditions for {κ i } (the κ-condition) by taking the limit λ 1 → 1 as follows. Divide the λcondition (Equation 33) into two groups, of which one group is the case of 2 ≤ i < j ≤ M, and the other group is the case of where we substituted λ i = κ i (1−λ 1 ), i ≥ 2 and divided by 1−λ 1 . Then, taking the limit λ 1 → 1 in Equation (66a), and taking that in Equation (66b), The condition (Equation 67a) expresses that the 1st-species cannot be a substrate, and the condition (Equation 67b) represents the desired κ-condition. Note that the CStoWTAS must be a limiting state corresponding to the PGFwoWTAS [see Equation (60)]. Therefore, the network structure {R ijk } must have definite {λ i }.
For example, let us consider the following three-species system ( Figure 1B): which implies κ 2 = 2/3 and κ 3 = 1/3. Obviously, this system is not weakly reversible, which means that the theorems (Theorem 4.1 and 4.2) in the previous study (Anderson et al., 2010) cannot be applied. This system should have a CStoWTAS because λ 1 , λ 2 , and λ 3 are definite from Equation (39) with 1 = 0, 2 = 2, and 3 = 4. As shown in Figure 9, this system certainly has the time-averaged concentrations x i represented by Equation (64) and the variances Var[x i ] represented by Equation (65). In this case also, the rank conservation law of concentrations holds.

M-Species 2-Molecules Systems with Non-catalytic Reactions
The 2mTESM (Equation 20) becomes the closed equation of the second-order moments x l x m if the first-order moments x i are substituted by the second-order moments according to the PRE (Equation 15). In this subsection, we consider catalytic-noncatalytic mixed reaction systems of M species, consisting of only 2 molecules in total.
We first focus on the second formula in the 2mTESMs (Equation 20). It can be seen that each variance of time-averaged concentration Var[x i ] = x 2 i − x i 2 in the steady state depends solely on its concentration: where we supposed the ergodicity x i = x i . The above equation expresses that the fluctuation of the concentration x i becomes larger as its time-averaged concentration approaches at which the fluctuation takes the largest value M+1 4M ρ 2 .
Furthermore, the time-averaged concentrations potentially have the maximum value because the variance must be positive; Next, we focus on the first formula in the 2mTESMs (Equation  20). Let us consider the steady state of Equation (20a) and suppose the ergodicity x i = x i . Eliminating the first-order moments in Equation (20a) by using the SME (Equation 18) gives the determination equation of second-order moments in the 2 molecules system (2mDESM); The above 2mDESM can be rewritten in the M(M−1)/2×M(M− 1)/2 matrix form. We demonstrate the procedure for processing  Figure 1B. Empty symbols were numerically obtained using the Gillespie algorithm in the same way as Figure 3. Lines in each figure represent the theoretical expressions, Equations (64) and (65), for κ 2 = 2/3 and κ 3 = 1/3. FIGURE 10 | ε-dependence of (A) the second-order moments x i x j = n i n j /N 2 and (B) the time-averaged concentrations x i = n i /N in the three-species system of Figure 1A with non-catalytic reactions (N = 2, ρ = 1). Empty symbols are numerically obtained using the Gillespie algorithm in the same way as Figure 3. Lines in each figure represent the theoretical expressions, Equations (74) and (75). Note that there exist differences between the limit ε → 0 and the case ε = 0 for the second-order moments (solid black points) because the SME (Equation 18) breaks down at ε = 0. The solid black points are calculated from Equation (45b) as λ 1 = 2/11, λ 2 = 3/11, λ 3 = 6/11, N = 2, and ρ = 1. the 2mDESM by considering the case of M = 3. In this case, the 2mDESM (Equation 72) has a matrix of the following form: Moreover, we consider the specific three-species system of Figure 1A (ρ = 1) including non-catalytic reactions. Solving Equation (73) as ε > 0 and ρ = 1 in the case of Figure 1A, then which is consistent with the result obtained when solving the CME (Equation 7) directly [recall x 1 x 2 = 1 4 P(n 1 = 1, n 2 = 1, n 3 = 0), x 1 x 3 = · · · and so on]. By the SME (Equation 18), the concentrations become .
(75) Figure 10 shows the second-order moments x i x j and the timeaveraged concentrations x i as functions of ε. Note that in the case of ε = 0 (catalytic reactions only), we need to derive the secondorder moments x i x j from Equation (45b) corresponding to λ 1 = 2/11, λ 2 = 3/11, and λ 3 = 6/11. The non-catalytic reaction rate constant ε seems to be a singular perturbation against the second-order moments (not the concentrations). . If one regards the species 1, 3, and half of 5 as the species A (similarly, 2, 4, and half 5 as B), the behavior of n A = n 1 + n 3 + n 5 /2 and n B = n 2 + n 4 + n 5 /2 is similar to that of the 2TK model.

Non-autocatalyzation of Autocatalytic Reaction Networks
The framework we developed in this paper applies to nonautocatalytic reaction networks. However, our framework may be applicable to autocatalytic reaction networks if it were possible to convert autocatalytic to non-autocatalytic networks. Here, we show several examples of such conversions using a minimal autocatalytic reaction network "2TK model" (Ohkubo et al., 2008;Saito and Kaneko, 2015), which is a well-studied model in the context of discreteness-induced phenomena. The 2TK model consists of only two species, and includes both autocatalytic reactions (rate const. r) and non-catalytic reactions (rate const. ε ≪ r) (see Figure 11A); A ⇄ B (non-cat. react.).
If we suppose that each of the species A and B consists of two further species, we can convert autocatalytic reactions Equation (76) of the 2TK model to non-autocatalytic reactions as shown in Figure 11B. The catalysis of the species 5 between the species 1 and 3 (or, 2 and 4) is to establish an equilibrium between the concentrations of the species 1 and 3 (or, 2 and 4), that would make it possible to regard these as one species. The nonautocatalytic reactions of Figure 11B plus non-catalytic reactions i Prob. 1/5 − −−−− → j is the desired five-component model, of which the CME is simply described by Equation (7) with As shown in Figure 12, the behavior of the variables n A = n 1 + n 3 + n 5 /2 and n B = n 2 + n 4 + n 5 /2 in the five-component model is similar to that of the 2TK model (compare with Figures 1A,B in Saito and Kaneko, 2015). The λ-condition (Equation 33) of catalytic reactions in the five-component model can be written in matrix form, in which there exist six non-trivial eigenvectors corresponding to the 0-eigenvalue; hence, the general eigenvector corresponding to the 0-eigenvalue can be expressed as a linear summation of those six non-trivial eigenvectors: The above equation can be used to derive six types of non-trivial solutions (λ 1 , λ 2 , · · · , λ 5 ), which immediately correspond to the stationary states of catalytic reactions in the five-component model by using Equation (43): where each case (i-vi) corresponds to (i) c 3 = c 2 , c 6 = c(1 − 2c) (others 0), (ii) c 2 = c 2 , c 5 = c(1 − 2c) (others 0), (iii) c 1 = 2c(1 − 2c) (others 0), (iv) c 2 = 2c(1 − 2c) (others 0), (v) c 3 = 2c(1 − 2c) (others 0), and (vi) c 4 = 2c(1 − 2c) (others 0). The switching behavior of the five-component model may be explained in terms of transition processes between the above six types of steady states and the five trivial steady states λ i = 1 (others 0) of catalytic reactions in the five-component model, which are sometimes caused by non-catalytic reactions.  Figures 1A,B in Saito and Kaneko, 2015). (A) Time series of the total concentration of the species 1, 3, and half of 5 for N = 20 (red line) and N = 2000 (green line). (B) Stationary distributions of (n 1 + n 3 + n 5 /2)/N, obtained numerically from a long time series n i (t) using the Gillespie algorithm whose simulation conditions are the total number of reactions: 10 8 , the number of reactions for transient exclusion: 10 7 , and the initial value is randomly selected from W \ I such that the average per one-species is N/M. The unimodal distribution (blue), the flat distribution (green), and the bimodal distribution (red) indicate the stationary distribution for N = 500, 40, and 20, respectively. Furthermore, the marginal distribution of the species A, p A (n), can be derived by calculating the convolution of the marginal distributions p 1 (n), p 3 (n), and p 5 (n), if p 1 (n), p 3 (n), and p 5 (n) were obtained for the case including non-catalytic reactions.
Other non-autocatalytic reaction networks duplicating the 2TK model are shown in Figure 13A, and consist of 4 components only. Although readers would think that the four-component models do not function properly because the transitions between the species group (1, 3) [or (2, 4)] depend on the other species group (2, 4) [or (1, 3)], the four-component models actually generate much the same behavior with the 2TK model (see Figure 14). The switching behavior of the fourcomponent model may also be explained as a transition process between the following three types of steady states and the four trivial steady states λ i = 1 (others 0) of catalytic reactions in the four-component model, which are sometimes caused by non-catalytic reactions. We also confirmed . If one regards the species 1 and 3 as the species A (similarly, 2 and 4 as B), the behavior of n A = n 1 + n 3 and n B = n 2 + n 4 is almost equivalent to that of the 2TK model.
that the results are not changed even in the other fourcomponent model of Figure 13B, of which the catalysis role of each the species 1 and 2 is exchanged by that of 3 and 4, respectively.

SUMMARY AND DISCUSSIONS
The framework we presented in this paper facilitates the prediction of the effect of the small-number issue on the concentration of each species in catalytic reaction networks. This can be described in an extreme manner by comparing the concentrations between the continuous limit (N → ∞) and the case of 2 molecules (N = 2). If the reaction network does not include non-catalytic reactions (or, includes negligible non-catalytic reactions), we can use the formula (Equation 45a) to compare them. On the other hand, if the reaction network includes (non-negligible) non-catalytic reactions, we need to apply the formula, 2mDESM (Equation 72), with the SME (Equation 18) to obtain the concentrations in the case of 2 molecules. Although our theory has a presupposition referred to as entire ergodicity, the presupposition is intuitively verifiable if the system is specified, as in Section 2.4.2.4. We also demonstrated three examples for non-autocatalyzation conversions of autocatalytic reaction networks in Section 3.4. We consider this type of conversion to be generalized relatively easily such that our analytical framework can be applied to more general catalytic reaction networks including autocatalytic reactions. One might think that the analysis presented in this paper can be straightforwardly extended to the case including autocatalytic reactions [in fact, the CME (Equation 7) and GFE (Equation 11) themselves hold even in the case including autocatalytic reactions]. However, if autocatalytic reactions are included (i.e., the case R ikk > 0 is allowed), we cannot consider catalytic reactions and non-catalytic reactions to be completely separate. The reason is that, in the case including autocatalytic reactions,  Figures 1A,B in Saito and Kaneko, 2015). (A) Time series of the total concentration of the species 1 and 3 for N = 20 (red line) and N = 2000 (green line). (B) Stationary distributions of (n 1 + n 3 )/N, obtained numerically from a long time series n i (t) using the Gillespie algorithm in the same way as Figure 12. The unimodal distribution (blue), the uniform distribution (green), and the bimodal distribution (red) indicate the stationary distribution for N = 500, 40, and 20, respectively. the absence of non-catalytic reactions generally implies winnertakes-all steady states. Generally, solving the CME (or GFE) of catalytic-noncatalytic mixed reactions systems is more advanced and a more difficult task than that of catalytic reactions only. The proposed strategy, i.e., non-autocatalyzation conversions, is one of our ideas to address the problem.
The formulas obtained in the present work are specific and satisfactorily simple. Therefore, our theory has the capabilities to be developed into a general theory for catalytic reaction networks. On the other hand, there exists a mathematical theory for a certain class of catalytic reaction networks that are "weakly reversible" and "deficiency zero" (Anderson et al., 2010). Our formulas (Equations 27 and 43) are consistent with the main theorems (theorem 4.1 and 4.2 in Anderson et al., 2010) in the above-mentioned mathematical theory. One of the advantages of our theory compared to the above mathematical theory is understandably the probability generating function (PGF) approach, because the PGF is a major analytical tool for physicists, chemists, and mathematical biologists. Therefore, our theory is easily verifiable, and one can design a computer algorithm to calculate our analytical formulas. We also showed the extensibility of our theory by using applications (in Section 3), especially because of the CStoWTAS (Equation 61), which was not suggested by the above mathematical theory.
Actual biochemical pathways in the cell involve thousands of chemical species, and their chemical properties vary. Our theoretical framework is general and extensible to such complex reaction networks, if they can be represented by CMEs such as Equation (7). As our current model consists of simple twobody catalytic reactions, it is difficult to point out examples in actual biological systems that correspond exactly to our model. Biochemical reactions in reality may involve a number of intermediates. There are also autocatalytic processes such as autophosphorylation, and replication of templates such as DNA, in which the catalyst or template species is also a substrate or a product. Our framework is applicable to many such cases involving network conversion, as shown for simple autocatalytic cases.
Nevertheless, the reaction kinetics of each enzyme is not always simple. Enzymes are complex macromolecules and their reaction cycles may depend on their conformational states. Therefore, the prediction of biological phenomena caused by small-number effects in real biochemical reactions, would entail further analytical challenges for catalytic reaction networks including arbitrary higher-order mixed reactions (rather than first-and second-order reactions only) or internal dynamics of the enzymes (as modeled and analyzed in Togashi and Casagrande, 2015) as important issues.
Throughout this work, our primary intention is to approach small-number issues in biological systems. One might wonder how general these small-number issues appear, and how important they are, in living cells. Recently, absolute quantification of various proteins and mRNAs in the cell has become possible, and the integration of experimental results (e.g., the construction of a database Milo et al., 2010) is also underway. Taniguchi et al. investigated the copy number distribution for more than a thousand protein species in bacteria (Taniguchi et al., 2010). Li et al. further discussed the relationship between the copy number and synthesis rate, and also the role, of proteins (Li et al., 2014). According to the result, some transcription factors, particularly activators, are rare, of the order of 0.1 to 10 molecules per genome equivalent. Although stochastic gene expression has been intensively discussed for years (McAdams and Arkin, 1997;Thattai and van Oudenaarden, 2001;Elowitz et al., 2002;Raj and van Oudenaarden, 2008;Shahrezaei and Swain, 2008), the discrete small-number nature of transcription factors has often been ignored; hence, the finding may urge us to reconsider the issue. Synthetic approaches are also becoming popular to confirm small-number effects. Ma et al. reported that an additional stable state in a genetic bistable toggle switch attributable to the small-number effect, which was predicted by stochastic simulations, was indeed observed in bacteria containing the genetically engineered switch (Ma et al., 2012). These results suggest that such rare proteins, of the order of one molecule per cell, are common and affect regulatory function in bacteria.
Although eukaryotic cells are much larger than bacteria, they have complex membrane structures and cytoskeletons inside, and the small-number issues can be particularly significant in compartments or bottlenecks (e.g., if we consider the volume of a synaptic vesicle represented by a sphere 40 nm in diameter, then 1 molecule corresponds to ca. 50 µmol/L). Rare proteins are also involved in physiologically important signaling pathways in eukaryotes. In the Wnt signaling pathway, for example, the concentration of axin is reported to be 20 pmol/L in Xenopus eggs (Lee et al., 2003) (though suggested to be higher in mammalian cells Tan et al., 2012). Another example is the MAP-kinase cascade, where proteins of the order of merely 10 2 molecules (e.g., Ste5) exist in a yeast cell (Thomson et al., 2011). For scaffold proteins such as axin and Ste5, specifically, localization (locally high concentrations) of other chemical species around the scaffold may drastically change the reaction behavior, as spatial discreteness of the scaffolds becomes significant (Shnerb et al., 2000;Togashi and Kaneko, 2004). Further studies in which spatial structures (cf. reaction-diffusion equations) are considered are also expected to be important.
In the presented framework, we mainly focused on the steady-state solutions of GFE. Of course, temporal courses are biologically crucial in some cases. A well-studied example is oscillatory behavior in circadian clocks (Bell-Pedersen et al., 2005). In such oscillations, if a chemical factor is depleted down to a small number in a certain phase, then, the period can be susceptible to stochastic reactions involving the factor in that phase; on the other hand, sequestration of a factor may contribute to regular oscillations (Jolley et al., 2012). Again, the internal dynamics of enzymes can also be relevant in some systems. Although stochastic simulations are indeed powerful and many attempts are currently underway, further theoretical understanding as well as the experimental quantitative observation of rare factors would be required.
Note that a chemical "species" here can also be interpreted as a specific state of a molecule; e.g., we can consider proteins or genes, with and without modification, as separate species.
Moreover, a similar interpretation is also applicable to ecology and ethology (Biancalani et al., 2014), if the laws governing the system are analogous to reactions. We remain hopeful that theoretical frameworks including ours will facilitate the exploration of small-number issues at equally higher levels of biological systems in future.