Abstract
The Xlevel 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 Poissonlike. 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 [515].
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 Poissonlike 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 [911] 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 welldefined 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
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_{α }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
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
Definition 2. It is useful to introduce the order, or norm, of a vector
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
where ξ ≥ 1 is an arbitrary integer used to classify various theories.
Definition 4. In a similar way,
will denote the set of vectors
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
Definition 5. Finally, the set of all vectors that are not in Ω_{ξ }will be denoted as,
A vector from this set will be denoted by the bar character above the vector symbol,
e.g.
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
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
and it is assumed that a binomial coefficient
Definition 8. The factoriallike symbol applied to a vector
Definition 9. The product between two real numbers is generalized as
where
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
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
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
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
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
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 nonbase factorial moments
where the moment closure function is given by
and γ coefficients are defined as
with the matrix C specified as
Also, by construction one has that
Finally, by combining (24) with (19) gives the XARNES system of equations,
for
It is clear from the form of equation (19) that
An intriguing similarity
A similar set of equations as the one given by (2426) 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,
The goal was to determine the values of the γ 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 γ 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 zerocentered 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 Poissonlike 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
for any
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 γ 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
where
Definition 11. Let
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 γ 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
Proof. This identity can be proven as follows. First, one combines the XARNES ansatz (24) rewritten as
with the definition of the correlation forms (22) to obtain
By assuming that all correlation forms have predefined values given by
where
First, we focus on the left hand side of the above equation. From (37), and (22) with
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
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
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 γ coefficients.
Corollary 1. Let coefficients γ be given and let them satisfy the condition of Lemma 1. Then,
Proof. The proof is trivial. One only needs to apply the operator
Corollary 2. Let coefficients γ be given and let them satisfy the condition of Lemma 1. Then,
Proof. The identity can be obtained by evaluating the gradient of Eq. (34) with respect
to
Corollary 3. Let the coefficients γ satisfy the generating function equation (34). Then following identity holds
provided vectors
Proof. First, one has to assume that
where
By using the above expression in the generating function formula, and (17), one obtains
In the third step one has to apply operators
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) the equation becomes
and by a simple change of variables
To obtain the final form of the equation, and in particular the condition specified
in the sum, one has to use the fact that
Higher order generalizations of this corollary are possible. For example, by using
one can easily prove that
provided that vectors
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.
then nonbase moments also match:
for all
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,
and notation
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
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 all
Proof. The theorem can be easily proven by considering Eq. (58) with h = 1. By assumptions of the theorem one has that
Theorem 2. 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 Ω_{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,
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
This difference is zero provided
for every
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
and h is an integer and obeys h ≥ 1.
Theorem 3. 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 Ω_{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
By the use of the standard calculus one can show that
where
and
By using the XARNES equations, and identity (44) for γ coefficients, the
which vanishes provided
for any
What is left to show is that
then it is clear that
and this expression would vanish by Theorem 2.
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
By using the fact that
Eq (77) can be written in the most useful form as
The exact form of a coefficient
Let us use (79) in (74). This gives
and
for every
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
There are couple of reasons why such a conjecture might be valid. First, the structure of the proofs of theorems 13 (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 γ 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
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.
Competing interests
The authors declare that they have no competing interests.
References

Eldar A, Elowitz MB: Functional roles for noise in genetic circuits.

Lestas I, Vinnicombe G, Paulsson J: Fundamental limits on the suppression of molecular fluctuations.

Gillespie DT: Exact stochastic simulation of coupled chemicalreactions.

Elf J, Ehrenberg M: Fast evaluation of fluctuations in biochemical networks with the linear noise approximation.

Sasai M, Wolynes PG: Stochastic gene expression as a manybody problem.

Engblom S: Computing the moments of high dimensional solutions of the master equation.

Singh A, Hespanha JP: Approximate Stochastic Models for Chemically Reacting Systems. [http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.135.8468] webcite
2006.

Singh A, Hespanha JP: Lognormal moment closures for biochemical reactions.
Proceedings of the 45th Ieee Conference on Decision and Control, Vols 114 2006, 20632068.
IEEE Conference on Decision and Control

Singh A, Hespanha JP: A derivative matching approach to moment closure for the stochastic logistic model.

Lee CH, Kim KH, Kim P: A moment closure method for stochastic reaction networks.

Grima R: Investigating the robustness of the classical enzyme kinetic equations in small intracellular compartments.

Konkoli Z: Spontaneous noise reduction in a strongly cooperative reaction model.

Walczak AM, Sasai M, Wolynes PG: Selfconsistent proteomic field theory of stochastic gene switches.

Konkoli Z: Modelling reaction noise with a desired accuracy by using the X level Approach Reaction Noise Estimator (XARNES) method.
J Theor Biol 2011.
submitted

Kotomin E, Kuzovkov V: Modern aspects of diffusioncontrolled reactions: cooperative phenomena in bimolecular processes, Volume 34 of Comprehensive chemical kinetics. Amsterdam: Elsevier; 1996.

Konkoli Z: Application of Bogolyubov's theory of weakly nonideal Bose gases to the A + A, A + B, B + B reactiondiffusion system.