Mathematical explanation of the predictive power of the X-level approach reaction noise estimator method
© Konkoli; licensee BioMed Central Ltd. 2012
Received: 11 January 2012
Accepted: 13 April 2012
Published: 13 April 2012
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.
Noise is an integral part of the workings of the living cell biochemistry . There are many types of noise and this work focuses on the intrinsic noise. If reactant copy numbers are low they can fluctuate widely . These fluctuations can severely influence the dynamics of the cell and need to be carefully controlled .
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 . 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–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 well-defined 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  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  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  willbe closely followed. This report was formalized in much shorter form in . 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 . 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
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.
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 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.
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.
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 where the digit 1 appears on the i-the place.
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.
where ξ ≥ 1 is an arbitrary integer used to classify various theories.
will denote the set of vectors that are smaller than a vector 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 with . Sets of the type will be used to limit various sums.
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.
and it is assumed that a binomial coefficient in the product is zero if m2i> m1i.
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
for all . A vector 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 or the sum over the contraction vectors in (19) is restricted.
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 there are terms that involve higher order moments since it may happen that . 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
where denotes the correlation form labeled by a vector . The detailed motivation behind this equation is given in . It is assumed that that correlation forms with the orders above a given threshold ξ can be assumed small, . Since a moment is not exactly equal to the related moment , two different symbols, v and ρ, 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.
Also, by construction one has that .
It is clear from the form of equation (19) that does not depend on time. For the sum over contraction vectors can only contain one vector, , and since the time derivative of vanishes. This expresses the fundamental probability conservation requirement for any reasonable theory. In contrast to the moment, a correlation moment with 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 γ are same regardless whether the zero vector is singled out or not.
An intriguing similarity
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 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
for any .
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,
where are arbitrary real valued coefficients. No other restrictions are imposed on the sum.
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.
where and 1 ≡ (1,1,...,1) are vectors in R T and . Eq. (34) will be referred to as the generating function equation.
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.
Proof. The proof is trivial. One only needs to apply the operator on the both sides of the generating function equation (34).
Proof. The identity can be obtained by evaluating the gradient of Eq. (34) with respect to at the point .
provided vectors and satisfy .
In the third step one has to apply operators and 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.
To obtain the final form of the equation, and in particular the condition specified in the sum, one has to use the fact that . Finally, from the last equation one can see that, indeed, the application of the operator gives zero provided that . 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 γ coefficients given in (25).
provided that vectors , , and satisfy .
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.
for all .
where (45) was used in the last step.
and notation with h = 1,2,3,... is implied where φ(t) is any function, and φ ≡ φ(t0). 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.
for all .
Proof. The theorem can be easily proven by considering Eq. (58) with h = 1. By assumptions of the theorem one has that which eliminates the sum over in (58). From Lemma 2 it follows that which eliminates the sum over . This finally proves the theorem.
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, which eliminates the sum over in (58). Thus what is left to show is that a difference with vanishes.
for every . Please note that this condition is almost identical to the equation that characterize the γ coefficients of the XARNES ansatz. In one replaces with , 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 γ 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.
and h is an integer and obeys h ≥ 1.
Proof. The theorem can be proven by considering Eq. (58) with h = 3. By theorem 2 one has that . What is left to show is that all differences with vanish as well. Unfortunately this is a highly nontrivial task.
for any . This is indeed true by corollary 3 under the assumptions of the theorem.
and this expression would vanish by Theorem 2.
In fact, it is easier to start from (77) and obtain (76). First, the manipulation requires that the sum of and is changed, into the sum over and After that the Vandermonde identity needs to be used which consumes the sum over , resulting finally in (77).
The exact form of a coefficient 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.
for every . 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 γ 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.
Eventual. Inductive proof for general h could be used. However, the problems is that one would need to inspect a difference 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 . However, as demonstrated for the h = 2 (D = 3) case the computation of 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  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 2X1 → 0 and 2X1 → X1, both with h = 0, 1, 2, 3, 4, 5 and ξ = 2, 4, 6, 8, 10 respectively, and one multi particle reaction 3X1 → 2X1 with h = 0, 1, 2, 3 and ξ = 3, 6, 9.
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 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.
- Eldar A, Elowitz MB: Functional roles for noise in genetic circuits. Nature. 2010, 467: 167-173.PubMed CentralView ArticlePubMedGoogle Scholar
- Paulsson J: Summing up the noise in gene networks. Nature. 2004, 427: 415-418.View ArticlePubMedGoogle Scholar
- Lestas I, Vinnicombe G, Paulsson J: Fundamental limits on the suppression of molecular fluctuations. Nature. 2010, 467: 174-178.PubMed CentralView ArticlePubMedGoogle Scholar
- Gillespie DT: Exact stochastic simulation of coupled chemical-reactions. J Phys Chem. 1977, 81: 2340-2361.View ArticleGoogle Scholar
- Elf J, Ehrenberg M: Fast evaluation of fluctuations in biochemical networks with the linear noise approximation. Genome Res. 2003, 13: 2475-2484.PubMed CentralView ArticlePubMedGoogle Scholar
- Nasell I: An extension of the moment closure method. Theor Popul Biol. 2003, 64: 233-239.View ArticlePubMedGoogle Scholar
- Sasai M, Wolynes PG: Stochastic gene expression as a many-body problem. Proc Natl Acad Sci USA. 2003, 100: 2374-2379.PubMed CentralView ArticlePubMedGoogle Scholar
- Engblom S: Computing the moments of high dimensional solutions of the master equation. Appl Math Comput. 2006, 180: 498-515.View ArticleGoogle Scholar
- Singh A, Hespanha JP: Approximate Stochastic Models for Chemically Reacting Systems. 2006, http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.135.8468http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.135.8468Google Scholar
- Singh A, Hespanha JP: Lognormal moment closures for biochemical reactions. Proceedings of the 45th Ieee Conference on Decision and Control, Vols 1-14. 2006, 2063-2068. IEEE Conference on Decision and ControlView ArticleGoogle Scholar
- Singh A, Hespanha JP: A derivative matching approach to moment closure for the stochastic logistic model. Bull Math Biol. 2007, 69: 1909-1925.View ArticlePubMedGoogle Scholar
- Lee CH, Kim KH, Kim P: A moment closure method for stochastic reaction networks. J Chem Phys. 2009, 130: 134107-View ArticlePubMedGoogle Scholar
- Grima R: Investigating the robustness of the classical enzyme kinetic equations in small intracellular compartments. BMC Syst Biol. 2009, 3: 101-PubMed CentralView ArticlePubMedGoogle Scholar
- Konkoli Z: Multiparticle reaction noise characteristics. J Theor Biol. 2011, 271: 78-86.View ArticlePubMedGoogle Scholar
- Konkoli Z: Spontaneous noise reduction in a strongly cooperative reaction model. J Theor Biol. 2011, 285: 96-102.View ArticlePubMedGoogle Scholar
- Walczak AM, Sasai M, Wolynes PG: Self-consistent proteomic field theory of stochastic gene switches. Biophys J. 2005, 88: 828-850.PubMed CentralView ArticlePubMedGoogle Scholar
- Konkoli Z: Modelling reaction noise with a desired accuracy by using the X level Approach Reaction Noise Estimator (XARNES) method. J Theor Biol. 2011, submittedGoogle Scholar
- Kotomin E, Kuzovkov V: Modern aspects of diffusion-controlled reactions: cooperative phenomena in bimolecular processes, Volume 34 of Comprehensive chemical kinetics. 1996, Amsterdam: ElsevierGoogle Scholar
- Konkoli Z: Application of Bogolyubov's theory of weakly nonideal Bose gases to the A + A, A + B, B + B reaction-diffusion system. Phys Rev E. 2004, 69: 011106-View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.