Mathematical explanation of the predictive power of the X-level approach reaction noise estimator method

The X-level Approach Reaction Noise Estimator (XARNES) method has been developed previously to study reaction noise in well mixed reaction volumes. The method is a typical moment closure method and it works by closing the infinite hierarchy of equations that describe moments of the particle number distribution function. This is done by using correlation forms which describe correlation effects in a strict mathematical way. The variable X is used to specify which correlation effects (forms) are included in the description. Previously, it was argued, in a rather informal way, that the method should work well in situations where the particle number distribution function is Poisson-like. Numerical tests confirmed this. It was shown that the predictive power of the method increases, i.e. the agreement between the theory and simulations improves, if X is increased. In here, these features of the method are explained by using rigorous mathematical reasoning. Three derivative matching theoremsare proven which show that the observed numerical behavior is generic to the method.


Introduction
Noise is an integral part of the workings of the living cell biochemistry [1]. There are many types of noise and this work focuses on the intrinsic noise. If reactant copy numbers are low they can fluctuate widely [2]. These fluctuations can severely influence the dynamics of the cell and need to be carefully controlled [3].
Describing intrinsic noise has attracted a lot of effort. A range of theoretical methods have been developed to study intrinsic noise. However, an accurate characterization of the reaction noise is not easy. A direct solution of the chemical master equation for the system is often not possible since the number of configurations can be exponentially large. Numerical simulation methods can be used to avoid this problem, and are often implemented by using the Gillespie algorithm [4]. However, to obtain accurate prediction for moments of the particle number distribution function, e.g. the variance, sampling with a relatively large number of runs (simulations) is needed. This becomes impractical if the number of particle types is very large. A range of methods have been suggested to complement or replace these techniques. The focus of this work is on moment closure techniques [5][6][7][8][9][10][11][12][13][14][15].
The main idea behind moment closure approaches is to construct the equation system that can describe various moments of the particle number distribution function. In such a way there is no need to directly solve the master equation or perform a largenumber of computer simulations. The problem is that the equation system that describes these moments is, in principle, infinite. The main issue is to cut (or close) the hierarchy. This is done in two ways. First, one can try to make some specific assumptions about the reacting system which can be used to express higher order moments in terms of lower order ones. This procedure defines the moment closure function for the problem. Second possibility is to take the distribution function centered approach and assume that the particle number distribution function has a welldefined form, parameterized by a finite set of parameters. Variational calculus can be used to obtain the parameters [7,16]. In both cases the complicated many body dynamics is reduced to a set of ordinary differential equations that govern quantities of interest.
Two moment closure methods that were suggested previously will be of a particular interest. The PARNES method [14,15] is based on the assumption that pair effects dominate the dynamics. The method has been generalized into the XARNES method [17] so that higher order correlation effects can be included in the description. In fact, various choices for X result in a series of methods, e.g. X = P (pair effects), T (triplets), Q (quadruples), etc. Thus the PARNES method is the special case of the XARNES method with X = P. Both methods are based on the rather generic formalism of correlation forms which is used in statistical physics to model spatially extended many body effects. Each correlation form describes a particular correlation effect (single, pair, triple, and so forth). Correlation forms are used to perform a cluster like expansion of relevant quantities of interest. In the previous studies, the original correlation form formalism [18,19] used to model spatially extended diffusion controlled reactions has been adapted for describing well mixed reaction volumes.
A moment closure method works well in the intended domain of application, but it might easily fail if used elsewhere. Thus for a moment closure method to be useful it is necessary to, if possible, specify precisely in which situations the method is expected to work well. In [17] it was argued that the XARNES method should describe well systems with a Poisson-like particle number distribution function. This was confirmed by numerical studies. In here, it will be shown rigorously that such behavior is generic to the method. Also, there are some ambiguities while adapting the correlation form formalism to the well mixed situation, and the procedure is not unique. The problem is that a large number of spatial degrees offreedom need to be projected onto a much small number of variables, and there are many ways how this can be done. In here, it will be shown that the choice made in the previous studies is in some sense optimal.
To show all of the above the derivative matching procedure introduced in [9-11] will be used. In particular, the procedure from the technical report [9] willbe closely followed. This report was formalized in much shorter form in [10]. A very general multiplicative ansatz was suggested for the moment closure function. The precise form of the ansatz was found by using the derivative matching procedure. The function was parameterized by a finite set of parameters which were found by matching time derivatives of exact and approximate moments. This was done for the system in the pure state. The condition for a good match was derived in the form of a generic formula that involves the moment closure function parameters.
The same formula was obtained in a different context where the XARNES method was suggested [17]. This fact motivated the present work. The derivative matching procedure will be applied to investigate the accuracy of the XARNES method. There are some important differences in the setups in [10,11] and the present work. In here, the original procedure from [10,11] will be carried in the opposite direction. The XARNES method is based on the well-defined multiplicative ansatz for factorial moments. It will be shown that the ansatz implies good derivative matching properties. Also, while the original work focused on the pure state in here the Poisson state will be of interest.

The mathematical setup
A reaction system is defined as follows. It is assumed that reacting particles mix well and that particle positions are irrelevant. A configuration of the system can be specified by listing how many particles of each type there are in the system n = (n 1 , n 2 , · · · , n i , · · · , n T ) where variables n i , i = 1, 2,..., T, are positive definite integers used to denote the particle numbers. A configuration of the system changes in time due to the presence of reactions.
The full list of reactions is given by r 1 , r 2 ,..., r R and a reaction r α is formally defined in the usual chemical notation as The positive definite vectors with α = 1,2,..., R contain the stoichiometry coefficients for the reactions. It is assumed that the dynamics can be modeled as the Markov process where λ α is the reaction rate for a reaction r a having the unit ofinverse time. The dynamics defined in such a way is stochastic.
One possible way to describe the dynamics and characterize noise is to construct the master equation, solve it, and obtain the particle number distribution function P( n, t) which describes all stochastic properties of the system. However, both the computation and the direct inspection of the particle number distribution function is often tedious. It is more useful to investigate certain properties of the distribution function.
In practice this is done by computing various observables, or ensemble averages, as where f is the function that suitably parameterizes an observable of interest. A typical observable could be some low level moment such as the mean or the variance. Clearly, the direct use of (5) is not practical for large systems and one needs to avoid this route somehow. The idea is to project the details of the dynamics and obtain a coarse grained description by monitoring a selected set of observables instead of all configurations. These observables will be classified, i.e. labeled, in a precise mathematical way. For a fixed number of particle types, it is useful to introduce the vector space of all admissible labels. Formally, this is done through the following set of definitions.
Definition 1. The space of vectors that are used to label various correlation effects or, depending on the context, observables, will be denoted by Ω, where a variable m i is a positive definite integer. For example, the average number of particles of a type i will be labeled as e i = (0, 0, · · · , 0, 1, 0, · · · , 0) where the digit 1 appears on the i-the place.
Definition 2. It is useful to introduce the order, or norm, of a vector m ∈ as the sum of its components, This order will be used to classify various correlation effects. It is clear from the definition that The full set of observables will be split into two disjoint subsets. The first subset contains the observables that are being included in the projection and the second set contains the observables that are omitted. The observables in the second set need to be expressed somehow in terms of the observables in the first set. The following definitions will be used to classify such sets.
Definition 3. The set of vectors m ∈ that contains all vectors with orders 0,1,2,...,ξ will be denoted by Ω ξ , where ξ ≥ 1 is an arbitrary integer used to classify various theories. Definition 4. In a similar way, will denote the set of vectors m ∈ that are smaller than a vector m in the lexical sense defined above.
The set Ω ξ will be used to denote the set of observables being explicitly considered in the theory. For example, the mean field theory (the classical chemical kinetics) that neglect effects of fluctuation works with the first order effects where ξ = 1 (X = Singles). Thus the mean field theory would work by constructing equations for all ρ ei with e i∈ 1 . Sets of the type m will be used to limit various sums.
Definition 5. Finally, the set of all vectors that are not in Ω ξ will be denoted as, Further, several definitions will prove useful that generalize well known operations on numbers to similar operations in Ω. These generalizations greatly compactify some of the mathematical expressions that will be discussed. Definition 6. Let ω 1 and ω 2 be two arbitrary tuples of real numbers with rank T, such as ω = (ω 1 , ω 2 , · · · , ω T ) . The set of all such vectors will be denoted by R T . For any pair of such vectors the generalized power will be defined as Obviously, Ω⊂R T and this definition can be also used for vectors in Ω. It is clear from the definition that Definition 7. The binomial like coefficient that involves a pair of vectors m 1 and m 2 from Ω is formally defined as and it is assumed that a binomial coefficient The factorial-like symbol applied to a vector m from Ω is generalized as Definition 9. The product between two real numbers is generalized as where ω 1 and ω 2 are two arbitrary vectors from R T . Please note that the result of the operation is an element in the same set, ω 1 ω 2 ∈ R T . Also, the following identity related to this definition will be useful later on,

Exact equations of motion
There are several types of observables one could choose to work with. In this work factorial moments will be used. A factorial moment will be labeled by the related positive definite vector m = (m 1 , m 2 , · · · , m T ) ∈ and is defined as It was shown in [17] that the equation system for the exact factorial moments is given by and the structure of the equations will be explained in the following. Γ coefficients in the sum are given by for all m ∈ . A vector c ∈ R is a positive definite vector to be referred to as a contraction vector. Contraction vectors emerge during the field theoretic derivation of the equations of motion when one uses the Wick theorem. More details regarding the field theoretic setup can be found in [14,15]. One can derive the same set of equations in another way (not shown). Sums over contraction vectors will be of the central interest in the following.
The space of all possible contractions is defined by reactions that can occur in the system.
From the definition of Γ coefficients in (20) one can see that for a fixed u α or v α the sum over the contraction vectors in (19) is restricted. This is suggestive of the following formal definition of the space Ω R ⊂Ω: Despite the fact that the sum over contraction vectors is finite, the equation system for exact moments represents, in principle, an infinite hierarchy of equations since on the right hand side of the equation for a given ρ m there are terms that involve higher order moments since it may happen that m − c + u α > m . The hierarchy of equations for the exact moments appears not particularly useful. However, it can be used to devise approximation schemes. If higher order moments can be expressed in terms of few lower order ones then the equation system closes down.

A way of closing the hierarchy: the XARNES method
The XARNES method is based on the closure ansatz which is constructed by using the concept of correlation forms [14,15,17]: where wm denotes the correlation form labeled by a vectorm ∈ . The detailed motivation behind this equation is given in [17]. It is assumed that that correlation forms with the orders above a given threshold ξ can be assumed small, is not exactly equal to the related moment ρ m(t) , two different symbols, v and r, are used to denote them. However, if the approximation above works well, possibly when the number of vectors in Ω ξ becomes large, their values should be close.
Please note that if a factorial moment is zero, the left hand side of Eq. (22) becomes infinite owing to the singularity of the logarithmic function. This implies that the related correlation form becomes also infinite. Clearly, by construction, the XARNES ansatz is somewhat ambiguous for cases when some factorial moments vanish. In what follows it will be assumed that all factorial moments are strictly larger than zero but they can be arbitrary close to zero.
It was shown previously [14,15,17] that the assumption (22) is an implicit ansatz for defining the function that expresses a higher order non-base factorial moments νm (t) withm ∈¯ ξ in terms of the base moments ν(t) ≡ (· · · , νm (t),··· ) withm ∈ ξ : where the moment closure function is given by and g coefficients are defined as with the matrix C specified as Also, by construction one has that ψm( ν(t)) = vm (t) .
Finally, by combining (24) with (19) gives the XARNES system of equations, It is clear from the form of equation (19) that ρ 0 = 1 does not depend on time. For m = 0 the sum over contraction vectors can only contain one vector, c = 0 , and since α 0 = 0 the time derivative of ρ 0 vanishes. This expresses the fundamental probability conservation requirement for any reasonable theory. In contrast to the ρ 0 moment, a correlation moment with m = 0 is a real dynamic quantity. Based on this one might partition the Ω ξ space in two spaces. The first space should consist of the null vector, while the second space would consist of all other vectors in Ω ξ . It is possible to show that the values of g are same regardless whether the zero vector is singled out or not.

An intriguing similarity
A similar set of equations as the one given by (24-26) has been obtained in a slightly different context [10,11] where the multiplicative ansatz given in (24) was the starting point in developing a moment closure method for zero centered moments, η m , which in the notation of the present work would be given by The goal was to determine the values of the g coefficients in Eq (24) from the requirement that time derivatives of exact and approximate moments match at the time instance when the distribution function resembles the distribution function of the pure system. It was shown that the derivatives match best if the coefficients g are given exactly by (25). Thus the equations would be identical if not for the fact that the related equations in [10,11] are for zero-centered moments.
This remarkable coincidence where the same set of coefficients is obtained in two different ways is rather intriguing. In particular, this strongly suggests that the XARNES method might have advantageous properties as it comes to the derivative matching between the exact and the approximate moments computed for a suitable initial condition. It will be shown by employing a strict mathematical analysis, in the same way as done in [10,11], that this is indeed the case for the Poisson initial condition. This in turn explains the previous numerical observation from [14,15,17] that the XARNES method performs better when the particle number distribution function is Poisson-like as compared to the situation when the distribution resembles the one of the pure system.

Derivative matching setup
In here a similar procedure as in [10,11] will be carried out to show that the XARNES ansatz expressed in (24) and (25) has some advantageous derivative matching properties. The original procedure will be implemented as follows. Consider the Eqs. (19) and (27) at time t = t 0 where the particle number distribution function is strictly given by the uncorrelated multivariate Poisson distribution. In such a case all correlation forms are zero except the ones specified by vectors e i , i = 1, . . . , T . Thus the multivariate (uncorrelated) Poisson distribution is parameterized by parameters μ = (μ 1 , μ 2 , · · · , μ T ) ∈ R T where μ i = exp(w ei) . By a trivial application of Eq. (22) one can see that the exact factorial moments of the reacting system computed at t = t 0 are given by for any m ∈ . At t = t 0 the base factorial moments that are used in the XARNES method have to be chosen. This choice specifies the boundary condition for the XARNES equation of motion in (27). Naturally, for the purposes of comparing time derivatives it will be assumed that since if the values of the exact and the approximate base moments do not match their derivatives will not likely match either. The question is: given that the values of the base and the exact moments are same, can one expect that time derivatives of these quantities match as well?
In order to make some progress in answering this question a couple of identities involving g coefficients will be needed. These identities are hard to prove by using the explicit form of these coefficients given in (25). The formalism of generating functions will be used instead. Before doing that the following two definitions need to be stated, Definition 10. Symbol P ε ( ω) will be used to denote a polynomial where A m are arbitrary real valued coefficients. No other restrictions are imposed on the sum. Definition 11. Let P ε ( ω) be a polynomial as defined previously. Symbol ω m will be used to denote the operator that extracts the coefficient in front of a particular ω m term in the polynomial, i.e. A m : for any Ω * ⊂ Ω. Likewise, The following lemma is extremely useful for proving a couple of identities that will be needed later. The lemma has the form of a generating function identity. Lemma 1. Let coefficients g be defined by Eqs. (22) and (35). In other words, the only requirement imposed on these coefficients is that they are used to parameterize higher order factorial moments in terms of the base moments with the XARNES ansatz implied. In such a case these coefficients obey the following identity where ω and 1 ≡ (1,1,...,1) are vectors in R T andm ∈¯ ξ . Eq. (34) will be referred to as the generating function equation.
Proof. This identity can be proven as follows. First, one combines the XARNES ansatz (24) rewritten as By assuming that all correlation forms have predefined values given by where ω ∈ R T is arbitrary but fixed, leads to the equation that is central for the proof of the lemma, which by the use of the binomial theorem can be recognized as This turns the left hand side of (38) into The sum over the vectors in Ω ξ in the right hand side of (38) can be extended to include all vectors in m . Naturally, these terms have to be subtracted afterwards. By doing that one obtains By the use of the binomial theorem the first term on the right hand side of the equation can be recognized as the first term in the right hand of (34). Likewise, it is trivial to see that the second term in the equation can be characterized as P ξ +1 ( ω) . Thus the equation above becomes, Finally, the Lemma follows by using (38), (41), and (43). A couple of useful identities will be proven that follow from this Lemma and are stated as three corollaries. The first two corollaries can be proven easily without using the Lemma, e.g. as in [10,11]. In here they are proven in a different way to illustrate how to use the Lemma. The third corollary is a highly non trivial statement that would be very hard to prove by direct use of the explicit form of g coefficients. Corollary 1. Let coefficients g be given and let them satisfy the condition of Lemma 1. Then, Proof. The proof is trivial. One only needs to apply the operator ωm 1   Proof. First, one has to assume that where ω 1 , ω 2 ∈ R T are arbitrary but bound by the above constraint. Also, it is useful to realize that By using the above expression in the generating function formula, and (17), one obtains m∈ ξ γm m (1 + ω 1 )m(1 + ω 2 )m = In the third step one has to apply operators [ ω c 1 1 ] and [ ω c 2 1 ] to the generating function expression above. Applying the operators to the left hand side of the equation above gives the left hand side of Eq. (46). Likewise by applying the operators to the first term on the right hand side of the equation gives the right hand side of Eq. (46). Thus what is left to show is that the action of the operator onthe remaining second term results in zero.
The second term has the following structure where A coefficients can be found but their exact form is not relevant. Next, by using (17) and (13)  and by a simple change of variables p + s = c 1 and q + s = c 2 one arrives at To obtain the final form of the equation, and in particular the condition specified in the sum, one has to use the fact that c 1 + c 2 = p + q + 2 s ≥ p + q + s ≥ ξ + 1.
Finally, from the last equation one can see that, indeed, the application of the operator gives zero provided that c 1 + c 2 ≤ ξ . This proves the corollary. Please note that it would be very hard, if not impossible, to obtain such a result from the explicit expression for g coefficients given in (25).
Higher order generalizations of this corollary are possible. For example, by using one can easily prove that provided that vectors c 1 , c 2 , and c 3 satisfy c 1 + c 2 + c 3 ≤ ξ .
Now we are ready to prove some derivative matching results. The first question is whether the XARNES ansatz and the related moment closure function are consistent for the Poisson initial condition. This is answered in a form of the following Lemma.
Lemma 2. If the base and the exact factorial moments match at t = t 0 , and the particle number distribution function at this time instance is the Poisson distribution, i.e.
Proof. The direct use of the XARNES ansatz gives where (45) was used in the last step. Finally, we arrive at the central discussion of this work, i.e. the discussion of various time derivatives of exact and approximate moments and when and how they differ. For doing this it is instructive to investigate the difference between the exact and the XARNES equation systems. In this context one can easily show that the following equation system is valid, is any function, and [0] ≡ (t 0 ). The equation system above will be useful for proving a series of derivative matching theorems, in the similar vein as done in [10,11].

Three derivative matching theorems
Three derivative matching theorems will be proven, one theorem per derivative. The first two theorems have been proven in [10,11] for the pure state. In here they are proven for the Poisson state. The third theorem is entirely new.
The structure of the proofs is somewhat different than in [10,11] since in here the focus is on factorial moments. It seems that the equation system for factorial moments is more compact than for other types of moments. As an artifact of that, the theorems do not contain the error terms ε that were used in [10,11]. The theorems proven in here are more generic since they hold even for multi particle reactions, not just binary reactions. Again, as stated previously, all components of μ are taken strictly larger than zero. If one of the components is zero the XARNES ansatz does not work. Theorem 1. If the base and the exact factorial moments match at t = t 0 , and the particle number distribution function is the Poisson distribution, i.e., if then the first order derivatives in time also match for the base factorial moments: for allm ∈ ξ .
Proof. The theorem can be easily proven by considering Eq. (58) with h = 1. By assumptions of the theorem one has that p and if Ω R ⊂ Ω ξ , then the second order derivatives in time also match: Proof. The proof of this theorem is somewhat lengthier. To prove the theorem one needs to consider Eq. (58) with h = 2. If the assumptions of the theorem are valid, by theorem 1, p [1] m − ψ [1] m = 0 which eliminates the sum overm ∈ ξ in (58). Thus what is left to show is that a difference ρ [1] m − ψ [1] m withm ∈ r ξ vanishes. A straight forward application of the time derivative on the moment closure function leads to the following expression, Also, the use of the exact, and the XARNES equation systems to evaluate ρ [1] m and ψ [1] m gives ρ [1] This difference is zero provided for every c ∈ R . Please note that this condition is almost identical to the equation that characterize the g coefficients of the XARNES ansatz. In one replaces c withm , and Ω R with Ω ξ in the equation above, then the equation obtained in such a way would be identical to Eq. (25) or (44). Thus if Ω R ⊂ Ω ξ then the equation above is contained in the condition that defines the g coefficients, and the equation is automatically valid. This proves the theorem.
The third order derivatives will be investigated in the same vein as the first and the second order derivatives. The result will be formulated in a precise mathematical theorem. However, before stating the next theorem, it is useful to generalize the space of contraction vectors as follows.
Definition 12. Vector space of sums of contraction vectors c 1 + c 2 + · · · + c h where each of the vectors in the sum is from Ω R , will be denoted as and h is an integer and obeys h ≥ 1. and if Ω 2⊗R ⊂ Ω ξ , then the third order derivatives in time also match Proof. The theorem can be proven by considering Eq. (58) with h = 3. By theorem 2 one has that ρ [2] m − ψ [2] m = 0. What is left to show is that all differences ρ [2] m−ψ [2] m with m ∈¯ ξ vanish as well. Unfortunately this is a highly nontrivial task.
Somewhat naive recursive application of the exact equations of motion results in By using a tedious manipulation of the binomial coefficients the expression above can be converted into a more useful form In fact, it is easier to start from (77) and obtain (76). First, the manipulation requires that the sum of c 2 and d is changed, into the sum over c 2 = c 2 + d and c 3 = c 2 . After that the Vandermonde identity needs to be used which consumes the sum over c 3 , resulting finally in (77). By using the fact that Eq (77) can be written in the most useful form as The exact form of a coefficient α,β ( μ, c 1 , c 2 ) can be found if needed and, in fact, it is a series in μ . However, the exact form of these coefficients is not relevant for the discussion that follows. Let us use (79) in (74). This gives and m is zero if the following condition is met, for every c 1,2 ∈ 2⊗R . Eq. (81) is satisfied by the assumption of the theorem which states that Ω 2⊗R ⊂ Ω ξ . In such a case equations in (81) are a subset of the equations satisfied by the g coefficients which are given in (44) and are automatically valid. The three theorems proven so far are suggestive of the fact that one might try to prove the following conjecture.
Conjecture 1. If the base and the exact factorial moments match at t = t 0 , and the particle number distribution function is the Poisson distribution, i.e., if and if Ω h⊗R ⊂ Ω ξ , where h is an arbitrary integer such that h ≥ 1, then the time derivatives with orders D = 0, 1, 2, ..., h + 1 will also match Eventual. Inductive proof for general h could be used. However, the problems is that one would need to inspect a difference ρ m for arbitrary h and show that it vanishes under some conditions. Presumably, the main condition, apart from the standard requirements, e.g. such as having the Poisson initial condition, would be that Ω h⊗R ⊂ Ω ξ . For example, one can easily see that an expressions such as the one shown in (53) will appear if one tries to calculate higher time derivatives of ψm . However, as demonstrated for the h = 2 (D = 3) case the computation of ψ [h] m for higher values of h is a rather cumbersome and technical procedure. Generalization to higher orders is apparently very hard but not impossible.
There are couple of reasons why such a conjecture might be valid. First, the structure of the proofs of theorems 1-3 (the cases h = 0, 1, 2) suggests such a possibility. Second, there is numerical evidence from a previous study [17] that increase in ξ improves the accuracy of the XARNES method. In the context of the theorems discussed in here, increase in ξ implies that the initial Ω ξ set becomes larger. This in turncan make the condition Ω h⊗R ⊂ Ω ξ valid for a larger value of h. Finally as a result of that a larger number of derivatives would match which would explain the observed accuracy improvements in the studied benchmark cases. Finally, third, this conjecture has be checked by Mathematica for the T = 1 case and two binary reactions 2X 1 0 and 2X 1 X 1 , both with h = 0, 1, 2, 3, 4, 5 and ξ = 2, 4, 6, 8, 10 respectively, and one multi particle reaction 3X 1 2X 1 with h = 0, 1, 2, 3 and ξ = 3, 6, 9.

Conclusions
The three theorems explain the mechanism behind the numerically observed fact that the XARNES method works well if the particle number distribution function is close to the Poisson distribution. Thus if all correlation forms are in some sense small, the XARNES method will provide a very accurate result.
For very small molecule counts one should expect problems with moment closure formulas like (24) with negative g coefficients. The XARNES ansatz can describe such situation but one needs to consider a limiting processwhere factorial moments approach zero from the above. For example, one might want to start the systems from a state with the Poisson distribution, which is natural in many biologically relevant cases, but with some components of μ equal to zero. This cannot be done directly since the XARNES ansatz breaks down. In more practical terms, any decent numerical Ordinary Differential Equation (ODE) solver should issue a warning for such an initial state. To start the system from a state wheresome copy numbers are zero it is necessary to consider increasingly smaller values for such copy numbers.
Previous numerical studies showed that there are systems for which the XARNES ODE system develops a singularity and the numerical solver has to stop [14,15,17]. Unfortunately, the theorems that have been proven cannot say anything about such singularities since the particle number distribution function becomes increasingly different from the Poisson distribution.