Chemical Kinetics Roots and Methods to Obtain the Probability Distribution Function Evolution of Reactants and Products in Chemical Networks Governed by a Master Equation

In this paper first, we review the physical root bases of chemical reaction networks as a Markov process in multidimensional vector space. Then we study the chemical reactions from a microscopic point of view, to obtain the expression for the propensities for the different reactions that can happen in the network. These chemical propensities, at a given time, depend on the system state at that time, and do not depend on the state at an earlier time indicating that we are dealing with Markov processes. Then the Chemical Master Equation (CME) is deduced for an arbitrary chemical network from a probability balance and it is expressed in terms of the reaction propensities. This CME governs the dynamics of the chemical system. Due to the difficulty to solve this equation two methods are studied, the first one is the probability generating function method or z-transform, which permits to obtain the evolution of the factorial moment of the system with time in an easiest way or after some manipulation the evolution of the polynomial moments. The second method studied is the expansion of the CME in terms of an order parameter (system volume). In this case we study first the expansion of the CME using the propensities obtained previously and splitting the molecular concentration into a deterministic part and a random part. An expression in terms of multinomial coefficients is obtained for the evolution of the probability of the random part. Then we study how to reconstruct the probability distribution from the moments using the maximum entropy principle. Finally, the previous methods are applied to simple chemical networks and the consistency of these methods is studied.


Introduction
The early studies on stochastic kinetics and the chemical master equation (CME) started in the early sixties with the works of McQuarrie and Montroll [1,2], followed two year later by the paper of Nicolis and Babloyantz [3]. Then, Nicolis and Prigogine [4], Nicolis and Malek-Mansour [5], Gillespie [6,7] and Van Kampen [8] obtained the complete form of the CME and, at the same time, these researchers developed tools for solving this equation at different degrees of approximation, for several simple chemical network cases.
The chemical reaction networks can be considered as Birth-and-Death processes without memory, which are Markov processes in a multidimensional space, this point of view was studied by in detail in Section 3.1, is the method of the probability generating function or z-transform. It is obtained in this section the equation obeyed by the probability generating function for the reactant and product population of the network, while Section 3.2 is devoted to show the method of the expansion of the master equation in one order parameter, general equations are obtained at different orders of approximation in terms of multinomial coefficients.
In Section 4, we study the methods to reconstruct the probability distribution function (PDF) from the moments using the maximum entropy principle. Finally, Section 5 is devoted to the application of the previous methods to study several simple reaction networks with non-linear effects in open and closed systems. Also, in this section we check the consistency of some of the approximations performed, finishing this section with the closure problem and discussing how to perform the closure in an optimal way.

Introductory Concepts on Chemical Reactions Networks and the Chemical Master Equation
In this section, we give a brief introduction about the main concepts on chemical reaction kinetics, chemical reaction networks, chemical reaction propensities, to finish with the chemical master equation and its limitations.

Some Introductory Concepts about Chemical Reactions Networks and Entropy Production
Let us consider a stirred tank of given volume Ω containing n S chemical compounds or species X j . Let us denote by N j the number of molecules of the j-th chemical compound. The number of molecules of all the chemical species will be represented as a point N belonging to a N n s multidimensional space, being N the set of the natural numbers, so that is possible to write: Let us consider now a set of n r chemical reactions and their inverses: 2 X 2 + . . . r = 1, 2, . . . , n r (2) l (r) 2 X 2 + . . . r = n r + 1, . . . 2n r That specifies, in general, a chemical reaction network. The coefficients of Equations (2) and (3) i.e., s (r) j , l (r) j for a given reaction r are known as the stoichiometric coefficients of the reaction. In general, each reaction of an arbitrary type r that takes place inside the tank, change the state of the system, which moves from the point N ∈ N n s to the point N ∈ N n s given by: where in Equation (4), the vector ν (r) = l which take place when one chemical reaction of the r-th type occurs inside the system volume Ω. If the inverse reaction of r happens inside the system then the final state is N = N − ν (r) . Now, if we denote by ξ r the number of reactions of the r-th type that take place inside the system during a given time interval [t 0 , t = t 0 + ∆t], and by . ξ r the reaction rate. Then, it follows that if at time t 0 , the system state was N 0 , then the state of the system at time t = t 0 + ∆t will be given by the following expression: The total number of reactions that can take place inside a closed tank is limited by the available population of each species of atoms, denoted by the superscript α. Then if the symbol n for α = 1, 2, . . . ., n ατ and r = 1, 2, . . . , 2n r , (6) being n ατ the total number of different elements contained in the tank. Next, Equation (6) can be rewritten as follows: (r) ·n (α) = 0 for α = 1, 2, . . . , n ατ and r = 1, 2, . . . .2n r , where ν (r) ·n (α) denotes the scalar product in N n s . Equation (7) is true for all the reactions types and all the kinds of atoms. Because the number of atoms of each element inside the tank is conserved, then in a closed system we can express this fact as follows: for α = 1, 2, . . . , n ατ , (8) where N(t) ≡ (N 1 (t), N 2 (t), . . . , N n s (t)), is the number of molecules of each chemical compound inside the tank at a given time t, and N (α) 0 is the total number of atoms of the element α, which are available initially at time t = 0. If the system is open this number can change with time and the conservation law (8) transform in: Equation (9) expresses the fact that the number of atoms of each element can change with time in an open chemical system in which specific types of molecules can be injected or extracted along the time.
Equations (8) and (9) define a set of n ατ hyperplanes in the space N n s . The part of the space N n s , which is accessible starting from the initial point N 0 for a closed system is in the intersection of the hyperplanes defined by Equation (8). Notice that not all points are accessible due to the restrictions imposed by the finite number of atoms of each element.
The degree of advance dξ r of a chemical reaction r during an infinitesimal time interval dt can be expressed in terms of the change in the number of molecules of the chemical compounds as follows [35]: Therefore, the change dN (r) j in the number of molecules of the chemical compound X j , produced by a reaction r with degree of advance dξ r is obviously given, because of Equation (10), by: From the Gibbs relation, the total change in the system entropy produced by the reactions taking place in the tank with volume Ω, when not external exchange of the chemical compounds (closed system) occurs, is given by the expression: where µ j is the chemical potential for the compound X j , and dS = d i S + d e S, is the sum of the internal and external entropy changes. Notice that the external change can be produced by external exchanges of chemical compounds and/or energy. Using Equation (12) in Equation (11) and defining the chemical affinity for the r-th reaction in the form: It is obtained the rate of change of the entropy with time: For non-equilibrium system thermodynamics processes, it is important to know their thermodynamic fluxes and their generalized thermodynamic forces. To do this, we consider the rate of change of the entropy per unit volume in an isolated system with constant volume and constant internal energy. This production rate of entropy denoted by σ, is given because of Equation (14) by: From Equation (15), it is obtained that thermodynamic forces are the affinities and the thermodynamic fluxes are the reaction rates per unit volume, which are given because of Equation (15) by: If the tank is well stirred, we can consider the system as homogeneous because matter flux of chemical compounds produced by the system in-homogeneities does not exist. Because of the thermodynamic fluxes depend, in first order approximation, linearly on the affinities i.e., J (r) = ∑ 2n r r=1 L r,s A (s) , then it follows that at chemical equilibrium (σ = 0), i.e., for not evolution of all the chemical reactions inside the system we must have: From conditions (17) and the expression for the chemical potentials given by [35][36][37]: where, a j is the chemical activity of the j-th component and µ is the chemical potential of the compound in a pure substance. Then, it is easily deduced because of Equation (13) the mass action law for the direct chemical reaction: where K (r) (p, T) is the equilibrium constant for the reaction r, at given pressure p and temperature T, being a 1, a 2 , . . . , a n s the activities of the chemical compounds at equilibrium conditions denoted by the overbar. From Equations (13), (18) and (19), it is deduced that the equilibrium constant is given by the expression: It is important to realize that a j = γ j x j are the chemical activities for non-ideal systems, being γ j the activity coefficients, and x j the molar fractions of the chemical compounds.

Probability Distribution Function (PDF) at Equilibrium
In this section we study the PDF of the number of molecules of each chemical compound at equilibrium conditions. Einstein was the first to invert the classical Boltzmann relation between entropy and probability, obtaining the probability in terms of the system entropy [4,38]. Let us consider a system that depends on a set of thermodynamic magnitudes as the internal energy U, the volume Ω and the number of molecules N j of each chemical compound X j , in this case we can write S = S(U, Ω, N 1 , N 2 , . . . , N n s ). In addition, according to the second principle at equilibrium conditions, the variables take the values that maximize the system entropy, and using statistical mechanics it is possible to obtain the probability of fluctuations of the variables around the maximum entropy state, described by S = S U, Ω, N 1 , N 2 , . . . N n s . If we assume that the tank volume remains constant, then from Einstein relation, we can write: being k the Boltzmann constant, and P the probability at equilibrium conditions. Next, we use the Gibbs relation (12) at constant U and Ω, that yields, if the system is closed: where j is the change in the number of molecules belonging to j-th species produced by all the reactions occurring inside the stirred tank. Expanding the entropy in Taylor series up to second order around the maximum entropy state, it is found that: Next from the Gibbs relation it is deduced that the chemical potential is given by: Finally, from Equations (21), (23) and (24), it is obtained that the probability of fluctuations in the number of chemical compounds around the equilibrium state is given by:

Roots of Chemical Reaction Kinetics
We start studying the case of the reaction rate of two different chemical compounds, X i and X j , which can react producing a chemical compound, X k . If we have at time t, N i and N j molecules of these chemical compounds inside the tank with volume Ω and temperature T, and we denote by σ ij the collision cross section of both molecules, and by Then, if the velocity of the individual particles follows a Maxwell-Boltzmann distribution, it is easy to prove that the relative velocity of the particles follows also a Maxwell-Boltzmann distribution, which has the form: where ε CM = 1 2 mv 2 r is the center of mass (CM) energy of both particles, being m = m i m j m i +m j the reduced mass. To produce the reaction, we need at least a minimum threshold energy in the CM system denoted by ε CM,0 = 1 2 mv 2 r,0 , so that, if the CM energy of the colliding compounds is smaller than ε CM,0 , then the reaction does not take place. So, if the tank is well stirred and at temperature T then we can compute the collision R c,ij and reaction rates R r,ij per unit volume by means of the following expressions: where we have considered that d 3 v r = v 2 r dv r dΩ solid = v 2 r dv r 2π sin θdθ. Substituting Equation (26) into Equation (27) and performing the integrals in Equation (27), (see the Handbook of Mathematical Functions section 7 on Error functions and Fresnel integrals [39]), yields after some calculus: where, v r,ij = 8πkT m 1/2 is the average relative velocity between the molecules at temperature T. But not all the collisions between the molecules, X i and X j produce reactions X i + X j → X k . Because only those collisions with energy in the center of mass system bigger than the critical energy will yield reactions. Therefore, to obtain the reaction rate per unit volume one must consider the lower limit of the velocity in the integral of Equation (27) to be equal to the velocity corresponding to the critical energy. In addition, according to Gillespie [6], the reaction will take place only if the collisional point of contact takes place inside the sensitive area of each colliding molecule, this fact introduces an additional factor, that will be denoted by p (s) i,j . This factor will depend on the subtended solid angle of the sensitive area for each molecule. So, the reaction rate per unit volume for bimolecular reactions will be given by: where p (s) i,j is given by the expression: being ω (s) i and ω (s) j the solid angles which are subtended by the reaction sensitive areas of molecules i and j respectively. So, the reaction rate for bimolecular collisions can be expressed because of Equations (28) and (29) as follows: where η ij is the probability that one collision between two chemical compounds X i and X j yields one reaction event, this probability is given by the expression: The propensity function γ r (N) for this type of bimolecular reaction r is defined in such a way that γ r (N)dt gives the probability to have one reaction event of type r in the time interval between (t, t + dt) inside the stirred tank of volume Ω assuming that the system was at time t in the state characterized by the state vector N. On account of this definition and expression (31), we can compute the propensity for a bimolecular reaction of different species as: Therefore K r dt gives the probability to have one reaction of type r between t and t + dt per i, j pair and per unit volume. Notice that the total number of i, j couples per unit volume is If the chemical compounds that collide are identical then, m i = m j , and the number of different couples that can collide inside Ω is N i (N i − 1)/2, then proceeding as previously it is obtained that the collision rate per unit volume is: where σ col,ii denotes the collision cross-section between two identical molecules, m = m i 2 is the reduced mass and v r,ii denotes the average relative velocity between two identical particles. The rate of collisions per unit volume with energy in the CM system bigger than the critical energy and that produce reactions X i + X i → X l is: where the expression for η ii is given by Equation (32). The probability to have one reaction X i + X i → X l in the differential interval (t, t + dt), if the system was at the state N at time t, is given on account of Equations (34) and (35) by: where γ r (N) denotes the propensity for this chemical reaction, and K r dt gives the probability to have one reaction in the tank per couple of identical chemical compounds per unit volume.
Generalizing these previous expressions to a generic reaction r: ∑ n s j=1 s the propensity can be expressed by the following equation: For a bimolecular reaction of identical particles, we have that s (r) j = 2, and Equation (37) yields Equation (36) for this case. For a bimolecular reaction of different chemical compounds X i , and X j s (r) i = 1 and s (r) j = 1, and Equation (37) gives for this case Equation (33).

The Chemical Master Equation
Let us consider the set of reactions ∑ n s j s j , with n r = 1, 2, . . . , n r , with kinetic constants denoted by K r and their inverses ranging from r = n r+1 + 1, . . . , 2n r with kinetic constants denoted by K n r +1 = K 1 , K n r +2 = K 2 , . . . K 2n r = K n r . We notice that the change vectors ν (r) of the population of the different compounds produced by the inverse reactions are related to their corresponding change vectors of the direct ones by the set of equations: Next one defines a continuous time Markov process {N(t), t ≥ 0} with discrete state space N(t) ∈ N n s . The main characteristic of these processes is that the system state depends only on the last previously known state and not on the previous ones at earlier times. Therefore, the n-times conditional probability for a Markov process reduces to [40]: P n (N n , t n |N n−1 , t n−1 ; . . . ; N 1 , t 1 ) = P 2 (N n , t n |N n−1 , t n−1 ) ∀t n ≥ t n−1 ≥ . .
The following equation, known as the Chapmann-Kolmogorov equation is true for all the Markov processes [40], and it is easily deduced because of the definition of conditional probability: where the summation extend over all possible states N at an intermediate time t . To deduce the Master equation one consider a time instant t = t − ∆t, arbitrarily close to t, and denoting by T t (N|N ) ≥ 0 the transition probability per unit time from the state N to the state N at time t, then it is possible to write down, dropping the sub index 2 in the probability: where o(∆t) means that there are terms that goes to zero faster than ∆t, and therefore for small ∆t can be neglected. Because of Equations (41) and (42) in Equation (40), it is obtained: After some calculus, expanding P(N, t − ∆t|N 0 , t 0 ) and P(N , t − ∆t|N 0 , t 0 ) in Taylor series and retaining only the terms of order ∆t it is obtained: Equation (44) is known as the Master equation, the first term of this equation gives the increase in the probability per unit time due to the transitions from other states N to the state N at time t, while the second term gives the probability diminishing rate due to the transitions from the state N to other states N .
Let us discuss now the usual form of the Master Equation for chemical kinetics known as the chemical master equation (CME). This equation can be obtained directly from Equation (44), but it is more instructive to deduce it by a probability balance method.
Let us write first a probability balance for P(N, t + ∆t|N 0 , t 0 ) in terms of the probability at an earlier time t, and the propensities using the propensities defined in the previous section. Because of the propensities only depend on the system collective state N when the tank system it is well stirred and the properties like the temperature are uniforms. Then, we can write because of the propensity definition and the change ν (r) produced by each reaction type r on the number of molecules of each chemical compound: Passing P(N, t|N 0 , t 0 ) to the left hand side, dividing then both sides by ∆t, and computing the limit when ∆t → 0 , it is obtained, on account that lim ∆t→0 o(∆t)/∆t → 0 , the following equation for P(N, t), where we have omitted the N 0 , t 0 because we consider only those states that can be attained from the chemical system initial conditions: This is the Chemical Master Equation that expresses the probability evolution of P(N, t) in terms of the propensities for each network reaction type and the probabilities P N − ν (r) , t to be at the states N − ν (r) that can reach the state N in one reaction. The last term of this equation gives the probability to be removed from the state N by the different kinds of reactions. Now we consider the following identity that holds for the inverse reactions: Because of Equation (47), the CME (46) can be written in the following form: This is an alternative way of writing the CME because of the existing relations (47) between the direct and the inverse reactions.

Limitations of the Chemical Master Equation
To finish this section, we will discuss deeply the limitations of the CME formalism developed in the previous sections. As pointed out by Nicolis and Prigogine [4], the fact that the propensities depend on system collective variables implies that important fluctuation phenomena associated to local parameters of the system, as for instance local variations of the concentration, and correlation lengths over which two part of the system influence mutually are not considered in the approach of the CME. Also, Gillespie et al. [41], have discussed the conditions for validity of the CME in diffusion-limited systems, especially when bimolecular reactions occur, the two fundamental conditions which must hold are that the reactant molecules in the reacting system are at the same time "dilute" and "well-mixed". These authors have proved that the satisfaction of both conditions is necessary for the predictions of the CME to be true in a system where bimolecular reactions take place. Notice that the requirement for diluteness in a solution of reacting chemical compounds applies only to the reactant molecules (solute), and not to the chemically inert solvent molecules as pointed out by Gillespie et al. [41]. However spatial effects are present in many chemical and biological systems, and therefore for these systems do not hold the assumption of spatially uniform or well mixed distribution of chemical reactants. In his seminal paper, Turing [42,43] demonstrated, theoretically, that a system of reacting and diffusing chemicals compounds could spontaneously evolve to spatially heterogeneous patterns when departing from an initially uniform state in response to very small perturbations. He considered a very simple chemical system formed by two chemical compounds, in which one was an activator and the other one an inhibitor. That is, one chemical compound, the activator, stimulated and enhanced the production of the other compound, the inhibitor, which, in turn, depleted or inhibited the formation of the activator. Turing showed that if the diffusion of the inhibitor was bigger than that of the activator, then a diffusion-driven instability could result. This is a non-intuitive behavior, because we are used to think that diffusion is a homogenizing process. The first experimental evidence of Turing spatial patterns was observed in the chlorite-iodide-malonic acid (CIMA) reaction with starch [43]. More recently Maini, Painter and Chau [44], have reviewed the formation of spatial patterns in chemical and biological systems.
In this way spatial stochastic effects, which it is well known play important roles in many systems are neglected in the CME approach. The situation for these cases can be amended using the Reaction Diffusion Master Equation (RDME), in which space Ω is partitioned discretely into an arbitrary number of voxel (Volume Elements) each one with a volume ∆Ω [44,45]. Diffusion can then occur between different voxels, and reactions can occur within voxels on the assumption that reactants are well-mixed.
However, even in this case the description is not complete for many system because if the reactants are dissolved in a solvent, the motion of the solvent affects to the reactants, and we can have in addition to the normal diffusion process, the one known as turbulent diffusion produced by the small eddies and also it is possible to have transport due to the solvent motion.
It is important to know the factors which influence the spatial correlations to understand the limitations of the CME, consider for instance a simple example of non-equilibrium non-linear chemical kinetics to understand the correlation origin. Nicolis, Prigogine, Van Kampen, and Van Den Broeck among others started the study of this issue [4,8,46]. Consider the following chemical reaction A K (+) 1 → 2X, and its inverse 2X → A between particle pairs which are close enough this reaction depletes the system from pairs of X particles and produce A molecules. Because the particles X created in the same reaction are both produced at the same spatial point → r and time t, then these particles are obviously correlated. Then, due to diffusion this correlation propagates spatially until one of the two particles suffers a reaction, if this reaction takes place on average at a rate κ, then the characteristic correlation length, denoted by l c is given by the expression [8,46]: being D F , the diffusion coefficient of the chemical species X. Partitioning the space into equal volume cubic cells with volume ∆Ω, and denoting the centers of this cells by their position vectors The time evolution of this multivariate probability is produced by two different mechanisms the first one are the chemical reactions occurring inside the cells, and the second one is the diffusion of particles between adjacent cells [8,46]. The equal time correlation function g (X) → r , ∆Ω , evaluated at these points in the continuous limit when ∆Ω → 0 , as follows: The equation obeyed by the correlation function at two different points can be calculated by the usual procedure of multiplying the MCME by N The steady state solution of the equation obeyed by the equal time correlation function is given by the following expression [8,46]: where the characteristic correlation length l c is given by Equation (49), being κ the rate at which the X molecules are destroyed by the inverse reaction, this rate is obviously given for this simple example by: and J (1) X is the chemical thermodynamic flux for the X chemical compound which is given by: If the system is at equilibrium conditions then the flux J (1) X = 0, and there is not local creation of correlations. Also notice that according to Equation (49) increases then the correlation length diminishes. Therefore, correlations can be important when the chemical system is far from equilibrium and therefore spatial effects can become important under this situation.

The Method of the Probability Generating Function or Z-Transform
The method known as the probability generating function or simply generating function method, or z-transform method consist in defining the z-transform of the function P(N, t) by means of the following expression [47]: where z ≡ (z 1 , z 2 , . . . , z n s ), denotes the corresponding points in the transformed space. From the multivariate generating function G(z, t), is easy to obtain various kinds of moments. For instance, the k-th factorial-moment can be obtained easily from G(z, t) as follows: The second order cross-moment N i N j , can also be obtained easily from G(z, t): To obtain the evolution equation for the multivariate generating function we proceeds as follows, first we multiply the chemical master Equation (46), by z n s , and then we sum up for all the N values i.e., N 1 from 0 to ∞, N 2 from 0 to ∞, and so on, and finally, we substitute the propensities by their expressions given by Equation (37). This set of operations yields the following result: To obtain Equation (57), we have considered Equation (46), and we have defined K (+) r which denotes the K r constant for the direct chemical reaction r, while K (−) r denotes this constant for the inverse of the direct r-th reaction. Considering that all the components of N = N − ν (r) are positive, because physically only can take positive values, then it follows that the probabilities P N − ν (r) , t must be zero when some component of N = N − ν (r) is negative. Because the components of the vector N can take only values between 0 and ∞, then we can perform a change of variables in the summations from N to N . Finally, we rename again N as N. Notice than when performing these operations for a given reaction r we have that z , so finally, Equation (57) can be rewritten in the form: Finally, we can recast Equation (58) in a more simplified form: Equation (59) can be further simplified as follows, first we consider: From Equation (60) the following identity is obtained: The identity (60), when used in the Equation (59), yields the following result: is common to both terms in Equation (62), so this equation can be further simplified to yield: Equation (63) is the general partial differential equation (PDE) for the generating function of the chemical reaction networks. This equation is complex to solve in general due to its non-linear character but can be solved exactly for some simple cases.
To see that this equation is correct let us consider the chemical reaction aA + bB which is the example used by Smadbeck and Kaznessis [13]. For this example, the equation obeyed by the multiplicity generating function of this simple reaction because of Equation (63) is: Defining, 1 Ω 1 Ω a Ω b , it yields the same result that the one obtained by Smadbeck and Kaznessis.

The Method of the Expansion of the Chemical Master Equation in an Order Parameter
As we have discussed previously, the fluctuations in the number of molecules of a given chemical compound X j , are produced because of chemical reactions originate by collisions between the individual chemical compounds that are discrete in nature. Each chemical reaction produces one jump ν (r) , in the multidimensional vector space N n s , the components of this jump ν (r) j are very smalls compared with the total number of molecules N j of each chemical species X j . In addition, the macroscopic characteristics are due to the total population of each chemical compound. It is well known that for a population of N j molecules, the fluctuations are of the order of N j , and its effect on the macroscopic properties is of order N j/N j ∼ N −1/2 j . The election of the parameter to perform the expansion is usually the system volume Ω, because for large volumes the jumps are very small compared with the system population. There are different forms to perform this expansion, Van Kampen [8] proposed to split the dynamic of the molecular concentration in a deterministic part denoted by φ j and a fluctuating component η j , so that we can write down: However, Samoilov and Arkin [48] proved that for some class of biochemical networks the mean behavior differs significantly from the rate equation of chemical classical kinetics (CCK). For this reason, especially for some class of biochemical networks involving bimolecular reactions, where the propensities according to the Equations (33) and (36) depend non-linearly on the populations, Andreychenko et al. [49] propose to perform an expansion around the true mean: where from Equations (65) and (66) it is obtained that η j , is a random variable centered around zero and that quantifies the fluctuations around the true mean value.
To perform the systematic expansion of the Chemical Master Equation, first we write the CME (48) in the following form: Next, we express the propensities for the direct reactions as follows: and for the reverse reactions: where K (+) r = K r , and K (−) i = K n r +i . Therefore, the chemical master equation can be expressed because of Equations (67)-(69) in the form: Next, one defines the translation operators in the multidimensional vector space N n s , as follows: where: Using the operators defined in equations (71) and (72) in (70), the CME can be written down in the compact form: Next, we split the stochastic variable N j using Equation (65), in this approach the translation operators can be expanded in Taylor series [50]: Therefore, assuming that ν , on account that the ν (r) j are integer numbers and retaining only the first three terms of the E j expansion can be written as follows: where the multinomial coefficients [51] are defined in the standard way: Also, it is observed that the partial derivative with respect to t, is computed at N constant for all its components, therefore we can write: where to compute the time derivative dη j dt , we have considered that at N j = cte, and because of N j Ω = φ j + Ω −1/2 η j , then it is deduced that dη j dt = −Ω 1/2 dφ j dt . Therefore, the final form for the probability of the fluctuations is obtained because of Equations (76), (77) and (79) in Equation (73), which yields the following result: j , therefore for this case from Equation (73) it is deduced: Obviously, the expansion in Equation (80) can be performed at higher orders and in this case higher order multinomial coefficients will appear in (80) and (81).
If we use an expansion around the true mean as in Equation (66), then the partial derivative with respect to t, is computed at N constant for all its components, therefore we can write: Also, in this approach the translation operators can be expanded in Taylor series: Then, if ν (r) j > 0, it is easy to prove that the probability Π(η, t) of the fluctuations around the mean obeys the following equation: A similar equation to (84) can be also deduced for the case ν (r) j < 0.

Methods to Obtain the Probability Distribution from the Moments
The CME has an analytical solution in a few cases. The handbook of stochastic methods by Gardiner provides several examples of exactly solvable models of the CME [40]. More examples can be found in references [52,53]. In the rest of cases the solution of the CME can be found by Monte-Carlo methods, by the finite state projection method [12], or reconstructing the PDF from the moments as it is shown below.

Reconstructing the Probability Distribution from the Moments for One Component Case
For this case, we assume for simplicity a single component, and that at some time instant t, the moments µ k k = 1, 2, . . . , N M, of the distribution are known, where NM is the number of moments. If P(N, t), denotes the probability to have N molecules at time t, then we can write the following set of conditions: µ k (t) = ∑ N N k P(N, t), k = 1, 2, . . . , N M, In addition, we have the normalization condition for the probability, so that Then following the maximum entropy principle [29,30,49], the probability distribution function P(N, t), which must be chosen, is the one that maximizes the Shannon information entropy and at the same time satisfies the restrictions (85) imposed by the moments and the probability normalization condition (86). This problem is a typical optimization problem with constraints that can be solved using the Lagrange multipliers method by building the functional: This optimization problem leads to the following solution, known as the maximum-entropy PDF: The Lagrange multipliers are obtained from the restriction conditions (85) and (86). Applying the normalization condition (86), to Equation (88), it is obtained: so, the probability distribution that maximizes the information entropy is given by the expression: where the unknown Z depends on the Lagrange multipliers. There exist different algorithms to obtain the Lagrange multipliers from the moments as explained in Andreychenko et al. and Woodbury [49,54]. The previous deduction can be generalized to the case of the moment evolution with time, in this case the time can be considered as a parameter, and therefore we can use variational calculus as in Caticha and Preuss [55]. Because we must maximize the entropy at each time instant then we look for the time function P (N, t), which makes maximum the information entropy along time and at the same time satisfy the restrictions imposed by the moments, notice that for this case the Lagrange multipliers should depend on time, so we write the functional in the form: Then, we perform the variations of the PDF and the Lagrange multipliers that are now function of the parameter i.e., P(N, t) + δP(N, t), and λ j (t) + δλ j (t), and we look for the probability distribution that makes δF = 0. Where the first variation is computed using standard rules of variational calculus as follows: Performing the operation (92) in the functional (91) one arrives to the following equation: Due to the arbitrary character of the functions δP(N, t) and δλ j (t), j = 0, 1, . . . , N M, one arrives to: where now the Lagrange multipliers are functions of t, also it is obtained the set of restrictions (85), and (86), with the probability P(N, t), given by (94).
We can use another moment basis, which is known as the factorial moment basis; in this case for one component we define the following set of population functions: Then, the factorial moments of order q denoted by {N q } f are defined by the following set of equations: If one knows the factorial moments up to a given order M, at time t, then using the MEP one can obtain the probability distribution that maximizes the information entropy and at the same time satisfy the restrictions imposed by the moments. In this case it is obtained the following result: The Lagrange multipliers λ q are obtained solving the system of equations:

Results of Applications to Chemical Reactions Networks
In this section we apply the methods explained previously to some specific chemical reaction networks. Beginning with very well-known cases and finishing with complex ones. Also, we study in this section the consistency of the results obtained by the expansion of the CME in an order parameter, by comparison with the results obtained using the PGF method. Finally, because of its importance in practical applications we have included a section devoted to the different methods to solve the closure problem.

Case of One Single Component in an Open System
One of the simplest chemical reaction networks, which contain non-linear terms, is the following one: where it is assumed that the reaction takes place in a well stirred tank, being the compound A supplied continuously, so the amount of A remains constant and the chemical compound B is also removed continuously so B is also assumed constant. In this example we denote the reaction (100) as reaction 1 and the reaction (101) as the reaction (2). For these reactions ∑ we have that the stochiometric coefficients and the jumps are: Then Equation (80) for the network formed by the direct chemical reactions (101) and (102) can be written down as follows: At order Ω 1/2 , i.e., collecting all the terms proportional to Ω 1/2 , and which are the biggest ones, we have the following equation: From Equation (105) it is deduced the macroscopic equation of chemical kinetic for the macroscopic concentration φ x for the open systems chemical network (100) and (101): We observe that in Van-Kampen book [8] the coefficient of the term φ 2 x is 2. The reason of this is that K (+) 2 = 2K, because the propensities were defined in Section 2 per pair of interacting particles, and the number of couples when the particles are identical is divided by 2, because it is given by the combinatorial number N 2 .
The next step is to collect all the terms of order Ω 0 , which after some simple algebra yields: From, Equation (107) at this order Ω 0 it is obtained the following equation for the probability Π(η x , t): According to Risken's book [56], Equation (108) is a linear Fokker-Plank partial differential equation with time dependent coefficients through φ x (t). This equation is the same one obtained by Van Kampen [8] at the same order. The solution of Equation (108) leads to a Gaussian distribution as shown by Van Kampen [8]. Therefore, we only need to obtain the first two moments of the distribution to know its PDF: where D is the distribution support. Obviously, if we go to higher orders of the expansion, the equation that governs the evolution of the probability is not Equation (108) and will contain additional terms and higher derivatives. The next step is to obtain the equations that govern the evolution equation for the first and second moments, in this case from Equation (108), because of (109) it is obtained: Let us study now the reaction network formed by (100) and (101) by the method of the moment generating function, in this case we assume also that the population of A is constant, and that B is also removed continuously to maintain B also constant. For this simple case the generating function evolution equation of this reaction network is given because of Equations (63), (102) and (103) by: Application to Equation (112)  yields for the first moment N x the following result: Defining the variation δN x of the number of molecules of the chemical compound X respect to the mean value as δN x = N x − N X , it is obtained: From (113) and (114), it is deduced that: Defining an average concentration of X as φ x = N x /Ω, it is obtained after dividing Equation (115) by the system volume Ω: Only if the probability distribution follows a Poisson PDF, we have that σ 2 = N x , and in this case Equation (116) simplifies to: But in general, this is not true, so that in the equation for the mean average appear terms that depend on the fluctuations. In addition, as we will see below, the evolution equation for any order moment depends on higher order moments.
Application of the operator ∂ 2 to the equation (112) of the moment generating function yields after some calculus: Equation (118) tell us that the evolution equation for the second factorial moments depends on the first, second and third factorial moment. Let us check now the consistency of this result, to do that let us assume that the PDF follows a Poisson distribution. Then as it is well known the factorial moments N q x f of a Poisson PDF verify the following conditions, see the book of Gardiner [40]: Then, because of (119), the right-hand side of Equation (118) yields if the PDF is a Poisson one: Then the left-hand side of Equation (118), yields: Apparently, the right-hand side of Equation (120) (120) and (121) by Ω 2 , and expressing the magnitudes in terms of average macroscopic concentrations it is obtained from (120) and (121) the following results: So, it is deduced that the difference between both equations is the term −K (+) 2 Ω −1 φ x 2 . It is observed that this term is much smaller than the other two terms because is divided by the order parameter Ω that is the system volume, that is generally large, and makes this term very small compared with the other two terms so the apparently contradiction for poissonian systems is solved.

Examples of Two Component Chemical Reaction Networks
One simple example of two component chemical reaction networks in a closed system with non-linear effect is the one studied by Smadbeck and Kaznessis [13] and consist in the reversible dimerization of two monomers. In this case we have only one reaction and its inverse, i.e., the reaction of two monomers A to form a dimer B and at the same time the decomposition of this dimer B in the two original monomers A. Each reaction has its own constant, so we write: For this case the stoichiometric coefficients and the jumps are: The next step is to apply Equation (63) to this chemical reaction with the stoichiometric coefficients (125), this gives the following equation for the evolution of the generating function G(z A , z B ; t): From Equation (126) one can deduce the equations obeyed by the factorial-moments, the q-th factorial moment of a single component is obtained by the following expression, as can be checked easily: Also, the bivariate factorial moments N q A N s B f can be obtained from the generating function G(z A , z B ; t) by the expression: In general application of the operator ∂ q+s to both sides of Equation (126) lead to the evolution equation obeyed by the factorial moment N q A N s B f . It is better to gain efficiency as recommended by Smadbeck and Kadnessis [13], to apply the derivative operator successively and systematically to generate the equations obeyed by all the factorial moments up to a given order. This can be performed with the symbolic operational algebra and calculus tools of MATLAB or MATHEMATICA. For instance, applying the operators ∂ ∂z A to Equation (126) it is obtained the equation: Setting z A = z B = 1, in Equation (129) gives the equation for the evolution of the first factorial moment in terms of higher order factorial moments: Applying again the operator ∂ ∂z A , to Equation (129) yields the equation: Setting again z A = z B = 1, in Equation (131) yields the evolution equation for the second factorial moment N 2 A f : The idea now is to define an ordered factorial moment column vector {N 1 N 2 . . . .N n s } f that contains the lower order moments in the lower index components and the higher order moments in the higher index components. Then, for the case of two components this column vector is defined as the transposed of the following row vector: where the upper script T denotes the transpose of the row vector to transform it in a column vector.
Where in practical application we need to truncate the number of components of this vector at some finite number NC. In general, one obtains an equation for the truncated factorial moment vector. This equation for the evolution of the vector {N 1 N 2 . . . .N n s } f was first deduced by Smadbeck and Kaznessis [13], and for the case of the reaction (124), the evolution equation for the column vector of ordered factorial moments, takes de form: where the matrix M f , which gives the evolution of the factorial moments for the dimerization reaction and its inverse in a closed system has the components displayed in Table 1. Notice that for the zero order factorial moments, the components are 0. The advantage of using the factorial moment expansion is that the matrix M f is sparser than if one uses an evolution equation for the polynomial moments [13,14]. Also, the bandwidth structure of the matrix for the factorial moments has fewer bands than for the polynomial case. Due to its importance in practical and numerical calculation we discuss in the last section the closure problem. Table 1. Components of the matrix M f , which determines the factorial moments evolution of the dimerization reaction network.

The Michaelis-Menten Chemical Reaction Network
In this section we study the Michaelis-Menten chemical network, that has been previously studied by several authors [13,57,58] and which is more complex than the previous examples. In this reaction network we have three reactions, in the first one a substrate (S) reacts with an enzyme (E), to form a complex denoted as (E:S). In the second one, this complex degrades back to its original components through the inverse reaction. In the third reaction the complex (E:S) is transformed in a product (P) plus the enzyme (E), so we can write: This network has apparently four components and three reactions (reaction 1 plus its reverse reaction and reaction 2). For this network the stoichiometric and the jumps coefficients are: The next step is to apply Equation (63) to this chemical reaction network with the stoichiometric coefficients (137), and (138), this calculation gives that the generating function G z S , z E , z E:S , z p ; t , obeys the equation: From Equation (139) it is possible to obtain the equations obeyed by the factorial moments, , . . . . to Equation (139), and then setting z E = z S = z E:S = z P = 1, proceeding in this way, we have obtained the following results for the evolution equations of the factorial moments: Next, because of the following relations that express the conservation in the total molecule number of substrate and enzyme in their different forms, we can write: Then, we substitute in Equations (140), (141) and (142), N E:S by N E T − N E , and after a little algebra the following results are obtained: where we recall that N 2 . Also, we have obtained the evolution equations for the factorial moments {N Sf } and N 2 S f not displayed here for simplicity. Therefore, the matrix M f , that gives the evolution of the factorial moments for the Michaelis-Menten reaction network has the components displayed at Table 2. Table 2. Components of the matrix M f , which determines the factorial moments evolution of the Michaelis-Menten reaction network.
We must remark that all the terms of the matrix M f , which determine the moments evolution for the Michaelis-Menten network are exactly the same that the ones obtained by Smadbeck and Kaznessis [12] in reference [13], and −2K (−) 1 in this paper.

The Closure Problem
In this section we discuss the closure problem that consist in approximating the high order moments, which appear in the equations for the moment evolution up to a given order M, in terms of lower order moments, this problem has been studied by different author as Ruess et allied [25], Gillespie [15], Grima [16], and many others using different approaches and by by Smalbeck in his PhD thesis [59], using the zero information closure, because this last method is based on the maximum entropy principle will be the first one to be explained. In addition, Ruess et al. have performed a systematic comparison of different closure methods [25], so this paper it is a good starting point for studying this subject.

The Zero-Information Closure
In this subsection we use one component system to simplify the explanation. The closure problem consists in finding a function F to express the higher order moments µ i>M in term of lower ones, so it is possible to write: where if one uses the polynomial moments, which are defined by the expression: and, if we know all the polynomial moments µ i (t) of order i ≤ M at a given time t, then we can build a functional with the information entropy, and the known moments up to order M as restrictions as was done in Section 4. In this case we can obtain an approximate PDF that maximizes the information entropy with the M restrictions imposed by the known moments at time t. If we denote this PDF based in the known M polynomial moments with the subscript S, then we may write: The Lagrange multipliers λ k (t), are obtained solving the following set of equations for the previously known moments µ i (t): Once we obtain the Lagrange multipliers, we know an approximate expression for the PDF, which is obtained from the first M-th moments, and that we denote P S (N, t) and which is a function of µ 1 , µ 2 , . . . , µ M . Then, the higher order moments can be expressed in term of lower order ones by means of the expression: The estimation of the higher moments in this way is denoted by the superscript S, that means that the MEP has been used, this closure scheme is known as the zero-information (ZI) closure, and has been recently used by Vlysidis and Kaznessis, using factorial moments for the Slogi chemical reaction network [60]. It does not matter to use one polynomial basis or the factorial moment basis because both polynomial and factorial basis are related by a transformation matrix.

Different closure methods
In this section we discuss some additional closure methods that have been used in chemical reaction networks. We perform this discussion for a multivariate chemical system with n s different chemical species, as we have seen in previous sections. Fist we define the so called ordinary or normal moments of order k for the population of the different chemical species as follows: where we have adopted the convention that the indexes are ordered from the smaller ones to the larger ones, i 1 ≤ i 2 ≤ . . . ≤ i k , and that any index can be repeated, for instance µ 1,1,2 = N 2 1 N 2 . The second type of moments that will be used are the central moments, these moments are defined as: and the convention index adopted is the same one that for the ordinary moments, which some authors denote as raw moments. The third type of moments introduced are the cumulants [25,40], which are defined as follows: where G κ (s 1, s 2,..., s n s ) is the cumulant generating function, which is defined by the following expression [19,40]: For instance, to show how Equation (154) works, one can use this expression to obtain the expressions for some low order cumulants in terms of ordinary or normal moments: From the previous equations, it is easily deduced that we can express normal or raw moments in terms of cumulants for instance it is immediately deduced from Equations (156) and (158) that: With these previous definitions it is easy to stablish several of the most popular closure schemes found in the literature for chemical reaction networks. The first one is the central moment neglect (CMN) moment closure approximation (MCA), which will be denoted by CMN-MCA. If this closure scheme is performed at order M, then all the central moments of order larger than M, are neglected. Thus, we can write: The second closure scheme, which it is perhaps the most popular one is the "cumulant-neglect MCA". In this closure scheme if the selected order is M, then all cumulant moments of order larger than M are neglected. Thus, the following expression can be written: For instance, in the cumulant neglect closure of order 2, all cumulants of order larger than two are set to zero, this immediately leads to the result that it is possible to find expressions for the third and higher order normal moments in terms of first and second order ones [19,40], for instance one finds for a two component system that: and similar equations are obtained for the rest of third order moments. The next step to closure the problem is to substitute the thirds order moments in the moment evolution equations, by their expression in terms of lower order moments and in this way, we close the hierarchy of moment equations. The "cumulant-neglect MCA" has been used previously in the field of biochemical kinetics by different authors as Gomez-Uribe and Verghese [61], and also by Schnoerr et al. [19].
There are many more closure methods as the Poisson MCA, the log-normal MCA and so on, also there are several software packages as MOCA, StochDynTools, MomentClosure, which contain several methods to close the moments evolution equations [19,61].

Study of the Consistency of the Development of the CME in an Order Parameter
In this section we are going to check if the results obtained using the expansion of the CME in an order parameter are consistent with the results obtained using the probability generating function method.
We start with the chemical reaction network formed in an open system by reactions (100) . It is obtained after some calculus that the first and second polynomial moments obey the equations: The next step is to substitute N x = Ωφ x + Ω 1/2 η x , as performed by Van-Kampen in Equations (163) and (164) or N x = N x + Ω 1/2 η x as performed by Andreychenko et al. [30].
Performing first the substitution N x = Ωφ x + Ω 1/2 η x in (163) and collecting all the terms of order Ω 0 , yields the following macroscopic kinetic equation: Collecting now the terms of order Ω 1/2 , gives the result: Equations (165) and (166) are the same ones that the ones obtained from Van Kampen expansion of the Master Equation. If we perform now the substitution N x = Ωφ x + Ω 1/2 η x , in (164) and we use the following identities: Because of Equations (167) and (168), the following result is obtained from Equation (164), collecting the terms of order Ω: The terms in this case are the same ones that in Van Kampen approach but now some of the coefficients are different as it is deduced comparing Equations (169) and (111).
Let us use now, the alternative recently proposed by Andrechenko et al. [49], in this case they propose to use N x = Ωφ x + Ω 1/2 η x with φ x = N x Ω . The great advantage of this formulation is that η x = 0, as was mentioned previously in Section 5.1. Because of this property we can write: Performing first the substitution N x = Ωφ x + Ω 1/2 η x in (163) and collecting all the terms of order Ω 0 , yields the following macroscopic kinetic equation: Collecting the terms of order Ω 1/2 gives: In this case, because of (170) and (171), the following result, is obtained from equation (164) collecting the terms of order Ω: We observe that the macroscopic Equations (165) and (172) are the same ones that the ones obtained by Van-Kampen method of expansion of the master equation. Also, Equation (166) for the first moment of the fluctuations when we use N x = Ωφ x + Ω 1/2 η x , is also the same. But for the second moment evolution equation of the fluctuations, Equation (174), the terms are the same in Andreychenko approach N x = N x + Ω 1/2 η x , but some coefficients are different if N x = Ωφ x + Ω 1/2 η x , as it is observed in Equation (169) for the second moment of the fluctuations.

Discussion and Conclusions
The application of methods based on the theory of stochastic processes has become very popular in chemistry and biological systems, especially since the development of the Chemical Master Equation (CME), and the Monte-Carlo methods to solve properly the CME by sampling over the probability distributions [4,[6][7][8]. In addition, the introduction of the finite state projection method (FSP) by Munsky and Khammash, which provides a direct solution of the CME at a desired degree of approximation without generating Monte-Carlo simulations, has supplied an alternative way to Monte-Carlo methods [12]. Since the pioneering works of Van Kampen [8] and Nicolis and Prigogine [4], the developed methods and the range of applications of the CME have experienced a noteworthy increase. This is especially true for the methods based on Monte-Carlo as the Sampling Stochastic Approach (SSA) of Gillespie, which performs a sampling over the underlying probability distributions. It is necessary to mention the methods based on the Chemical Fokker-Planck (CFP) equation, which can be obtained performing a Kramers-Moyal expansion in the Chemical Master equation and truncating this expansion to second order, in this way one obtains the CFP Equation.
However as pointed out by Gardiner (p. 267, [40]) this is only an approximation whose asymptotic expansion is identical to order Ω −1 with that of the CME.
More recently people as Smadbeck, Kaznessis, Sotiropoulos, Vlysidis and others [13,14,59] have tried to solve the moment equations and to reconstruct the PDF from the moments applying the maximum entropy principle. The problem of this approach is that the moment order must be truncated at some value, so one needs a method to compute and truncate in a consistent way the hierarchy of moments. Because if not the number of differential equations to solve becomes too large, and many moments are necessary to obtain a good approximation to the probability distribution. There have been many attempts in the past to find a consistent methodology to truncate the moment hierarchy; some of these methods as the "cumulant-neglect MCA" are reviewed in this paper [20,40,62]. At this point the zero-information truncation or ZI truncation provides a methodology to obtain in a consistent way the truncated moments in terms of lower ones; this method gives similar results than the SSA method with the improvements to this method performed by Gibson and Bruck [63].
In this review we have tried to present a unified vision of the chemical reaction networks and the existing methods to compute the probability distribution of the different chemical compounds that enter in the network. We have started by the root basis of chemical kinetics from a microscopic point of view to compute properly the reaction rates of the different chemical reactions that are necessary to obtain the chemical propensities that enter the chemical master equation as a fundamental cornerstone. We have centered the analysis in the bimolecular reactions computing from a microscopic basis the reaction rates, that are necessary to obtain the chemical propensities.
The next step has been to obtain the chemical master equation (CME), traditionally this equation is deduced from the Chapman-Kolmogorov equation, and we have used this way to obtain the master equation. In addition, we have considered very illustrative to deduce the chemical master equation from a probability balance that capture better the physical roots of the problem. Then we have discussed the limitations of the CME, first we have reviewed the two fundamental conditions, which must hold to be valid the application of the CME in diffusion-limited systems, especially when bimolecular reactions occurs. These conditions were studied by Gillespie et al. [41], and are that the reactant molecules in the reacting system are at the same time "dilute" and "well-mixed". In addition, it is important to know the factors which influence the spatial correlations to understand the limitations of the CME; this influence is studied through the two point correlation function in Section 2.5, it is found that correlations can be important when the chemical system is far from equilibrium and in addition the correlation length is large.
Then, we present two alternative methods that are used to obtain the probability distribution of the chemical species of the reaction network, one is the expansion of the chemical master equation in an order parameter that is the system volume and the other method is the probability generating function method. To see if both methods are equivalent we have performed a consistency analysis of both methods in Section 5.5, and we have found that although the results of both methods for the macroscopic evolution equations and for the first moment of the fluctuations are similar, some differences can exists in the coefficients of the second moment equation depending on the way that is performed the splitting of the random variables.
Then, we have deduced the expansion of the chemical master equation in the order parameter using a multinomial expansion of the jump operators, E ν (r) j j . The equation obtained gives the same results than Van-Kampen expansion method of the Master Equation for some simple reaction network studied in Sections 5.1 and 5.2.
Because we can compute the moments of the distribution of the chemical species then it is important to reconstruct the PDF from these moments. The method to obtain the unknown PDF from the moments is based on the application of the maximum information entropy principle (MEP). The main drawback of this methods is that requires to solve a set of non-linear equations to obtain the values of the Lagrange multipliers that are present in the exponent of the PDF expression. Three types of basis can be used to reconstruct the PDF, the polynomial basis, the factorial moment basis and the cumulant basis. The factorial moments evolution equations can be obtained directly from the generating function evolution equation as was shown in Sections 5.1-5.3. Performing some algebra is possible to obtain the polynomial moments from the factorial moments. The polynomial moments and the factorial moments as shown by Smadbeck and Kaznessis [13] are both related by a similarity transform matrix. Other important basis, is the cumulant basis, the interest of this basis is that the cumulants have cluster properties, in the sense that for independent random variables all the cumulant vanish, notice that the factorial moments do not enjoy this property [40]. Also, it is possible to express the polynomial moments in terms of cumulants. So, we can pass from one basis to the other.
The expansion of the master equation in an order parameter and the multiplicity generating function method are applied to several simple univariate and multivariate chemical reaction networks. First the multinomial expansion formula derived in Section 3, is applied to a non-linear simple reaction network and the resulting equation for the evolution of the macroscopic variable and the moments of the fluctuation are compared with Van-Kampen results [8]. It is found that the results are the same. Then this same case is studied by the probability generating function method (PGFM), by application of this method, see also Section 5.1, we obtained the equations obeyed by the factorial moments, at this point a study of the self-consistency of the PGFM is performed in the equations obeyed by the factorial moments when the system is Poissonian and the method is self-consistent as explained at the end of Section 5.1. The next step was to check if the two previous methods, the PGFM and the expansion in an order parameter of the CME, were cross-consistent i.e., give the same results. This study was performed in Section 5.5. Instead to perform a systematic expansion in an order parameter we obtain first the equations for the evolution of the factorial moments, of first order, second order and so on. And at posteriori, we perform the substitution N x = Ωφ x + Ω 1/2 η x , as performed by Van-Kampen in Equations (163) and (164) or N x = ΩN x + Ω 1/2 η x as performed by Andreychenko et al [49]. Then we collect the terms of the same order in the expansion, and the result yields the correct macroscopic equation, the same equation for the evolution of the first moment of the fluctuations, but for the second moment η 2 x , the equation have the same terms than in Van-Kampen results for this same case [8], but two of the terms has a different coefficient in the approach N x = Ωφ x + Ω 1/2 η x . However, the equation for the second moment was the same in Andreychenko et al. approach [49]. This cross checking shows that both methods are not completely equivalents, and that further research is necessary.
Then in Section 5.2, we apply the PGFM to a two-component chemical reaction network studied by Smadbeck, and Kaznessis [12], that consisted in the reversible dimerization of two monomers. Also, in Section 5.3 we have applied the PGFM method to the Michaelis-Menten reaction network. We used the factorial moment basis for both cases, and the results were the same ones that the ones obtained by these authors, except one small difference in one term for the Michaelis-Menten network. Finally, in Section 5.4 we have reviewed the closure problem studying some classical closure methods as the "cumulant neglect MCA" and the recent advances in the ZI closure method, as a systematic and consistent approach to perform the closure of the moment hierarchy at the desired truncation order.