A Symbolic Computation Approach to a Problem Involving Multivariate Poisson Distributions^{∗}

Eduardo D. Sontag and Doron Zeilberger Department of Mathematics Rutgers University Hill Center-Busch Campus 110 Frelinghuysen Rd., Piscataway, NJ 08854-8019, USA.

              [sontag,zeilberg] at math dot rutgers dot edu

June 13, 2009

Abstract

Multivariate Poisson random variables subject to linear integer constraints arise in several application areas, such as queuing and biomolecular networks. This note shows how to compute conditional statistics in this context, by employing Wilf-Zeilberger Theory and associated algorithms. A symbolic computation package has been developed and is made freely available. A discussion of motivating biomolecular problems is also provided.

1 Introduction

In application areas such as queuing and biomolecular networks, one is often interested in the study of independent Poisson random variables subject to side information represented by linear integer constraints. We show how to reduce the computation of conditional statistics for this problem to the evaluation of coefficients of generating functions. These coefficients can, in turn, be computed using Wilf-Zeilberger (WZ) theory. We discuss this reduction, and make available a symbolic computation package developed for that purpose.

We next provide a formulation of the problem, and briefly indicate its motivations. In Section 2, we explain the reduction to exponential type generating functions, and in Section 3 we discuss the fact that recurrences can be obtained for their coefficients. Section 4 discusses the special case of just two side constraints, which is considerably simpler. Section 5 illustrates the use of the symbolic package through a number of examples, all of which arise from the biomolecular networks discussed in Section 6. An Appendix includes a proof of the basic representation theorem which enables application of this techniques to certain reaction networks.

Suppose that we have n independent Poisson random variables, X_{j} (j = 1 … n), with parameters λ_{j} respectively. In other words

Pr (X_{1} = k_{1}, X_{2} = k_{2}, …, X_{n} = k_{n}) = e^{−(λ_{1}+...+λ_{n})} \frac{λ_{1}^{k_{1}}}{k_{1}!} \frac{λ_{2}^{k_{2}}}{k_{2}!} … \frac{λ_{n}^{k_{n}}}{k_{n}!} .  (1)

∗Accompanied by Maple package MVPoisson downloadable from http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/mvp.html

Suppose that we can’t observe the X_{j}’s directly, but only a certain number, m, of linear combinations of them:

Y_{i} = ∑_{j=1}^{n} a_{ij}X_{j} , (i = 1, …, m) ,

where A = (a_{ij}) is a certain m × n matrix with non-negative coefficients.

We are interested in the following questions:

1. Can one compute (fast!), for any given vector (b_{1}, …, b_{m}) (possibly with large coordinates), the probability

F (b_{1}, …, b_{m}) := Pr (Y_{1} = b_{1}, …, Y_{m} = b_{m}) .

2. Can one compute (fast!), for any given vector (b_{1}, …, b_{m}), (possibly with large coordinates) the conditional expectation

G_{j}(b_{1}, …, b_{m}) := E[X_{j} | Y_{1} = b_{1}, …, Y_{m} = b_{m}] , (1 ≤ j ≤ n).

3. More generally, can one compute (fast!), the higher moments

G_{j}^{(r)}(b_{1}, …, b_{m}) := E[X_{j}^{r} | Y_{1} = b_{1}, …, Y_{m} = b_{m}] , (r ≥ 2) ,

that would immediately allow us to compute the moments about the mean. Can we compute (fast!) mixed moments, in particular the covariances?

For example, suppose that X_{i} is Poisson with parameter λ_{i}, i = 1, 2, X_{1} and X_{2} are independent, and A = (1 1). Thus, Y = X_{1} + X_{2} is Poisson with parameter λ_{1} + λ_{2}. Fix a non-negative integer b. The probability that X_{1} = k given that Y = X_{1} + X_{2} = b is:

e^{−(λ_{1}+λ_{2})} \frac{λ_{1}^{k}}{k!} \frac{λ_{2}^{b−k}}{(b − k)!} / e^{−(λ_{1}+λ_{2})} \frac{(λ_{1} + λ_{2})^{b}}{b!}

which equals

\binom{b}{k} p^{k}(1 − p)^{b−k}

with p = \frac{λ_{1}}{λ_{1}+λ_{2}} . It follows that (X_{1}|Y = b) is a binomial random variable B(b, p), and similarly (X_{2}|Y = b) is a binomial random variable B(b, 1 − p). Statistics for binomial variables (means, variances, and all moments) are well-known and easy to compute. On the other hand, for more complicated linear constraints, and especially if more than one such constraint is imposed, statistics become considerably harder to obtain.

A simple example of where this type of problem might arise is as follows. Suppose that the random variables X_{i} count the number of calls placed, during a typical time period, to an international service center and originating from a specific country or geographical area and in a specific customer language. For example, X_{1} may represent the number of English-speaking callers from the USA, X_{2} the number of Spanish-speaking callers from the USA, X_{3} the number of English-speaking callers from Latin America, X_{4} the number of Spanish-speaking callers from Latin America, X_{5} the number of English-speaking callers from the UK, and X_{6} the number of Spanish-speaking callers from the UK. It is natural to assume that each of the random variables is Poisson-distributed. Now, suppose that we want to know what are the statistics of the variable X_{1}, for example, the variance in the number of English-speaking callers from the USA, subject to the additional information that the total number of Spanish-language calls received was 100 and that the number of calls received from the US was 50. That is, we are interested in the statistics of X_{1} conditioned on Y_{1} = 100, Y_{2} = 50 with Y_{1} = X_{2} + X_{4} + X_{6} and Y_{2} = X_{1} + X_{2}. (More interestingly, one might have mixed information, represented by more general linear combinations.) We were originally motivated in this work by applications in molecular biology; we defer to Section 6 a detailed discussion and examples.

2 The generating function

Fix a matrix A = (a_{ij}) (1 ≤ i ≤ m, 1 ≤ j ≤ n), once and for all. Let

F_{0}(b_{1}, …, b_{m}) = ∑_{{k_{1},...,kn≥0 \atop a_{11}k_{1}+...+a_{1n}kn=b_{1},...,a_{m1}k_{1}+...+amnkn=bm}} \frac{λ_{1}^{k_{1}}}{k_{1}!} \frac{λ_{2}^{k_{2}}}{k_{2}!} … \frac{λ_{n}^{k_{n}}}{k_{n}!}

(value is zero if the sum is empty). Thus, our focus will be on computing F_{0}, from which we can easily obtain F, since

F (b_{1}, …, b_{m}) = e^{−(λ_{1}+...+λ_{n})}F_{0}(b_{1}, …, b_{m}) .

Let f_{0} be the (multivariable) generating function of F_{0}, in other words

f_{0}(z_{1}, …, z_{m}) = ∑_{b_{1}≥0,...,b_{m}≥0} F_{0}(b_{1}, …, b_{m})z_{1}^{b_{1}} … z_{m}^{b_{m}} .

Our quantity of interest, F_{0}(b_{1}, …, b_{m}), is the coefficient of z_{1}^{b_{1}} … z_{m}^{b_{m}} in the multivariable Taylor expansion about the origin of f_{0}(z_{1}, …, z_{m}).

We have:

f_{0}(z_{1}, …, z_{m}) = ∑_{b_{1}≥0,...,b_{m}≥0} ( ∑_{{k_{1},...,kn≥0 \atop a_{11}k_{1}+...+a_{1n}kn=b_{1},...,a_{m1}k_{1}+...+amnkn=bm}} \frac{λ_{1}^{k_{1}}}{k_{1}!} \frac{λ_{2}^{k_{2}}}{k_{2}!} … \frac{λ_{n}^{k_{n}}}{k_{n}!} ) z_{1}^{b_{1}} … z_{m}^{b_{m}} .

By changing the order of summation, this equals

∑_{k_{1}≥0,...,k_{n}≥0} \frac{λ_{1}^{k_{1}}}{k_{1}!} \frac{λ_{2}^{k_{2}}}{k_{2}!} … \frac{λ_{n}^{k_{n}}}{k_{n}!} z_{1}^{a_{11}k_{1}+...+a_{1n}k_{n}} … z_{m}^{a_{m1}k_{1}+...+a_{mn}k_{n}}

= ∑_{k_{1}≥0,...,k_{n}≥0} \frac{(λ_{1}z_{1}^{a_{11}}z_{2}^{a_{21}} … z_{m}^{a_{m1}})^{k_{1}}}{k_{1}!} … \frac{(λ_{n}z_{1}^{a_{1n}}z_{2}^{a_{2n}} … z_{m}^{a_{mn}})^{k_{n}}}{k_{n}!}

= ( ∑_{k_{1}≥0} \frac{(λ_{1}z_{1}^{a_{11}}z_{2}^{a_{21}} … z_{m}^{a_{m1}})^{k_{1}}}{k_{1}!} ) … ( ∑_{k_{n}≥0} \frac{(λ_{n}z_{1}^{a_{1n}}z_{2}^{a_{2n}} … z_{m}^{a_{mn}})^{k_{n}}}{k_{n}!} )
= exp(λ_{1}z_{1}^{a_{11}}z_{2}^{a_{21}} … z_{m}^{a_{m1}}) … exp(λ_{n}z_{1}^{a_{1n}}z_{2}^{a_{2n}} … z_{m}^{a_{mn}})
= exp (λ_{1}z_{1}^{a_{11}}z_{2}^{a_{21}} … z_{m}^{a_{m1}} + … + λ_{n}z_{1}^{a_{1n}}z_{2}^{a_{2n}} … z_{m}^{a_{mn}}) .

We have just derived Theorem 1:

f_{0}(z_{1}, …, z_{m}) = exp (∑_{j=1}^{n} λ_{j} ∏_{i=1}^{m} z_{i}^{a_{ij}} ) .

The conditional probability

Pr (X_{1} = k_{1}, X_{2} = k_{2}, …, X_{n} = k_{n} | Y_{1} = b_{1}, …, Y_{m} = b_{m})

is the same as the expression in (1) divided by F (b), provided that ∑_{j=1}^{n} a_{ij}k_{j} = b_{i} for all i, and is zero otherwise. Recall that the rth factorial moment of a random variable W, E[W^{(r)}], is, by definition, the expectation of W!/(W − r)!. We are interested in the conditional factorial moments of X_{j} given Y = b, which we will denote as E[X_{j}^{(r)} | Y]. By definition, E[X_{j}^{(r)} | Y] is the following expression divided by F_{0}(b):

∑_{{k_{1},...,kn≥0 \atop a_{11}k_{1}+...+a_{1n}kn=b_{1},...,a_{m1}k_{1}+...+amnkn=bm}} k_{j}(k_{j} − 1) … (k_{j} − r + 1)\frac{λ_{1}^{k_{1}}}{k_{1}!} \frac{λ_{2}^{k_{2}}}{k_{2}!} … \frac{λ_{n}^{k_{n}}}{k_{n}!} .  (2)

Now, expression (2) is the same as the result of applying the operator λ_{j}^{r}(\frac{∂}{∂λ_{j}})^{r} to F_{0}(b_{1}, …, b_{m}) when viewing the λ’s as variables and not as constants. On the other hand,

λ_{j}^{r}(\frac{∂}{∂λ_{j}})^{r}f_{0}(z_{1}, …, z_{m}) = ∑_{b_{1}≥0,...,b_{m}≥0} λ_{j}^{r}(\frac{∂}{∂λ_{j}})^{r}F_{0}(b_{1}, …, b_{m})z_{1}^{b_{1}} … z_{m}^{b_{m}}

and therefore expression (2) is the same as the coefficient of z_{1}^{b_{1}} … z_{m}^{b_{m}} in λ_{j}^{r}(\frac{∂}{∂λ_{j}})^{r}f_{0}(z_{1}, …, z_{m}). Since, as formal power series, we have the representation in Theorem 1, we conclude that expression (2) is the same as the coefficient of z_{1}^{b_{1}} … z_{m}^{b_{m}} in (∏_{i=1}^{m} z_{i}^{a_{ij}})^{r}f_{0}(z), which is the same as F (b_{1} −ra_{1j}, b_{2} −ra_{2j}, …, b_{m} −ra_{mj}) when all b_{i} −ra_{ij} ≥ 0 and zero otherwise. In conclusion, E[X_{j}^{(r)} | Y] equals λ_{j}^{r} · F_{0}(b_{1} − ra_{1j}, b_{2} − ra_{2j}, …, b_{m} − ra_{mj}) divided by F_{0}(b). We have proved:

Theorem 2: The conditional factorial moments E[X_{j}^{(r)} | Y] are given in terms of the F_{0}(b_{1}, …, b_{m}) by

λ_{j}^{r} · \frac{F_{0}(b_{1} − ra_{1j}, b_{2} − ra_{2j}, …, b_{m} − ra_{mj})}{F_{0}(b_{1}, …, b_{m})}

when all b_{i} − ra_{ij} ≥ 0 and zero otherwise.

So everything depends on a fast computation of the coefficients F_{0}(b_{1}, …, b_{m}), of f_{0}(z_{1}, …, z_{m}).

By taking mixed partial derivatives, we can easily derive analogous expressions for mixed moments, in particular, the covariances.

3 Recurrences

From now on, let’s assume that the entries of A, (a_{ij}), are non-negative integers. In that case, we can write

f_{0}(z) = exp(Q(z_{1}, …, z_{m})) ,

where Q(z_{1}, …, z_{m}) is the polynomial

Q(z_{1}, …, z_{m}) := ∑_{j=1}^{n} λ_{j} ∏_{i=1}^{m} z_{i}^{a_{ij}} .

By Cauchy’s theorem, we can express F (b_{1}, …, b_{m}) as a multi-contour integral:

F (b_{1}, …, b_{m}) = (\frac{1}{2πi})^{m} ∫_{|z_{1}|=c} … ∫_{|z_{m}|=c} \frac{exp(Q(z_{1}, …, z_{m}))}{z_{1}^{b_{1}+1} … z_{m}^{b_{m}+1}} dz_{1} … dz_{m} .

By the celebrated Wilf-Zeilberger theory ([16]), F (b_{1}, …, b_{m}) satisfies pure linear recurrences with polynomial coefficients in each of its arguments. This means that for each i between 1 and m, there exists a positive integer R_{i} (the order) and polynomials P_{r}^{(i)}(b_{1}, …, b_{m}) (0 ≤ r ≤ R_{i}) such that the following holds, for all (b_{1}, …, b_{m}):

∑_{r=0}^{R_{i}} P_{r}^{(i)}(b_{1}, …, b_{m})F (b_{1}, …, b_{i−1}, b_{i} + r, b_{i+1}, …, b_{m}) = 0 .

Once these recurrences are known, one can compute F (b_{1}, …, b_{m}) in time linear in b_{1} +…+b_{m} and with constant memory allocation (one only needs to remember, at each stage, a constant number of values).

In rare cases, the leading term of the recurrence would vanish, in which case, we would encounter a (discrete) “singularity”, and would not be able to go on, since we would have to divide by 0, but in that case one can show that there is an alternative route, using another order of applying the recurrences.

The Apagodu-Zeilberger[3] multi-variable extension of the Almkvist-Zeilberger[1] algorithm can find such recurrences explicitly. Unfortunately, for matrices A with more than three rows, the time taken to find such recurrences is prohibitive, but many matrices of interest have two or three rows.

4 Two-Rowed matrices

If the matrix A only has two rows, and its entries are only {0, 1}, then one can express F (b_{1}, b_{2}) as a single sum. Indeed, let

• c_{10} be the sum of the λ_{j}’s for which a_{1,j} = 0, a_{2,j} = 1,

• c_{01} be the sum of the λ_{j}’s for which a_{1,j} = 1, a_{2,j} = 0,

• c_{11} be the sum of the λ_{j}’s for which a_{1,j} = 1, a_{2,j} = 1.

Then, we have

Q(z) = c_{01}z_{1} + c_{10}z_{2} + c_{11}z_{1}z_{2},

and so

f_{0}(z_{1}, z_{2}) = e^{Q(z)} = ∑_{k=0}^{∞} \frac{Q(z)^{k}}{k!}

= ∑_{α≥0,β≥0,γ≥0} \frac{(c_{01}z_{1})^{α}(c_{10}z_{2})^{β}(c_{11}z_{1}z_{2})^{γ}}{α!β!γ!}

= ∑_{α≥0,β≥0,γ≥0} \frac{c_{01}^{α}c_{10}^{β}c_{11}^{γ}z_{1}^{α+γ}z_{2}^{β+γ}}{α!β!γ!} .

To get F_{0}(b_{1}, b_{2}), we must extract the coefficient of z_{1}^{b_{1}}z_{2}^{b_{2}} which entails α = b_{1} − γ, β = b_{2} − γ, and we have the single-sum binomial coefficient (hypergeometric) sum (replacing γ by k)

F_{0}(b_{1}, b_{2}) = ∑_{k=0}^{min(b_{1},b_{2})} \frac{c_{11}^{k}c_{01}^{b_{1}−k}c_{10}^{b_{2}−k}}{k!(b_{1} − k)!(b_{2} − k)!} .

Using the Zeilberger Algorithm ([17, 13]) , we get the following linear recurrence:

(c_{10}b_{1}^{2} + 4c_{10} − 2c_{10}b_{2} + 4c_{10}b_{1} − c_{10}b_{1}b_{2})F_{0}(b_{1} + 2, b_{2})+

(−c_{11}b_{1}−c_{11}+b_{2}c_{11}+b_{2}c_{10}c_{01}−2b_{1}c_{10}c_{01}−3c_{01}c_{10})F_{0}(b_{1}+1, b_{2})+(c_{11}c_{10}+c_{01}c_{10}^{2})F_{0}(b_{1}, b_{2}) = 0

and an analogous formula holds for a recursion on b_{2}.

5 The Maple package MVPoisson

All this is implemented in the Maple package MVPoisson accompanying this article. It is available from the webpage of this article

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/mvp.html ,

where one can also find sample input and output.

We next discuss several examples of matrices A and computations using MVPoisson. These examples, of interest in themselves, are motivated by the biochemical networks discussed in Section 6.

Mostly, we illustrate the use of the command “RecsV”, which provides the recurrences satisfied by the coefficients F_{0}, but we also show a few examples of other commands that compute moments.

5.1 A one-row example

The matrix A is:

A = (1 1) .  (3)

As discussed in the Introduction, the conditional random variables (X_{i}|Y = b) are binomial. With the notations of this paper,

F_{0}(b) = ∑_{i+j=b} \frac{λ_{1}^{i}}{i!} \frac{λ_{2}^{j}}{j!} = \frac{1}{b!}(λ_{1} + λ_{2})^{b} .

This function F_{0} clearly satisfies the following recurrence:

F_{0}(b + 1) = \frac{λ_{1} + λ_{2}}{b + 1} F_{0}(b)

with F_{0}(0) = 0 and F (1) = λ_{1} + λ_{2}. Indeed, for the matrix A in (3), the “RecsV(A, λ, b)” command provides the following recurrence:

F_{0}(b_{1} + 1) = \frac{λ_{1} + λ_{2}}{1 + b_{1}} F_{0}(b_{1})

with initial condition F_{0}(1) = λ_{1} + λ_{2}.

5.2 A two-row example

The matrix A is:

A = ( \begin{matrix}1 & 0 & 1 \\ 0 & 1 & 1 \end{matrix} ) .  (4)

For the matrix A in (4), the “RecsV(A, λ, b)” command provides the following two-dimensional recurrence:

F_{0}(b_{1} + 2, b_{2}) = − \frac{−b_{2}λ_{3} + b_{1}λ_{3} + λ_{3} − λ_{1}λ_{2}}{λ_{2}(2 + b_{1})} F_{0}(b_{1} + 1, b_{2})
+ \frac{λ_{1}λ_{3}}{λ_{2}(2 + b_{1})}F_{0}(b_{1}, b_{2})

on b_{1} and

F_{0}(b_{1}, b_{2} + 2) = \frac{−λ_{3} + b_{1}λ_{3} − b_{2}λ_{3} + λ_{1}λ_{2}}{λ_{1}(b_{2} + 2)} F_{0}(b_{1}, b_{2} + 1)
+ \frac{λ_{2}λ_{3}}{λ_{1}(b_{2} + 2)}F_{0}(b_{1}, b_{2})

on b_{2}, with the following initial conditions:

( \begin{matrix}F_{0}(1, 2) & F_{0}(2, 2) \\ F_{0}(1, 1) & F_{0}(2, 1) \end{matrix} ) =

( \begin{matrix}λ_{3} + λ_{1}λ_{2} & λ_{2}λ_{3} + \frac{1}{2}λ_{2}^{2}λ_{1} \\ λ_{1}λ_{3} + \frac{1}{2}λ_{2}λ_{1}^{2} & \frac{1}{2}λ_{3}^{2} + λ_{2}λ_{1}λ_{3} + \frac{1}{4}λ_{2}^{2}λ_{1}^{2} \end{matrix} ) .

5.3 Another two-row example

The matrix A is:

A = ( \begin{matrix}1 & 0 & 1 \\ 0 & 2 & 1 \end{matrix} ) .  (5)

For the matrix A in (5), the “RecsV(A, λ, b)” command provides the following two-dimensional recurrence:

F_{0}(b_{1} + 3, b_{2}) = \frac{λ_{1}}{6 + 2b_{1}} F_{0}(b_{1} + 2, b_{2})
− \frac{(λ_{3}^{2} − b_{2}λ_{3}^{2} + 2λ_{2}λ_{1}^{2} + b_{1}λ_{3}^{2})}{2λ_{2}(3 + b_{1})(2 + b_{1})}F_{0}(b_{1} + 1, b_{2})
+ \frac{λ_{1}λ_{3}^{2}}{2λ_{2}(3 + b_{1})(2 + b_{1})}F_{0}(b_{1}, b_{2})

on b_{1} and

F_{0}(b_{1}, b_{2} + 3) = \frac{λ_{3}(b_{1} − 2 − b_{2})}{λ_{1}(3 + b_{2})}F_{0}(b_{1}, b_{2} + 2) + \frac{2λ_{2}}{b_{2} + 3}F_{0}(b_{1}, b_{2} + 1)
+ \frac{2λ_{2}λ_{3}}{λ_{1}(3 + b_{2})}F_{0}(b_{1}, b_{2})

on b_{2} with the initial conditions:

( \begin{matrix}F_{0}(1, 3) & F_{0}(2, 3) & F_{0}(3, 3) \\ F_{0}(1, 2) & F_{0}(2, 2) & F_{0}(3, 2) \\ F_{0}(1, 1) & F_{0}(2, 1) & F_{0}(3, 1) \end{matrix} ) =

( \begin{matrix}λ_{3} & λ_{1}λ_{2} & λ_{2}λ_{3} \\ λ_{1}λ_{3} & \frac{1}{2}λ_{3}^{2} + \frac{1}{2}λ_{2}λ_{1}^{2} & λ_{2}λ_{1}λ_{3} \\ \frac{1}{2}λ_{1}^{2}λ_{3} & \frac{1}{2}λ_{1}λ_{3}^{2} + \frac{1}{6}λ_{2}λ_{1}^{3} & \frac{1}{6}λ_{3}^{3} + \frac{1}{2}λ_{2}λ_{1}^{2}λ_{3} \end{matrix} ) .

See Figure 1.

[figure: 1; 2; 3; 4; 5; 1; 3; 4; 2]

Figure 1: Two-dimensional recursion fills-in the values of F_{0}(i, j) at the locations indicated by the open circles, using the initial data given at the locations indicated by the filled circles. For programming convenience, indices are positive integers: in the example shown, the initial conditions are specified for i, j = 1, 2, 3.

5.4 A two-row example with five columns

The matrix A is:

A = ( \begin{matrix}0 & 0 & 1 & 1 & 1 \\ 1 & 1 & 0 & 1 & 1 \end{matrix} ) .  (6)

For the matrix A in (6), the “RecsV(A, λ, b)” command provides the following two-dimensional recurrence:

F_{0}(b_{1} + 2, b_{2}) = − \frac{(−b_{2}λ_{5} − b_{2}λ_{4} + λ_{5} + λ_{4} − λ_{1}λ_{3} − λ_{2}λ_{3} + b_{1}λ_{5} + b_{1}λ_{4})}{(λ_{1} + λ_{2})(2 + b_{1})}F_{0}(b_{1} + 1, b_{2})
+ \frac{λ_{3}(λ_{5} + λ_{4})}{(λ_{1} + λ_{2})(2 + b_{1})}F_{0}(b_{1}, b_{2})

on b_{1} and

F_{0}(b_{1}, b_{2} + 2) = \frac{(−λ_{5} − λ_{4} + b_{1}λ_{5} + b_{1}λ_{4} − b_{2}λ_{5} − b_{2}λ_{4} + λ_{2}λ_{3} + λ_{1}λ_{3})}{λ_{3}(b_{2} + 2)}F_{0}(b_{1}, b_{2} + 1)
+ \frac{(λ_{5} + λ_{4})(λ_{1} + λ_{2})}{λ_{3}(b_{2} + 2)}F_{0}(b_{1}, b_{2})

on b_{2}, with the initial conditions:

( \begin{matrix}F_{0}(1, 2) & F_{0}(2, 2) \\ F_{0}(1, 1) & F_{0}(2, 1) \end{matrix} ) =

( \begin{matrix}λ_{5} + λ_{4} + (λ_{1} + λ_{2})λ_{3} & (λ_{5} + λ_{4})(λ_{1} + λ_{2}) + \frac{1}{2}(λ_{1} + λ_{2})^{2}λ_{3} \\ λ_{3}(λ_{5} + λ_{4}) + \frac{1}{2}(λ_{1} + λ_{2})λ_{3}^{2} & \frac{1}{2}(λ_{5} + λ_{4})^{2} + (λ_{1} + λ_{2})λ_{3}(λ_{5} + λ_{4}) + \frac{1}{4}(λ_{1} + λ_{2})^{2}λ_{3}^{2} \end{matrix} ) .

The command “CorMf(A, λ, b)” provides the correlation matrix for the X_{i}’s subject to Ax = b and assuming that the parameters are λ. For the matrix A considered here we obtain, for example with λ = (1, 1, 1, 1, 1) and b = (5, 5), the following result:

( \begin{matrix}1.0 & −.3647053019 & .5636021195 & −.2407443460 & −.2407443460 \\ −.3647053019 & 1.0 & .5636021195 & −.2407443460 & −.2407443460 \\ .5636021195 & .5636021195 & 1.0 & −.4271530174 & −.4271530174 \\ −.2407443460 & −.2407443460 & −.4271530174 & 1.0 & −.6350805992 \\ −.2407443460 & −.2407443460 & −.4271530174 & −.6350805992 & 1.0 \end{matrix} ) .

Note the negative entry for the correlation between X_{1} and X_{2}. This corresponds to the fact that Y_{2} = X_{1} + X_{2} + X_{4} + X_{5} = 5, so increases in X_{1} should be expected to result in decreases in X_{2}. Similar interpretations apply to the other entries.

5.5 A two-row example with six columns

The matrix A is:

A = ( \begin{matrix}1 & 0 & 1 & 0 & 1 & 1 \\ 0 & 1 & 1 & 1 & 0 & 1 \end{matrix} ) .  (7)

For the matrix A in (7), the “RecsV(A, λ, b)” command provides the following two-dimensional recurrence:

F_{0}(b_{1} + 2, b_{2}) = − \frac{λ_{3} + λ_{6} − λ_{5}λ_{4} − λ_{5}λ_{2} + b_{1}λ_{3} + b_{1}λ_{6} − λ_{1}λ_{4} − λ_{1}λ_{2} − b_{2}λ_{3} − b_{2}λ_{6}}{(λ_{4} + λ_{2})(2 + b_{1})} F_{0}(b_{1} + 1, b_{2})
+ \frac{(λ_{3} + λ_{6})(λ_{5} + λ_{1})}{(λ_{4} + λ_{2})(2 + b_{1})}F_{0}(b_{1}, b_{2})

on b_{1}, and

F_{0}(b_{1}, b_{2} + 2) = \frac{b_{1}λ_{6} + b_{1}λ_{3} − b_{2}λ_{6} − b_{2}λ_{3} + λ_{1}λ_{4} + λ_{5}λ_{4} + λ_{1}λ_{2} + λ_{5}λ_{2} − λ_{6} − λ_{3}}{(b_{2} + 2)(λ_{5} + λ_{1})} F_{0}(b_{1}, b_{2} + 1)
+ (λ_{3} + λ_{6})(λ_{4} + λ_{2})/(b_{2} + 2)(λ_{5} + λ_{1})F_{0}(b_{1}, b_{2})

on b_{2}, with the initial conditions:

F_{0}(1, 2) = λ_{3} + λ_{6} + (λ_{4} + λ_{2})(λ_{5} + λ_{1})
F_{0}(2, 2) = (λ_{3} + λ_{6})(λ_{4} + λ_{2}) + \frac{1}{2}(λ_{4} + λ_{2})^{2}(λ_{5} + λ_{1})
F_{0}(1, 1) = (λ_{5} + λ_{1})(λ_{3} + λ_{6}) + \frac{1}{2}(λ_{4} + λ_{2})(λ_{5} + λ_{1})^{2}
F_{0}(2, 1) = \frac{1}{2}(λ_{3} + λ_{6})^{2} + (λ_{4} + λ_{2})(λ_{5} + λ_{1})(λ_{3} + λ_{6}) + \frac{1}{4}(λ_{4} + λ_{2})^{2}(λ_{5} + λ_{1})^{2} .

5.6 A three-row example

The matrix A is:

A = ( \begin{matrix}0 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 1 \\ 1 & 1 & 0 & 0 & 1 & 1 \end{matrix} ) .  (8)

For the matrix A in (8), the “RecsV(A, λ, b)” command provides the following two-dimensional recurrence:

F (b_{1} + 2, b_{2}) = − \frac{−b_{2}λ_{6} − b_{2}λ_{3} − λ_{1}λ_{2} − λ_{1}λ_{4} − λ_{5}λ_{2} − λ_{5}λ_{4} + b_{1}λ_{6} + b_{1}λ_{3} + λ_{6} + λ_{3}}{(λ_{2} + λ_{4})(2 + b_{1})} F (b_{1} + 1, b_{2})
+ \frac{(λ_{6} + λ_{3})(λ_{1} + λ_{5})}{(λ_{2} + λ_{4})(2 + b_{1})}F (b_{1}, b_{2})

on b_{1}, and

F (b_{1}, b_{2} + 2) = \frac{b_{1}λ_{3} + b_{1}λ_{6} − λ_{3} − λ_{6} + λ_{1}λ_{2} + λ_{5}λ_{2} + λ_{1}λ_{4} + λ_{5}λ_{4} − b_{2}λ_{3} − b_{2}λ_{6}}{(b_{2} + 2)(λ_{1} + λ_{5})} F (b_{1}, b_{2} + 1)
+ \frac{(λ_{6} + λ_{3})(λ_{2} + λ_{4})}{(b_{2} + 2)(λ_{1} + λ_{5})}F (b_{1}, b_{2})

on b_{2}, with the initial conditions:

F_{0}(1, 2) = λ_{6} + λ_{3} + (λ_{2} + λ_{4})(λ_{1} + λ_{5})
F_{0}(2, 2) = (λ_{6} + λ_{3})(λ_{2} + λ_{4}) + (1/2)(λ_{2} + λ_{4})^{2}(λ_{1} + λ_{5})
F_{0}(1, 1) = (λ_{1} + λ_{5})(λ_{6} + λ_{3}) + (1/2(λ_{2} + λ_{4}))(λ_{1} + λ_{5})^{2}
F_{0}(2, 1) = (1/2)(λ_{6} + λ_{3})^{2} + (λ_{2} + λ_{4})(λ_{1} + λ_{5})(λ_{6} + λ_{3}) + (1/4)(λ_{2} + λ_{4})^{2}(λ_{1} + λ_{5})^{2} .

5.7 A four-row example

The A matrix is:

A = ( \begin{matrix}1 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 1 & 1 \end{matrix} ) .  (9)

For 4-row matrices as this one, the package MVPoisson is not able to return recurrences in a reasonable amount of time. However, one can now use the generating functions directly to compute the relevant quantities of interest, except that it is no longer possible to treat large inputs.

The command “SipurD” is used to generate averages and variances (“SipurD2f” implements a more efficient algorithm specifically for matrices with two rows). For the matrix A in (9) and, for example, λ = (1, 1, 1, 1, 1, 1, 1, 1) we obtain that E[X_{1}|Y = b] ≈ 1.897 when b = (10, 10, 10, 10) and ≈ 2.813 when b = (20, 20, 20, 20) (the value may be obtained to arbitrary precision), and that the variance of X_{1} conditioned on this same b is ≈ 1.112 when b = (10, 10, 10, 10) and ≈ 1.379 when b = (20, 20, 20, 20). The program also guesses asymptotic formulas for these quantities as a function of the entries of b, and as such is a useful tool in research, suggesting possible general formulas that one could attempt to prove.

6 Biochemical applications

We now explain how the problem studied here arises in the context of systems described by chemical network theory, and in particular chemical kinetics. There are two fundamentally different ways to mathematically model chemical reactions. One of them is based on differential equations modeling, and the other one on stochastic models. Our problem arises from this second approach. However, to understand its interest, it is important to first discuss the differential equation case. Differential equation models are useful when the number of molecules is very large, so that a continuous approximation is appropriate.

Suppose that n “species” interact through a network of reactions. The term species is used to refer to the elementary objects participating in the interactions: in molecular biology, these are typically ions, atoms, or molecules; in population biology and ecology, they may represent distinct animal or plant populations, particular age groups, and so forth. It is natural to describe such a network by a system of n differential equations which constrains the time evolution of the concentrations (or average populations) of the various species. These sets of differential equations take the following general form:

ẋ = ΓR(x)

(dot indicates time derivative) where x = x(t) is an n-vector of species concentrations (nonnegative real numbers) and Γ is an n × m matrix, called the “stoichiometry” matrix, whose columns describe how many units of each species are created or destroyed by each of m possible reactions. The components of the m-vector R(x) quantify the reaction rates for each of the m reactions, as a function of the current concentrations as well as parameters (reaction constants) that reflect physical and chemical information.

Chemical reactions are often described by graphs whose nodes are the “complexes” (the species, or combinations of species, that participate in the reactions) and whose edges are labeled by reaction rate information. Often, a mass-action kinetics model is used, which means that the reaction rate is proportional to the product of the concentrations of the reactants, and only the proportionality constant, called the kinetic constant associated to the corresponding reaction, is displayed on an edge. There is a systematic and simple way to map graph descriptions to differential equations.

Some of the main results in chemical network theory were obtained by Horn, Jackson, and Feinberg (see [7, 8] and also [14] for an exposition using a somewhat different formalism). These results guarantee that solutions of the system of differential equations are well-behaved (stability of equilibria, uniqueness of equilibria modulo stoichiometric constraints), provided that certain structural properties are satisfied by the network. The main such theorem is valid for what are called complex balanced networks. A sufficient (though not necessary) condition for complex balancing is that the network be “weakly reversible” and have “deficiency zero”. The deficiency is computed as c − ℓ − r, where c is the number of complexes, r is the rank of the matrix Γ, and ℓ is the number of “linkage classes” (connected components of the reaction graph). Weak reversibility means that each connected component of the reaction graph must be strongly connected. We refer the reader to the citations for details on deficiency theory. Our examples are all complex balanced.

When the numbers of molecules are very small, as is sometimes the case in molecular biology, a discrete stochastic model may be more appropriate than a continuous differential equation model. Indeed, fluctuations cannot be ignored when dealing with genes (usually one or two copies), mRNA’s (in the tens), ribosomes and RNA polymerases (up to hundreds) or certain proteins that have low concentrations.

Stochastic models fully account for the probabilistic nature of reactions. The number of individual copies of each species at (continuous) time t is viewed as a random process X_{i}(t), i = 1, …, n. The Chemical Master Equation (CME), which is the differential form of the Chapman-Kolmogorov forward equation, is a linear first-order differential equation that describes the time evolution of the joint probability distribution of the X_{i}(t)’s. Often, the interest is in long-time behavior, after a transient, that is to say in the probabilistic steady state of the system: the joint distribution of the random variables X_{i} = X_{i}(∞) that result in the limit as t → ∞ (provided that such a limit exists in an appropriate technical sense). This joint distribution is a solution of the steady state CME (ssCME), the infinite set of linear equations obtained by setting the right-hand side of the CME to zero.

A very beautiful recent observation, made in [2] (basically a rewording of classical results in queuing theory in Chapter 8 of [11], see also [12] for a discussion) is that the complex balancing condition, introduced originally for deterministic differential equation models, also guarantees that there is a solution π of the ssCME that is formally the joint distribution of n (the number of species) independent Poisson random variables. One associates to each deterministic steady state x̄ ∈ ℝ_{≥0}^{n} (that is, ΓR(x̄) = 0, in other words, a zero of the vector field ΓR(x)), a vector π that is a solution of the ssCME. The vector π is indexed by the n-dimensional lattice of nonnegative integers, N = (N_{1}, …, N_{n}) ∈ ℤ_{≥0}^{n}. We write the Nth entry of π as P (N) (thought of as a probability, in steady state, of the event (X_{1}, X_{2}, …, X_{n}) = (N_{1}, …, N_{n})). Let us write the product x̄_{1}^{N_{1}} … x̄_{n}^{N_{n}} as “x̄^{N}” and N_{1}! … N_{n}! as “N!”. Then, the assertion is that the vector π whose components are

P (N) = \frac{x̄^{N}}{N!}

(as well as any scalar multiple of this vector) is a solution of the ssCME. We provide a selfcontained proof of this fact in an Appendix to this paper.

However, the existence of this product form distribution does not mean that the joint distribution of the variables X_{i} will be independent Poisson, because the solution of the ssCME is not, in general, unique. The lack of uniqueness stems from conservation laws. Because of possible conservation laws, things are a bit subtle.

As an example, suppose that two molecules of species A and B can reversibly combine through a bimolecular reaction to produce a molecule of species C: A + B ↔ C. Let us denote the number of molecules of species A, B, and C at time t by X_{i}(t), i = 1, 2, 3 respectively. The count of A molecules goes down by one every time that a reaction takes place, at which time the count of C molecules goes up by one. Thus, the sum of the number of A molecules plus the number of C molecules remains constant: X_{1}(t) + X_{3}(t) = b_{1}. Similarly, X_{2}(t) + X_{3}(t) = b_{2}, because the total count of B and C molecules is also constant. This holds for all t, so taking limits as t → ∞ (ignoring technicalities!), we have that, for the steady state random variables, still X_{1} + X_{3} = b_{1} and X_{2} + X_{3} = b_{2}. Let us introduce Y_{1} = X_{1} + X_{3} and Y_{2} = X_{2} + X_{3}. Thus, depending on the initial conditions b_{1} = X_{1}(0) + X_{3}(0) and b_{2} = X_{2}(0) + X_{3}(0), the limiting distribution will be that of X_{1} and X_{2} conditioned on Y_{1} = b_{1} and Y_{2} = b_{2}. Once we collect this information into a matrix A, in this case

A = ( \begin{matrix}1 & 0 & 1 \\ 0 & 1 & 1 \end{matrix} ) ,

we are back to the situation where we want to study the behavior of the conditioned variables X_{i}|Y_{j}, where the X_{i}’s are Poisson distributed.^{∗} The rest of this section discusses various examples. To make the notations compatible with usage in probability theory, we use “λ” for the Poisson rates (instead of x̄_{i}) and k for multiindices (instead of N).

∗Our discussion is incomplete from a probabilistic viewpoint, as we have not addressed questions of uniqueness

6.1 A simple reversible reaction

Consider the following reaction:

X_{1} −⇀^{k_{1}} ; ↽−_{k_{2}} X_{2}  (10)

in which one molecule of substance X_{1} reversibly transforms to X_{2}.

This reaction system is complex-balanced, because it is weakly reversible and it has 2 complexes, 1 strongly connected component, and rank 1, and hence deficiency zero.

The steady states of this reaction network are given by the solutions λ = (λ_{1}, λ_{2}) of the equation k_{1}λ_{1} = k_{2}λ_{2}. We may pick, for example, λ = (1, k_{1}/k_{2}).

Every time that the forward reaction takes place, the count of molecules of X_{1} gets decreased by one and the count of molecules of X_{2} gets increased by one; the converse happens for the backward reaction. Thus, the total number of molecules of X_{1} and X_{2} remains constant. The corresponding A matrix is given in (3).

6.2 A bimolecular reaction

Consider the following reaction:

X_{1} + X_{2} −⇀^{k_{1}} ; ↽−_{k_{2}} X_{3}  (11)

in which one molecule of X_{1} combines reversibly with one molecule of X_{2} in order to produce one molecule of X_{3}.

This reaction system is complex-balanced, because it is weakly reversible and it has 2 complexes, 1 strongly connected component, and rank 1, and hence deficiency zero.

The steady states of this reaction network are given by the solutions λ = (λ_{1}, λ_{2}, λ_{3}) of the equation

k_{1}λ_{1}λ_{2} = k_{2}λ_{3} .

We may pick, for example, λ = (1, 1, k_{1}/k_{2}).

Every time that the forward reaction takes place, the counts of molecules of X_{1} and X_{2} each gets decreased by one and the count of molecules of X_{3} gets increased by one; the converse happens for the backward reaction. Thus, the total number of molecules of X_{1} and X_{3} remains constant, as does the total number of molecules of X_{2} and X_{3}. The matrix A is as in (4).

6.3 A more interesting bimolecular reaction

Consider the following reaction:

2X_{1} + X_{2} −⇀^{k_{1}} ; ↽−_{k_{2}} 2X_{3}  (12)

and convergence. These questions require a careful study of irreducibility properties of the associated Markov chains. We are only interested here in the computational problem of obtaining statistics for the conditioned variables.

which may represent, when X_{1} = H_{2}, X_{2} = O_{2}, and X_{3} = H_{2}O, the reversible creation of a molecule of water, when two molecules of the diatomic hydrogen gas combine with one molecule of the diatomic oxygen gas to produce two molecules of water. (The forward reaction produces energy, and the reverse reaction, breaking water to form hydrogen and oxygen, requires energy, for instance through electrolysis. The chemical reaction formalism used here does not account for energy production or consumption.)

This reaction system is complex-balanced, because it is weakly reversible and it has 2 complexes, 1 strongly connected component, and rank 1, and hence deficiency zero.

The steady states of this reaction network are given by the solutions λ = (λ_{1}, λ_{2}, λ_{3}) of the equation

k_{1}λ_{1}^{2}λ_{2} = k_{2}λ_{3}^{2} .

We may pick, for example, λ = (1, 1, \sqrt{k_{1}/k_{2}}).

The total sum of hydrogen and water molecules remains constant, and for each two molecules of oxygen there is one of water produced and viceversa. The matrix A is as in (5).

6.4 A receptor-ligand model

Receptor-ligand interactions play an important role in the understanding of the biochemical mechanisms that initiate cellular signaling, and their study is central to pharmacology. A “two-state” model for such interactions studied in [5] is shown, pictorially, in Figure 2.

[figure: C_{2}; k_{21}; k_{12}; R_{2} + L; k_{43}; k_{34}; k_{13}; k_{31}; k_{42}; k_{24}; C_{1}; R_{1} + L]

Figure 2: A two-state receptor-ligand network

The species participating in this reaction are: R_{1} and R_{2}, which represent the free receptors in an inactive and active conformational state respectively, the free ligand L, and the respective receptor-ligand complexes C_{1} = R_{1}L and C_{2} = R_{2}L.

The steady-states λ = (λ_{1}, λ_{2}, λ_{3}, λ_{4}, λ_{5}) = (R_{1}, R_{2}, L, C_{1}, C_{2}) of this system must satisfy the following polynomial equations:

−(k_{21} + k_{31})R_{1}L + k_{12}C_{1} + k_{13}R_{2}L = 0
−(k_{13} + k_{43})R_{2}L + k_{31}R_{1}L + k_{34}C_{2} = 0
−k_{21}R_{1}L − k_{43}R_{2}L + k_{12}C_{1} + k_{34}C_{2} = 0
−(k_{12} + k_{42})C_{1} + k_{21}R_{1}L + k_{24}C_{2} = 0
−(k_{34} + k_{24})C_{2} + k_{42}C_{1} + k_{43}R_{2}L = 0

For example, when all kinetic constants are k_{i} = 1 (this is not a realistic biological choice of constants, but is picked simply for illustration), then λ = (1, 1, 1, 1, 1) is a steady-state.

This reaction system is complex-balanced, because it is weakly reversible and it has 4 complexes, 1 strongly connected component, and rank 3, and hence deficiency zero.

The conservation of L + C_{1} + C_{2} (total amount of ligand) and R_{1} + R_{2} + C_{1} + C_{2} (total amount of receptors) leads to the matrix in (6).

6.5 A two-component signaling system in bacteria

The next example is from [4]. It models the “EnvZ/OmpR system” in E.coli bacteria. This system regulates the production of certain transport proteins (porins OmpF and OmpC) which act as pores allowing molecules to diffuse through the cell membrane. The system includes a kinase, EnvZ, which phosphorylates and dephosphorylates the response regulator OmpR, and is a particularly well-studied “two-component signaling system” in bacteria. The model is shown, pictorially, in Figure 3, where, for simplicity, we omit labeling each arrow by a reaction constant. We are using the following short-hand notations for the respective notations in [4]:

[figure: R + Z; R + ZP; ERP; EPR; RP + Z]

Figure 3: A two-component signaling system

X_{1} = R = OmpR, X_{2} = ZP = EnvZ-P (phosphorylated form), X_{3} = ERP = (EnvZ-P)OmpR (complex), X_{4} = Z = EnvZ, X_{5} = RP = OmpR-P (phosphorylated form), and X_{6} = EP R = (EnvZ)OmpR-P (complex).

This reaction system is complex-balanced, because it is weakly reversible and it has 5 complexes, 1 strongly connected component, and rank 4, and hence deficiency zero.

With all reaction constants equal to one, λ = (1, 1, 1, 1, 1, 1) is a steady state.

The system is described by six differential equations, subject to two constraints. These constraints reflect that the total amount of each of OmpR and EnvZ should stay constant, respectively, and give the rows of the A matrix for this example as that shown in (7).

6.6 A receptor antagonist model

The paper [10] analyzes a model involving the cytokine Interleukin-1 (IL-1), which is produced in response to inflammatory stimuli. The species in the model are IL-1 (denoted as L for “ligand”), the IL-1 receptor (denoted by R), the human IL-1 receptor antagonist (denoted by A), a decoy receptor or “trap” (denoted by T) which, by binding to the ligand, helps block IL-1 signaling, and the four possible dimers RL, RA, AT, and LT. The model consists of four reversible reactions:

R + L −⇀^{k_{2}} ; ↽−_{k_{1}} RL
R + A −⇀^{k_{4}} ; ↽−_{k_{3}} RA
A + T −⇀^{k_{6}} ; ↽−_{k_{5}} AT
L + T −⇀^{k_{8}} ; ↽−_{k_{7}} LT

This reaction system is complex-balanced, because it is weakly reversible and it has 8 complexes, 4 strongly connected components, and rank 4, and hence deficiency zero.

The total amounts of R, L, A, and T are conserved, giving rise to a matrix A with 4 rows. Ordering the states as follows: R, L, A, T, RL, RA, AT, LT, the resulting matrix is as in (9).

6.7 A futile-cycle example

We now describe an example that motivates looking at a matrix A as in (8). In contrast to the previous examples, however, this one is not complex-balanced and thus does not fit the assumptions for the ssCME having a solution in product form. So the interest in the conditional statistics problem for Poisson variables is purely academic for this particular example. Nonetheless, it is worth seeing how such a matrix A arises.

“Futile cycles” involving phosphorylation and dephosphorylation are ubiquitous in molecular biology (see for example [15] for more discussion and references). In such processes, an enzyme E (a kinase) catalizes the transformation of a substrate S into a product P, passing through one or more intermediate complexes C. A different enzyme F (a phosphotase) catalizes the transformation of P back into S, also passing through one or mode intermediate complexes. The simplest model (just one intermediate) for such a reaction is as follows:

E + S −⇀^{k_{1}} ; ↽−_{k_{2}} C −⇀^{k_{3}} ; ↽−_{k_{4}} E + P

F + P −⇀^{k_{5}} ; ↽−_{k_{4}} D −⇀^{k_{7}} ; ↽−_{k_{8}} F + S

in which we used C and D to denote the intermediate complexes. (Usually, the backward reactions to complex dissociation, labeled by k_{4} and k_{8}, are not included in the model, since they are energetically very unfavorable.) This system has deficiency one (6 complexes, two classes, and rank 3). Thus, the basic deficiency zero theory does not apply. Interestingly, however, a variation, “deficiency one theory”, can be used to predict the existance of multiple steady states for this system; see [6].

There are three conservation laws, corresponding to the conservation of total kinase, phosphotase, and substrate/product. Ordering the variables as S, P, E, F, C, D, we obtain the matrix A as in (8).

Appendix

For completeness, we show here that complex balanced reactions admit product form equilibrium densities for their Chemical Master Equations. The proof is basically that in [2, 11, 12].

Setup

A chemical reaction network is specified by:

R = {1, …, m}, the set of reactions.

C ⊆ ℝ_{≥0}^{n}, a finite set of complexes.

Example: if there are two reactions 1 : A + B → C + D and 2 : 2A + C → B, then the set C will have four elements, listing the species participating in each: (1, 1, 0, 0), (0, 0, 1, 1), (2, 0, 1, 0), (0, 1, 0, 0).

S, T : ℝ → C are the source and target functions that describe which are the reactant and product complexes, respectively.

Example: with the above reactions, S(1) = (1, 1, 0, 0), T (1) = (0, 0, 1, 1), S(2) = (2, 0, 1, 0), T (2) = (0, 1, 0, 0).

We make the following notational convention: for vectors x, c ∈ ℝ_{≥0}^{n}, x^{c} := x_{1}^{c_{1}} … x_{n}^{c_{n}} (with 0^{0} = 1), and for nonnegative integer vectors N = (N_{1}, …, N_{n}), N! := N_{1}! … N_{n}!.

By definition, a vector π = (P (N), N ∈ ℤ_{≥0}^{n}) is a steady-state solution of the Chemical Master Equation associated to a given reaction network if it satisfies:

∑_{i∈R} P (N − T (i) + S(i)) A_{i}(N − T (i) + S(i)) = ∑_{i∈R} P (N) A_{i}(N)  (13)

for each N ∈ ℤ_{≥0}^{n}, where A_{i}(N) is the ith “propensity function” [9]: A_{i}(N)dt is the probability that reaction i will occur in a small time interval [t, t+dt] if the state of the system is N at time t. This function is proportional to the number of ways in which the N molecules can combine to form the ith complex:

A_{i}(N) = k_{i} \frac{N!}{(N − S(i))!} .

The constant k_{i} is related to the deterministic kinetic constant of the respective reaction through division by a power of the volume in which the reaction takes place.

A complex balanced steady state (CBSS) with respect to the given network and kinetic constants k is an x̄ ∈ ℝ_{>0}^{n} (which is thought of as a vector of species concentrations) such that the following property holds for each complex c ∈ C:

∑_{i∈T^{−1}(c)} k_{i}x̄^{S(i)} = ∑_{i∈S^{−1}(c)} k_{i}x̄^{S(i)}  (14)

(note that one can equally well write “x̄^{c}” and bring this term outside of the sum, in the right-hand side).

Complex balancing means that each “complex” is balanced in inflow and outflow. This is a Kirschoff current law (in-flux = out-flux, at each node) when one writes a chemical network.

A counter-example to complex-balancing is this reaction network:

A ⟶^{k_{1}} B, 2B ⟶^{k_{2}} 2A

(or, if one prefers reversible reactions, one may take instead an example due to Wegsheider, A ↔ B and 2A ↔ B). In steady state, k_{1}a − 2k_{2}b^{2} = 0. But complex-balancing would require that the outflow of “A” be zero (since there are no inflows into the “complex” A), which means k_{1}a = 0, and misses the nonzero steady states. (One could also argue with the complex 2A, or with B, or with 2B.)

A complex-balanced system is one with the property that every steady state is complex balanced. This concept was studied in detail by Horn and Jackson and by Feinberg in the early 1970s. Feinberg [7, 8] showed that for a special type of system (weakly reversible and deficiency zero), for any kinetic constants there is a steady state x̄ ∈ ℝ_{>0}^{n} satisfying (14) (and, in fact, every other steady state will also be complex balanced). Let us call a system with these properties a “Feinberg-like system”.

The Key Lemma

Suppose that the reaction network is a Feinberg-like system. Let x̄ satisfy (14), where the k’s are the proportionality factors in (13). Take any function α : C → ℝ on complexes. Then:

∑_{i∈R} k_{i}x̄^{S(i)−T (i)}α(T (i)) = ∑_{i∈R} k_{i}α(S(i)) .  (15)

Proof. Since

∑_{i∈R} = ∑_{c∈C} ∑_{i∈T^{−1}(c)} and ∑_{i∈R} = ∑_{c∈C} ∑_{i∈S^{−1}(c)}

it is enough to show that, for each fixed c:

∑_{i∈T^{−1}(c)} k_{i}x̄^{S(i)−c} α(T (i)) = ∑_{i∈S^{−1}(c)} k_{i} α(S(i)) .

Since T (i) = c and S(i) = c in the left-hand side and right-hand side respectively, this is the same as the CBSS condition upon multiplication by x̄^{−c}α(c).

Corollary. The vector Π with

P (N) = \frac{x̄^{N}}{N!}

is a steady state solution of the CME.

Proof. Obvious using α(c) = \frac{x̄^{N}}{(N−c)!}.

References

[1] G. Almkvist and D. Zeilberger, The method of differentiating under the integral sign, J. Symbolic Computation 10(1990), 571-591.

[2] D.F. Anderson, G. Craciun, and T.G. Kurtz. Product-form stationary distributions for deficiency zero chemical reaction networks, 2008. arXiv.org:0803.3042.

[3] M. Apagodu and D. Zeilberger, Multi-Variable Zeilberger and Almkvist-Zeilberger Algorithms and the Sharpening of Wilf-Zeilberger Theory , Adv. Appl. Math. 37 (2006)(Special Regev issue), 139-152

[4] E. Batchelor and M. Goulian. Robustness and the cycle of phosphorylation and dephosphorylation in a two-component regulatory system. Proc. Natl. Acad. Sci. U.S.A., 100:691–696, 2003.

[5] M. Chaves, E.D. Sontag, and R. J. Dinerstein. Steady-states of receptor-ligand dynamics: A theoretical framework. J. Theoret. Biol., 227(3):413–428, 2004.

[6] C. Conradi, J. Saez-Rodriguez, E.-D. Gilles, and J. Raisch. Using chemical reaction network theory to discard a kinetic mechanism hypothesis. IEE Proceedings Systems Biology, 152:243 – 248, 2005.

[7] M. Feinberg. Chemical reaction network structure and the stability of complex isothermal reactors - i. the deficiency zero and deficiency one theorems. Chemical Engr. Sci., 42:2229– 2268, 1987.

[8] M. Feinberg. The existence and uniqueness of steady states for a class of chemical reaction networks. Archive for Rational Mechanics and Analysis, 132:311–370, 1995.

[9] D. T. Gillespie. The chemical Langevin equation. Journal of Chemical Physics, 113(1):297– 306, 2000.

[10] G. Gnacadja, A. Shoshitaishvili, M.J. Gresser, B. Varnum, D. Balaban, M. Durst, C. Vezina, and Y. Li. Monotonicity of interleukin-1 receptor-ligand binding with respect to antagonist in the presence of decoy receptor. Journal of Theoretical Biology, 244:478–488, 2007.

[11] F. Kelly. Reversibility and Stochastic Networks. Wiley, New York, 1979.

[12] J. Mairesse and H.-T. Nguyen. Deficiency zero petri nets and product form, 2009. arXiv.org:0905.3158.

[13] M. Petkovsek, H.S. Wilf and D. Zeilberger, A=B, AK Peters, Wellesley, (1996). [available on-line from the authors’ websites.]

[14] E.D. Sontag. Structure and stability of certain chemical networks and applications to the kinetic proofreading model of T-cell receptor signal transduction. IEEE Trans. Automat. Control, 46(7):1028–1047, 2001.

[15] L. Wang and E.D. Sontag. Singularly perturbed monotone systems and an application to double phosphorylation cycles. J. Nonlinear Science, 18:527–550, 2008.

[16] H.S. Wilf and D. Zeilberger, An algorithmic proof theory for hypergeometric (ordinary and “q”) multisum/integral identities, Invent. Math. 108 (1992), 575-633.

[17] D. Zeilberger, The method of creative telescoping, J. Symbolic Computat. 11, 195-204 (1991).