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
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 ni, 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 r1, r2,..., rR 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 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 mi 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.
Definition 2. It is useful to introduce the order, or norm, of a vector 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 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 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.
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 and be two arbitrary tuples of real numbers with rank T, such as . The set of all such vectors will be denoted by RT. For any pair of such vectors the generalized power will be defined as
Obviously, Ω⊂RT 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 from Ω is formally defined as
and it is assumed that a binomial coefficient in the product is zero if m2i > m1i.
Definition 8. The factorial-like symbol applied to a vector from Ω is generalized as
Definition 9. The product between two real numbers is generalized as
where and are two arbitrary vectors from RT. Please note that the result of the operation is an element in the same set, . 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 and is defined as
It was shown in  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 . 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.
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 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.
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 with in terms of the base moments with :
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,
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
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, , which in the notation of the present work would be given by
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
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 = t0 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 . Thus the multivariate (uncorrelated) Poisson distribution is parameterized by parameters where . By a trivial application of Eq. (22) one can see that the exact factorial moments of the reacting system computed at t = t0 are given by
for any .
At t = t0 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 will be used to denote a polynomial
where are arbitrary real valued coefficients. No other restrictions are imposed on the sum.
Definition 11. Let be a polynomial as defined previously. Symbol will be used to denote the operator that extracts the coefficient in front of a particular term in the polynomial, i.e. :
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 and 1 ≡ (1,1,...,1) are vectors in RT and . 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
with the definition of the correlation forms (22) to obtain
By assuming that all correlation forms have predefined values given by
where is arbitrary but fixed, leads to the equation that is central for the proof of the lemma,
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 . 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 . 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 γ 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 on the both sides of the generating function equation (34).
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 at the point .
Corollary 3. Let the coefficients γ satisfy the generating function equation (34). Then following identity holds
provided vectors and satisfy .
Proof. First, one has to assume that
where 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
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.
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 and 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 . 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).
Higher order generalizations of this corollary are possible. For example, by using
one can easily prove that
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.
Lemma 2. If the base and the exact factorial moments match at t = t0, and the particle number distribution function at this time instance is the Poisson distribution, i.e.
then non-base 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 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.
Theorem 1. If the base and the exact factorial moments match at t = t0, 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 which eliminates the sum over in (58). From Lemma 2 it follows that which eliminates the sum over . This finally proves the theorem.
Theorem 2. If the base and the exact factorial moments match at t = t0, 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, which eliminates the sum over in (58). Thus what is left to show is that a difference with 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 and gives
This difference is zero provided
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.
Definition 12. Vector space of sums of contraction vectors where each of the vectors in the sum is from ΩR, will be denoted as
and h is an integer and obeys h ≥ 1.
Theorem 3. If the base and the exact factorial moments match at t = t0 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 . What is left to show is that all differences with vanish as well. Unfortunately this is a highly nontrivial task.
By the use of the standard calculus one can show that
By using the XARNES equations, and identity (44) for γ coefficients, the can be expressed as
which vanishes provided
for any . This is indeed true by corollary 3 under the assumptions of the theorem.
What is left to show is that . A strategy for proving this is as follows. If one could show that the following identity holds for the exact moments
then it is clear that would be true since one could write
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 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).
By using the fact that
Eq (77) can be written in the most useful form as
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.
Let us use (79) in (74). This gives
and is zero if the following condition is met,
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.
Conjecture 1. If the base and the exact factorial moments match at t = t0, 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 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.
The authors declare that they have no competing interests.
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
Proceedings of the 45th Ieee Conference on Decision and Control, Vols 1-14 2006, 2063-2068.
IEEE Conference on Decision and Control
J Theor Biol 2011.