# A general method for the computation of probabilities in

## Transcript Of A general method for the computation of probabilities in

THE JOURNAL OF CHEMICAL PHYSICS 122, 104101 ͑2005͒

A general method for the computation of probabilities in systems of ﬁrst order chemical reactions

Xueying Zhang,a͒ Katrien De Cock,b͒ Mónica F. Bugallo,c͒ and Petar M. Djurićd͒ COSINE Laboratory, Department of Electrical & Computer Engineering, Stony Brook University, Stony Brook, New York 11794-2350

͑Received 21 June 2004; accepted 9 December 2004; published online 7 March 2005͒

We present a general method for the computation of molecular population distributions in a system of ﬁrst-order chemical reactions. The method is based on the observation that the molecules in ﬁrst-order reactions do not interact with each other. Therefore, the population distributions are a convolution of densities for one molecule. With this method one can study reactions involving any number of molecules. Such analysis is demonstrated on several examples, including an enzyme catalyst reaction and a ﬁrst-order reaction chain. © 2005 American Institute of Physics. ͓DOI: 10.1063/1.1855311͔

I. INTRODUCTION

In small biological systems, it is often more appropriate to model the chemical reactions in a stochastic way rather than with the traditional differential equations for evolution of concentrations. The reasons are ͑a͒ the number of molecules involved in a biological system can be very small;1 ͑b͒ there are many situations in which molecules are not homogeneously distributed;2,3 ͑c͒ the change of molecular population levels is a discrete integer amount; and ͑d͒ the population change is not a deterministic process.4,5

There are mainly two approaches for the stochastic study of the number of molecules in biochemical reactions: the ﬁrst is based on the analysis of the master equation6–8 and the second relies on Monte Carlo simulation methods.9–11 Here, we adopt the ﬁrst approach.

In 1967, McQuarrie summarized the stochastic study using the master equation approach.7 For a reaction, it is assumed that in an inﬁnitesimal time interval, the probability of having one reaction per unit reactant molecule combination is proportional to the length of the time interval. A probability difference function is ﬁrst obtained based on the assumption, which leads to a differential-difference equation called the master equation of the reaction. The momentgenerating function is then employed to transform the master equation into a partial differential equation. It can sometimes be solved to obtain the analytic solution. More often, the mean and variance of the number of molecules can only be obtained.

In 2000, Laurenzi introduced a different way of solving the master equation.8 Instead of using the momentgenerating function, he applied the Laplace transform. In this way, solving the partial differential equations is avoided, which is important for more complicated reactions. Instead, one needs to solve a set of linear equations. We will show

a͒Electronic mail: [email protected] b͒Electronic mail: [email protected] c͒Electronic mail: [email protected] d͒Electronic mail: [email protected]

how to obtain the solution of the master equation for any system of ﬁrst-order reactions in yet another way.

In this paper we study systems of ﬁrst-order reactions, i.e., combinations of ͑coupled͒ reactions of the type Si → Sj. Because in such systems molecules do not chemically interact, the analysis is much simpler than when second-order reactions are involved. Although this situation might occur only rarely in practice, our analysis is also useful for systems where the ﬁrst-order reactions are dominant and for systems where second-order reactions can be approximated by ﬁrstorder reactions ͑see Sec. IV for an example͒.

The independence of the molecules in a ﬁrst-order reaction system is exploited to derive the population distributions of the molecules in the system. Instead of solving the master equation for the complete system, where the number of states grows with the number of molecules, we ﬁrst solve the master equation for only one molecule and use its solution to ﬁnd the population distributions for all species.

The structure of the paper is as follows. In Sec. II we solve the master equation for one molecule in a ﬁrst-order reaction system. The resulting probability distribution functions are used in Sec. III to obtain the population distributions when more molecules are introduced. In Sec. IV we illustrate the general results of Sec. III with speciﬁc examples. Finally, in Sec. V we give some concluding remarks.

II. ONE MOLECULE

In this section we study the master equation for one molecule in a system of ﬁrst-order reactions. The solution will then be used to solve more general cases in Sec. III.

Assume that there are M molecule species, S1 , S2 , . . . , SM in the reaction system, and let the probability rate constant of the reaction Si → Sj be denoted by cij. For a ﬁrst-order reaction the speciﬁc probability rate constant is such that

cij⌬t͑i j͒ ϵ the probability for the first-order reaction

0021-9606/2005/122͑10͒/104101/8/$22.50

122, 104101-1

© 2005 American Institute of Physics

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-2 Zhang et al.

J. Chem. Phys. 122, 104101 ͑2005͒

cij

Si→Sj to happen per molecule in an infinitesimally small time interval ⌬t.

The master equation describes the evolution of the molecules’ population distribution in the chemical system as a function of time. When there is only one molecule in a ﬁrstorder reaction system, the master equation gives the time evolution of the probability that this molecule has become a certain chemical species. Let p͑i͒͑t͒ be the probability that the molecule is an Si molecule at time t. For example, in the simple system of reactions

c12 c23

S1 S2→S3,

͑1͒

c21

the differential equation describing the time evolution of p͑2͒ is the following ﬁrst-order linear differential equation:

dpd͑2t͒͑t͒ = c12 p͑1͒͑t͒ − ͑c21 + c23͒p͑2͒͑t͒. ͑2͒

We assume in the rest of this paper that the initial values for the probabilities are equal to zero: p͑i͒͑0͒ = 0, i = 1 , . . . , M. The injection of the molecule in the system is modeled by a source probability density function. For the given example, if the molecule begins as an S2 molecule, the differential equation becomes

dpd͑2t͒͑t͒ = c12 p͑1͒͑t͒ − ͑c21 + c23͒p͑2͒͑t͒ + f͑2͒͑t͒, ͑3͒

where f ͑2͒͑t͒ is the source’s probability density, i.e., f ͑2͒͑t͒⌬t is the probability that a molecule from the source is injected in state S2 in the interval ͓t , t + ⌬t͔. If the injection time is known precisely, e.g., t = , the source’s probability density function is a delta function: f ͑2͒͑t͒ = ␦͑t − ͒. Otherwise, the density function for the injection time of the molecule can have forms such as ae−a͑t−͒u͑t − ͒, a␦͑t − 1͒ + ͑1 − a͒␦͑t − 2͒, and

͚ϱ ne−

␦͑t − n͒. n=0 n!

In general, the system can be represented by a Markov chain with a transition matrix

͚ c11 c12 ¯ c1M

c21 c22 ¯ c2M

M

C=

, where cii = −

ci j ,

]]

]

j=1,j i

cM1 cM2 ¯ cMM

͑4͒

and the set of ﬁrst-order linear differential equations has the following form:

dp͑1͒͑t͒

dt

p͑1͒͑t͒

f ͑1͒͑t͒

dp͑2͒͑t͒

p͑2͒͑t͒

f ͑2͒͑t͒

dt

= CT

+

,

͑5͒

]

]

] dp͑M͒͑t͒

p͑M͒͑t͒

f ͑M͒͑t͒

dt

where

M

͚ p͑i͒͑t͒ = 1,

i=1

and

Ά͵ · ϱ f ͑i͒͑t͒dt = 1 t=0

if the molecule from source is injected

into state i, f ͑i͒͑t͒ = 0 otherwise.

The equations can be solved by ﬁrst applying the Laplace transform, followed by solving the resulting algebraic equation, and ﬁnally using the inverse Laplace transform. More speciﬁcally, after taking the Laplace transform on both sides of ͑5͒, we obtain

sp = CTp + f,

͑6͒

where p is an M ϫ 1 vector whose elements p͑i͒͑s͒ are the Laplace transforms of p͑i͒͑t͒, the sought-after probability

density functions, and f is an M ϫ 1 vector with elements f ͑i͒͑s͒, which are the Laplace transforms of the source functions f ͑i͒͑t͒. Note that f has only one nonzero element.

After solving for p in ͑6͒, we get

p = ͑sI − CT͒−1f.

͑7͒

Let

L = ͑sI − CT͒−1

͑8͒

and

G = L−1͑L͒,

͑9͒

where L−1 is the inverse Laplace operator. If the elements of G are denoted by gij͑t͒, and f ͑m͒͑s͒ is the nonzero element of f, the solutions p͑i͒͑t͒ are given by

p͑i͒͑t͒ = gim͑t͒ f͑m͒͑t͒, i = 1,2, . . . , M ,

͑10͒

where the symbol denotes convolution. Finding the inverse Laplace transforms of the elements of L may be tedious but it is straightforward.

This method is illustrated in details by example 1. As already mentioned, the technique of using the Laplace transform for solving the master equation was introduced by Laurenzi in Ref. 8. The difference between his and our approach is that Laurenzi was solving the many-molecule master equation, while here we only solve the one-molecule master equation.

In the deterministic framework, the concentrations of the species are found by solving the same differential equations

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-3 Probabilities in systems of ﬁrst-order reactions

͑5͒, without the stochastic sources, but taking into account the initial concentrations. Furthermore, the so-obtained results there yield only the means of the stochastic distributions.

III. MORE MOLECULES

In this section we treat the case when several molecules are injected into the system. First, in Sec. III A, we assume that all molecules are inserted by the same source, so that they all start as the same chemical species. In Sec. III B, the molecules can originate from different sources.

A. One source

Suppose that there are x molecules of species Sk injected into the system with the same probability density f ͑k͒͑t͒, and no molecules of other species. Denote the number of molecules at state Si ͑i = 1 , . . . , M͒ at time t by N͑i͒͑t͒, and let K͑ji͒͑t͒ = 1 if molecule j is at state Si at time t and K͑ji͒͑t͒ = 0 if it is not. Thus, N͑i͒͑t͒ = ⌺xj=1K͑ji͒͑t͒. Note that K1͑i͒͑t͒, K2͑i͒͑t͒ , . . . , Kx͑i͒͑t͒ are independent Bernoulli random variables with the same probability p͑i͒͑t͒, which can be computed by solving ͑5͒. Therefore, N͑i͒͑t͒ has a binomial distribution. The probability that there are n molecules in state Si at time t is

ͩ ͪ Pn͑i͒͑t͒ = nx ͑p͑i͒͑t͒͒n͑1 − p͑i͒͑t͒͒x−n, ͑11͒

with mean i͑t͒ and variance i2͑t͒ equal to i͑t͒ = xp͑i͒͑t͒,

i2͑t͒ = xp͑i͒͑t͒͑1 − p͑i͒͑t͒͒,

respectively. Moreover, the joint probability of ﬁnding ni molecules in state Si ͑i = 1 , . . . , M͒ at time t is equal to

P͑n , . . . ,n ;t͒ = x! ͑p͑1͒͑t͒͒n1 ͑p͑M͒͑t͒͒nM,

1

M n1! ¯ nM!

¯

͑12͒

where ⌺iM=1ni = x. Equation ͑11͒ is one of the marginals of ͑12͒.

When the number of molecules is large, i.e., when xp͑i͒͑t͒͑1 − p͑i͒͑t͒͒ ӷ 1, the binomial distribution ͑11͒ can be

well approximated by a Gaussian distribution, or

ͱ Pn͑i͒͑t͒ Ӎ

1 e−͓͑n − i͑t͒͒2/2i2͑t͔͒ . 2i2͑t͒

The purpose for employing the Gaussian approximation is to save computation time when the number of molecules is very large, as well as to make the analysis more tractable. Numerical examples of the comparison between analytic distributions and their Gaussian approximations can be found in Ref. 12.

The same purpose can be achieved if the binomial distribution ͑11͒ is approximated by a Poisson distribution. Conditions for the validity of such approximation are stricter. First, the number of molecules should be large. Second, the probability for a single molecule in state i should be small

J. Chem. Phys. 122, 104101 ͑2005͒

͓e.g., p͑i͒͑t͒ Ͻ 0.1͔. Finally, the product of x and p͑i͒͑t͒ should

be a sizable number. Under such conditions, the variance of the binomial distribution roughly equals the mean. Then ͑11͒ can be approximated by the Poisson distribution with parameter i = xp͑i͒͑t͒, or

P͑i͒͑t͒ Ӎ in e−i.

n

n!

B. More sources

If there are several sources which are independent from

each other, the probability mass function for state Si is the convolution of binomial distributions. Suppose there are J

independent sources. Denote the number of molecules at state Si ͑i = 1 , . . . , M͒ at time t by N͑i͒͑t͒, and denote the number of molecules at state Si that are injected from source j by K͑ji͒͑t͒ ͑j = 1 , . . . , J͒. Then, N͑i͒͑t͒ = ⌺Jj=1K͑ji͒͑t͒. Since K͑ji͒͑t͒ ͑j = 1 , . . . , J͒ are independent random variables, the probability mass function of the random variable N͑i͒͑t͒ is the con-

volution of the probability mass functions of the random variables K͑ji͒͑t͒.13

For example, let x Sk molecules be injected by the source functions f ͑k͒͑t͒ and y Sl molecules by f ͑l͒͑t͒. Each of the x Sk molecules gives rise to the state probabilities pk͑i͒͑t͒ ͑i = 1 , . . . , M͒, whereas a molecule introduced by f ͑l͒͑t͒, to p͑i͒͑t͒ ͑i = 1 , . . . , M͒. Then, the number of molecules at state Sil that are injected from source k is Kk͑i͒͑t͒, which has a binomial distribution ͑ nx ͒͑pk͑i͒͑t͒͒n͑1 − pk͑i͒͑t͒͒x−n. Similarly, the number of molecules at state Si that are injected from source l is Kl͑i͒͑t͒, which has a binomial distribution ͑ ny ͒͑pl͑i͒͑t͒͒n͑1 − pl͑i͒͑t͒͒y−n. The probability mass function of N͑i͒͑t͒, the total number of molecules at state Si is equal to the convolution of two binomial distributions, i.e.,

ͩ ͪ Pn͑i͒͑t͒ = nx ͑pk͑i͒͑t͒͒n͑1 − pk͑i͒͑t͒͒x−n

ͩ ͪ ny ͑pl͑i͒͑t͒͒n͑1 − pl͑i͒͑t͒͒y−n, ͑13͒

where the convolution is along the n axis, that is

ϱ

͚ y1͑n͒ y2͑n͒ = y1͑m͒y2͑n − m͒. m=0

The mean i͑t͒ and variance i2͑t͒ of N͑i͒͑t͒ are

i͑t͒ = xpk͑i͒͑t͒ + ypl͑i͒͑t͒,

͑14͒

i2͑t͒ = xpk͑i͒͑t͒͑1 − pk͑i͒͑t͒͒ + ypl͑i͒͑t͒͑1 − pl͑i͒͑t͒͒,

͑15͒

respectively. As we mentioned before, when the number of molecules

is large, the binomial distribution can be well approximated by a Gaussian distribution with the same mean and variance. One property of Gaussian distributions is that the convolution of two Gaussian distributions yields a Gaussian distribution. Therefore, the convolution result in ͑13͒ can be approximated by a Gaussian distribution with mean ͑14͒ and variance ͑15͒. The Poisson distribution has the same prop-

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-4 Zhang et al.

erty: the convolution of Poisson distributions with parameters ␣ and , respectively, yields a Poisson distribution with parameter = ␣ + .

IV. EXAMPLES

In this section several examples of systems of ﬁrst-order reactions are discussed. For some systems, an analytic solution for the population distribution is derived. This section also includes two numerical examples: the ﬁrst is a onedirectional reaction chain and the second an enzyme substrate reaction.

A. Example 1: The ﬁrst-order irreversible reaction

c

Consider the ﬁrst-order reaction S1→S2, where only

molecules S1 are injected into the system with probability density function f ͑1͒͑t͒ = ␦͑t͒. The ﬁrst step is to rewrite ͑5͒ with M = 2, c12= c, c21= 0, and f ͑2͒͑t͒ = 0

ͩ ͪ ͩ ͪ ͩ ͪ dp͑1͒͑t͒

dt

− c c T p͑1͒͑t͒ ␦͑t͒

dp͑2͒͑t͒ = 0 0 p͑2͒͑t͒ + 0 .

͑16͒

dt

Taking the Laplace transform leads to

ͩ ͪ ͩ ͪ ͩ ͪ ͩ ͪ p͑1͒͑s͒

− c c T p͑1͒͑s͒

1

s p͑2͒͑s͒ = 0 0 p͑2͒͑s͒ + 0 .

The matrix L from ͑8͒ is given by

1

0

s+c

L=

.

1 11

−

s s+c s

After taking the inverse Laplace transform of L, we obtain

ͩ ͪ e−ct

0

G = 1 − e−ct u͑t͒ ,

where u͑t͒ is the unit step function. The ﬁnal solution is

p͑1͒͑t͒ = e−ctu͑t͒,

͑17͒

p͑2͒͑t͒ = ͑1 − e−ct͒u͑t͒.

͑18͒

In the rest of the paper we will drop the use of u͑t͒, and so

we point out that the given solutions for p͑t͒ are valid only for t ജ 0.When f ͑1͒͑t͒ is some other density function

͵ p͑1͒͑t͒ =

ϱ

f ͑1͒͑͒e−c͑t−͒d,

0

͵ p͑2͒͑t͒ =

ϱ

f ͑1͒͑͒͑1 − e−c͑t−͒͒d.

0

Note that if an S2 instead of an S1 molecule was injected into the system with f ͑2͒͑t͒, p͑1͒͑t͒ = 0 and p͑2͒͑t͒ = ͐0t f ͑2͒͑͒d.

By substituting the probabilities ͑17͒ and ͑18͒ in ͑11͒, we obtain the well-known result for f ͑͑t1͒͒ = ␦͑t͒

J. Chem. Phys. 122, 104101 ͑2005͒

ͩ ͪ Pn͑1͒͑t͒ = nx e−nct͑1 − e−ct͒x−n, ͩ ͪ Pn͑2͒͑t͒ = nx ͑1 − e−ct͒ne−c͑x−n͒t.

B. Example 2. The ﬁrst-order reversible reaction

c1

Consider the ﬁrst-order reversible reaction S1 S2 which

c2

starts with x S1 molecules and y S2 molecules. Each of the x

S1 molecules gives rise to the following probabilities, which are obtained by solving ͑5͒ with M = 2, f ͑1͒͑t͒ = ␦͑t͒, and f ͑2͒͑t͒ = 0:

p͑1͒͑t͒ = 1 ͑c + c e−͑c1+c2͒t͒,

1

c1 + c2 2 1

p͑2͒͑t͒ = c1 ͑1 − e−͑c1+c2͒t͒.

1

c1 + c2

͑19͒

For the effect of each of the y S2 molecules, we take f ͑1͒͑t͒ = 0 and f͑2͒͑t͒ = ␦͑t͒ and obtain

p͑1͒͑t͒ = c2 ͑1 − e−͑c2+c1͒t͒,

2

c2 + c1

p͑2͒͑t͒ = 1 ͑c + c e−͑c2+c1͒t͒.

2

c2 + c1 1 2

͑20͒

By substituting the probabilities of ͑19͒ and ͑20͒ in ͑13͒, we

obtain the distribution for the S1 molecules.

ͩ ͪͩ ͪ Pn͑1͒͑t͒ = nx

1

n

͑c + c e−͑c1+c2͒t͒

c1 + c2 2 1

ͩ ͪ ϫ 1 −

1

x−n

͑c + c e−͑c1+c2͒t͒

c1 + c2 2 1

ͩ ͪͩ ͪ y

n

c2

n

͑1 − e−͑c1+c2͒t͒

c2 + c1

ͩ ͪ ϫ 1 −

c2

y−n

͑1 − e−͑c1+c2͒t͒ ,

c2 + c1

͑21͒

and similarly for the S2 molecules. Note that when only S1 or only S2 molecules are present

at t = 0, the distribution reduces to the ﬁrst and second convolution factor in ͑21͒, respectively. The population distribu-

tion for this special case has already been given by McQuarrie.7

Also note that deriving the mean and variance for the

number of S1 and S2 molecules is straightforward, as explained in Sec. III B, and leads to the same results as re-

ported in Refs. 14 and 15.

C. Example 3: A combination of the ﬁrst-order irreversible and ﬁrst-order reversible reaction

Consider the reaction system

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-5 Probabilities in systems of ﬁrst-order reactions

J. Chem. Phys. 122, 104101 ͑2005͒

FIG. 1. Direct calculation of the probability mass function of S7 at time t = 1.5 s.

c1 c3

S1 S2→S3,

c2

͑22͒

and assume there are x S1 molecules and y S2 molecules at t = 0. Again, we ﬁrst solve ͑5͒ with f ͑1͒͑t͒ = ␦͑t͒, f ͑2͒͑t͒ = f ͑3͒͑t͒ = 0 to ﬁnd the probabilities due to one S1 molecule

and obtain

ͱ ͱ p1͑1͒͑t͒ =

m − c1 + c2 + c3 e−͑1/2͒͑−ͱm+c1+c2+c3͒t 2m

ͱ + m + c1 − c2 − c3 e−͑1/2͒͑ͱm+c1+c2+c3͒t, ͱ2 m

p1͑2͒͑t͒

=

ͱc1 e−͑1/2͒͑−ͱm+c1+c2+c3͒t

m

−

ͱc1 e−͑1/2͒͑ͱm+c1+c2+c3͒t,

m

͑23͒

FIG. 3. The probabilities for a single molecule to be in state Si ͑i = 1 , 2 , . . . , 7͒ at time t.

ͱ ͱ p1͑3͒͑t͒ = 1 −

m + c1 + c2 + c3 e−͑1/2͒͑−ͱm+c1+c2+c3͒t 2m

ͱ − m − c1 − c2 − c3 e−͑1/2͒͑ͱm+c1+c2+c3͒t, ͱ2 m

where

m = c12 + 2c1c2 − 2c1c3 + c22 + 2c2c3 + c32.

Next, we solve ͑5͒ with f ͑1͒͑t͒ = f ͑3͒͑t͒ = 0, f ͑2͒ = ␦͑t͒ to ﬁnd the probabilities due to one S2 molecule and ﬁnd that

p2͑1͒͑t͒

=

ͱc2 e−͑1/2͒͑−ͱm+c1+c2+c3͒t

m

−

ͱc2 e−͑1/2͒͑ͱm+c1+c2+c3͒t,

m

ͱ ͱ p2͑2͒͑t͒ =

m + c1 − c2 − c3 e−͑1/2͒͑−ͱm+c1+c2+c3͒t 2m

ͱ + m − c1 + c2 + c3 e−͑1/2͒͑ͱm+c1+c2+c3͒t, ͱ2 m

͑24͒

ͱ ͱ p2͑3͒͑t͒ = 1 −

m + c1 + c2 − c3 e−͑1/2͒͑−ͱm+c1+c2−c3͒t 2m

ͱ − m − c1 − c2 + c3 e−͑1/2͒͑ͱm+c1+c2+c3͒t. ͱ2 m

The substitution of the probabilities in ͑13͒ by the probabilities in ͑23͒ and ͑24͒ yields the population distribution for the S1, S2, and S3 molecules.

FIG. 2. The probability mass of S7 at time t = 1.5 s calculated using the SSA, with 500 realizations.

D. Example 4: A ﬁrst-order reaction chain

Consider a ﬁrst-order reaction chain in its general form

c1 c2

ci−1 ci

cN−1

S1→S2→ ¯ → Si→ ¯ → SN.

͑25͒

Assume that ci cj when i j ͑i = 1 , . . . , N and j = 1 , . . . , N͒ and that x molecules enter the system in state S1 at time t = 0. The input probability distribution of state S1 is f ͑1͒͑t͒

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-6 Zhang et al.

J. Chem. Phys. 122, 104101 ͑2005͒

FIG. 4. ͑a͒ A cell on which enzyme molecules are immobilized. The substrate molecules ﬂoat around freely. ͑b͒ Illustration of the chemical reaction. The substrate molecule AB approaches the enzyme molecule R ͑left͒. The enzyme recognizes the A molecule of the AB structure ͑middle͒. The enzyme destroys the bond between A and B, A remains attached to R, and B is released into the solution ͑right͒.

= ␦͑t͒. Then, using the method introduced in the previous sections, the state probabilities for one molecule are equal to

Ά ͚͚ · p͑i͒͑t͒ =

i

1 M͑1,i͒c e−cjt, ci j=1 j j

ci Ͼ 0

i−1

1 − M͑j1,i−1͒e−cjt,

j=1

ci = 0,i ജ 2,

1, ci = 0,i = 1

where

Ά ͟ M͑pk,m͒ =

m

cj ,

j=k,j p cj − cp 1, m = k.

mϾk·

The probability that n of the x molecules are in state Si at time t is given by

FIG. 5. Mean for the population of molecule species R ͑full line͒, ARB ͑dot-dashed line͒, and B ͑dashed line͒.

FIG. 6. Zoomed plots of the mean for the population of molecule species R ͑full line͒, ARB ͑dot-dashed line͒, and B ͑dashed line͒.

ͩ ͪ Pn͑i͒͑t͒ = nx ͑p͑i͒͑t͒͒n͑1 − p͑i͒͑t͒͒x−n. ͑26͒

If there are ci’s equal to cj’s, a closed-form solution is also possible.

E. Example 5: A numerical example for the ﬁrst-order reaction chain

Consider the ﬁrst-order reaction chain

c1 c2 c3 c4 c5 c6

S1→ S2→ S3→ S4→ S5→ S6→ S7 .

͑27͒

Suppose there are 1000 S1 molecules at time t = 0. The values of the reaction parameters c1 , . . . , c6 are 4.3, 16.6, 1.2, 2.8, 11.4, 11.9, respectively, with unit s−1. We used two methods to obtain the probability mass of S7 at time t = 1.5 s

͑1͒ by direct calculation ͑DC͒ using ͑26͒, and ͑2͒ by running 500 realizations of the SSA ͑stochastic simu-

lation algorithm͒, which is a Monte Carlo method for numerically computing the time evolution of a chemical system.9

The resulting probability mass functions are shown in

Figs. 1 and 2. From Fig. 2 it is clear that the 500 SSA

realizations are not enough for an accurate evaluation of the

probability distribution, even though the computation time it took is much longer than the direct calculation using ͑26͒. Obviously, methods based on analytic solutions are better

than those based on Monte Carlo simulations. However, in

cases of more complex reactions where neither analytic so-

lutions nor good approximations of population distributions

exist, the SSA is the method of choice. We have shown in Sec. III A that when x → ϱ and

p͑i͒͑t͒ → 0 ͓e.g., p͑i͒͑t͒ Ͻ 0.1͔, the binomial distribution ͑11͒

can be approximated by the Poisson distribution with mean = xp͑i͒͑t͒. The purpose of the approximation is to make the

analysis more tractable. In many situations, however, the condition p͑i͒͑t͒ → 0 is not valid. Using the reaction chain ͑25͒ and reaction parameter c1 , . . . , c6 of this example, the probabilities of a single molecule to be in state Si ͑i

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-7 Probabilities in systems of ﬁrst-order reactions

J. Chem. Phys. 122, 104101 ͑2005͒

FIG. 7. Variance of the population of molecule species R ͑full line͒, ARB ͑dot-dashed line͒, and B ͑dashed line͒.

= 1 , 2 , . . . , 7͒ at time t are drawn in Fig. 3, from which it can be seen that during the ﬁrst 1.5 s, many p͑i͒͑t͒’s have values

greater than 0.1, which implies that the Poisson approxima-

tion in that period is not suitable.

F. Example 6: An enzyme–substrate reaction

Let us consider an enzyme reaction where a soluble substrate with an A – B structure reacts with immobilized enzyme molecules located on the surface of cells ͓see Fig. 4͑a͔͒. The formulation of this chemical reaction is given by

c1͑2͒

c3

AB + R ARB→AR + B,

͑28͒

c2

where AB is the substrate molecule, R is the enzyme mol-

ecule, and ARB is the enzyme–substrate complex, known as

the intermediate product. This reaction is illustrated in Fig. 4͑b͒. The superscript of c1͑2͒ emphasizes that this speciﬁc probability rate constant is for a second-order reaction.

Several simpliﬁcations can be applied. First, AR can be

eliminated because the number of AR molecules has the

FIG. 8. Zoomed plots of the variance for the population of molecule species R ͑full line͒, ARB ͑dot-dashed line͒, and B ͑dashed line͒.

same probability distribution as the number of B molecules.

Such elimination obviously has no effect on the accuracy of

the obtained results. Furthermore, if the number of AB mol-

ecules is XAB, the product of XAB and the reaction parameter c͑2͒ represents the probability rate for an R molecule to transfe1r into state ARB. If c1 = c1͑2͒XAB, the reaction ͑28͒ can be simpliﬁed as

c1

c3

R ARB→B.

c2

This equation has the same structure as Eq. ͑22͒. The only difference is that the value of c1 varies with time.

A property for many enzyme reactions is that the number of substrate molecules is sufﬁciently large,16 and that the

number of substrate molecules consumed during the course

of the reaction is negligible in comparison to the total num-

ber of substrate molecules. Therefore, we can assume that

XAB is approximately equal to its initial value, and c1 can be seen as constant. Thus, after simpliﬁcation and approximation, the results in ͑23͒ can be used here.

As an example, we have calculated the mean and vari-

ance of the population of molecule species R, ARB, and B with the following initial values: XR͑0͒ = 10 000, XAB͑0͒ = 4.8176ϫ 1019, XAR͑0͒ = XARB͑0͒ = XB͑0͒ = 0, and reaction parameters c2 = 1 s−1, c3 = 100 s−1, and c1͑2͒ = 4.1930 ϫ 10−16 s−1. The volume in which the reactions are observed is equal to 0.01 l.

The obtained results are shown in Figs. 5–8. The ex-

pected population of molecule species are given in Figs. 5 and 6 ͑zoomed plot͒. The variance is shown in Figs. 7 and 8.

V. CONCLUSIONS

In this paper we have presented a general approach for studying systems of ﬁrst-order reactions. By exploiting the independence of the molecules in such systems, we can use this method for studying chemical reactions with any number of molecules and obtain the probability distributions of molecule populations.

With this approach, we have analyzed an enzyme catalyst reaction and have obtained results for both the transient and steady states of the reaction. In addition, we have also derived an analytic solution for the probability distribution of ﬁrst-order chain reactions. We have found that the Gaussian approximation is more appropriate than the Poisson approximation when the molecule population is large.

1A. M. Kierzek, J. Zaim, and P. Zielenkiewicz, J. Biol. Chem. 276, 8165 ͑2001͒. 2T. G. Liou and E. J. Campbell, J. Immunol. 157, 2624 ͑1996͒. 3B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter, Molecular Biology of the Cell ͑Garland, New York, 2002͒. 4S. Newlands, L. K. Levitt, C. S. Robinson, A. B. C. Karpf, V. R. M. Hodgson, R. P. Wade, and E. C. Hardeman, Genes Dev. 12, 2748 ͑1998͒. 5D. T. Gillespie, J. Phys. Chem. 81, 2340 ͑1977͒. 6N. G. Van Kampen, Stochastic Processes in Physics and Chemistry ͑North-Holland, Amsterdam, 1981͒. 7D. A. McQuarrie, J. Appl. Probab. 4, 413 ͑1967͒. 8I. J. Laurenzi, J. Chem. Phys. 113, 3315 ͑2000͒. 9D. T. Gillespie, J. Comput. Phys. 22, 403 ͑1976͒. 10J. J. Lukkien, J. P. L. Segers, P. A. J. Hilbers, R. J. Gelten, and A. P. J. Jansen, Phys. Rev. E 58, 2598 ͑1998͒.

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-8 Zhang et al.

11A. B. Stundzia and C. J. Lumsden, J. Comput. Phys. 127, 196 ͑1996͒. 12K. De Cock, X. Zhang, M. F. Bugallo, and P. M. Djurić, Proceedings of

the 12th European Signal Processing Conference ͑2004͒. 13G. Grimmett and D. Stirzaker, Probability and Random Processes ͑Oxford

University Press, Oxford, 2001͒. 14M. Rathinam, L. R. Petzold, Y. Cao, and D. T. Gillespie, J. Chem. Phys.

J. Chem. Phys. 122, 104101 ͑2005͒

119, 12784 ͑2003͒. 15K. De Cock, X. Zhang, M. F. Bugallo, and P. M. Djurić, J. Chem. Phys.

121, 3347 ͑2004͒. 16X. Zhang, K. De Cock, M. F. Bugallo, and P. M. Djurić, Proceedings of

the 2003 IEEE Workshop on Statistical Signal Processing ͑2003͒.

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

A general method for the computation of probabilities in systems of ﬁrst order chemical reactions

Xueying Zhang,a͒ Katrien De Cock,b͒ Mónica F. Bugallo,c͒ and Petar M. Djurićd͒ COSINE Laboratory, Department of Electrical & Computer Engineering, Stony Brook University, Stony Brook, New York 11794-2350

͑Received 21 June 2004; accepted 9 December 2004; published online 7 March 2005͒

We present a general method for the computation of molecular population distributions in a system of ﬁrst-order chemical reactions. The method is based on the observation that the molecules in ﬁrst-order reactions do not interact with each other. Therefore, the population distributions are a convolution of densities for one molecule. With this method one can study reactions involving any number of molecules. Such analysis is demonstrated on several examples, including an enzyme catalyst reaction and a ﬁrst-order reaction chain. © 2005 American Institute of Physics. ͓DOI: 10.1063/1.1855311͔

I. INTRODUCTION

In small biological systems, it is often more appropriate to model the chemical reactions in a stochastic way rather than with the traditional differential equations for evolution of concentrations. The reasons are ͑a͒ the number of molecules involved in a biological system can be very small;1 ͑b͒ there are many situations in which molecules are not homogeneously distributed;2,3 ͑c͒ the change of molecular population levels is a discrete integer amount; and ͑d͒ the population change is not a deterministic process.4,5

There are mainly two approaches for the stochastic study of the number of molecules in biochemical reactions: the ﬁrst is based on the analysis of the master equation6–8 and the second relies on Monte Carlo simulation methods.9–11 Here, we adopt the ﬁrst approach.

In 1967, McQuarrie summarized the stochastic study using the master equation approach.7 For a reaction, it is assumed that in an inﬁnitesimal time interval, the probability of having one reaction per unit reactant molecule combination is proportional to the length of the time interval. A probability difference function is ﬁrst obtained based on the assumption, which leads to a differential-difference equation called the master equation of the reaction. The momentgenerating function is then employed to transform the master equation into a partial differential equation. It can sometimes be solved to obtain the analytic solution. More often, the mean and variance of the number of molecules can only be obtained.

In 2000, Laurenzi introduced a different way of solving the master equation.8 Instead of using the momentgenerating function, he applied the Laplace transform. In this way, solving the partial differential equations is avoided, which is important for more complicated reactions. Instead, one needs to solve a set of linear equations. We will show

a͒Electronic mail: [email protected] b͒Electronic mail: [email protected] c͒Electronic mail: [email protected] d͒Electronic mail: [email protected]

how to obtain the solution of the master equation for any system of ﬁrst-order reactions in yet another way.

In this paper we study systems of ﬁrst-order reactions, i.e., combinations of ͑coupled͒ reactions of the type Si → Sj. Because in such systems molecules do not chemically interact, the analysis is much simpler than when second-order reactions are involved. Although this situation might occur only rarely in practice, our analysis is also useful for systems where the ﬁrst-order reactions are dominant and for systems where second-order reactions can be approximated by ﬁrstorder reactions ͑see Sec. IV for an example͒.

The independence of the molecules in a ﬁrst-order reaction system is exploited to derive the population distributions of the molecules in the system. Instead of solving the master equation for the complete system, where the number of states grows with the number of molecules, we ﬁrst solve the master equation for only one molecule and use its solution to ﬁnd the population distributions for all species.

The structure of the paper is as follows. In Sec. II we solve the master equation for one molecule in a ﬁrst-order reaction system. The resulting probability distribution functions are used in Sec. III to obtain the population distributions when more molecules are introduced. In Sec. IV we illustrate the general results of Sec. III with speciﬁc examples. Finally, in Sec. V we give some concluding remarks.

II. ONE MOLECULE

In this section we study the master equation for one molecule in a system of ﬁrst-order reactions. The solution will then be used to solve more general cases in Sec. III.

Assume that there are M molecule species, S1 , S2 , . . . , SM in the reaction system, and let the probability rate constant of the reaction Si → Sj be denoted by cij. For a ﬁrst-order reaction the speciﬁc probability rate constant is such that

cij⌬t͑i j͒ ϵ the probability for the first-order reaction

0021-9606/2005/122͑10͒/104101/8/$22.50

122, 104101-1

© 2005 American Institute of Physics

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-2 Zhang et al.

J. Chem. Phys. 122, 104101 ͑2005͒

cij

Si→Sj to happen per molecule in an infinitesimally small time interval ⌬t.

The master equation describes the evolution of the molecules’ population distribution in the chemical system as a function of time. When there is only one molecule in a ﬁrstorder reaction system, the master equation gives the time evolution of the probability that this molecule has become a certain chemical species. Let p͑i͒͑t͒ be the probability that the molecule is an Si molecule at time t. For example, in the simple system of reactions

c12 c23

S1 S2→S3,

͑1͒

c21

the differential equation describing the time evolution of p͑2͒ is the following ﬁrst-order linear differential equation:

dpd͑2t͒͑t͒ = c12 p͑1͒͑t͒ − ͑c21 + c23͒p͑2͒͑t͒. ͑2͒

We assume in the rest of this paper that the initial values for the probabilities are equal to zero: p͑i͒͑0͒ = 0, i = 1 , . . . , M. The injection of the molecule in the system is modeled by a source probability density function. For the given example, if the molecule begins as an S2 molecule, the differential equation becomes

dpd͑2t͒͑t͒ = c12 p͑1͒͑t͒ − ͑c21 + c23͒p͑2͒͑t͒ + f͑2͒͑t͒, ͑3͒

where f ͑2͒͑t͒ is the source’s probability density, i.e., f ͑2͒͑t͒⌬t is the probability that a molecule from the source is injected in state S2 in the interval ͓t , t + ⌬t͔. If the injection time is known precisely, e.g., t = , the source’s probability density function is a delta function: f ͑2͒͑t͒ = ␦͑t − ͒. Otherwise, the density function for the injection time of the molecule can have forms such as ae−a͑t−͒u͑t − ͒, a␦͑t − 1͒ + ͑1 − a͒␦͑t − 2͒, and

͚ϱ ne−

␦͑t − n͒. n=0 n!

In general, the system can be represented by a Markov chain with a transition matrix

͚ c11 c12 ¯ c1M

c21 c22 ¯ c2M

M

C=

, where cii = −

ci j ,

]]

]

j=1,j i

cM1 cM2 ¯ cMM

͑4͒

and the set of ﬁrst-order linear differential equations has the following form:

dp͑1͒͑t͒

dt

p͑1͒͑t͒

f ͑1͒͑t͒

dp͑2͒͑t͒

p͑2͒͑t͒

f ͑2͒͑t͒

dt

= CT

+

,

͑5͒

]

]

] dp͑M͒͑t͒

p͑M͒͑t͒

f ͑M͒͑t͒

dt

where

M

͚ p͑i͒͑t͒ = 1,

i=1

and

Ά͵ · ϱ f ͑i͒͑t͒dt = 1 t=0

if the molecule from source is injected

into state i, f ͑i͒͑t͒ = 0 otherwise.

The equations can be solved by ﬁrst applying the Laplace transform, followed by solving the resulting algebraic equation, and ﬁnally using the inverse Laplace transform. More speciﬁcally, after taking the Laplace transform on both sides of ͑5͒, we obtain

sp = CTp + f,

͑6͒

where p is an M ϫ 1 vector whose elements p͑i͒͑s͒ are the Laplace transforms of p͑i͒͑t͒, the sought-after probability

density functions, and f is an M ϫ 1 vector with elements f ͑i͒͑s͒, which are the Laplace transforms of the source functions f ͑i͒͑t͒. Note that f has only one nonzero element.

After solving for p in ͑6͒, we get

p = ͑sI − CT͒−1f.

͑7͒

Let

L = ͑sI − CT͒−1

͑8͒

and

G = L−1͑L͒,

͑9͒

where L−1 is the inverse Laplace operator. If the elements of G are denoted by gij͑t͒, and f ͑m͒͑s͒ is the nonzero element of f, the solutions p͑i͒͑t͒ are given by

p͑i͒͑t͒ = gim͑t͒ f͑m͒͑t͒, i = 1,2, . . . , M ,

͑10͒

where the symbol denotes convolution. Finding the inverse Laplace transforms of the elements of L may be tedious but it is straightforward.

This method is illustrated in details by example 1. As already mentioned, the technique of using the Laplace transform for solving the master equation was introduced by Laurenzi in Ref. 8. The difference between his and our approach is that Laurenzi was solving the many-molecule master equation, while here we only solve the one-molecule master equation.

In the deterministic framework, the concentrations of the species are found by solving the same differential equations

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-3 Probabilities in systems of ﬁrst-order reactions

͑5͒, without the stochastic sources, but taking into account the initial concentrations. Furthermore, the so-obtained results there yield only the means of the stochastic distributions.

III. MORE MOLECULES

In this section we treat the case when several molecules are injected into the system. First, in Sec. III A, we assume that all molecules are inserted by the same source, so that they all start as the same chemical species. In Sec. III B, the molecules can originate from different sources.

A. One source

Suppose that there are x molecules of species Sk injected into the system with the same probability density f ͑k͒͑t͒, and no molecules of other species. Denote the number of molecules at state Si ͑i = 1 , . . . , M͒ at time t by N͑i͒͑t͒, and let K͑ji͒͑t͒ = 1 if molecule j is at state Si at time t and K͑ji͒͑t͒ = 0 if it is not. Thus, N͑i͒͑t͒ = ⌺xj=1K͑ji͒͑t͒. Note that K1͑i͒͑t͒, K2͑i͒͑t͒ , . . . , Kx͑i͒͑t͒ are independent Bernoulli random variables with the same probability p͑i͒͑t͒, which can be computed by solving ͑5͒. Therefore, N͑i͒͑t͒ has a binomial distribution. The probability that there are n molecules in state Si at time t is

ͩ ͪ Pn͑i͒͑t͒ = nx ͑p͑i͒͑t͒͒n͑1 − p͑i͒͑t͒͒x−n, ͑11͒

with mean i͑t͒ and variance i2͑t͒ equal to i͑t͒ = xp͑i͒͑t͒,

i2͑t͒ = xp͑i͒͑t͒͑1 − p͑i͒͑t͒͒,

respectively. Moreover, the joint probability of ﬁnding ni molecules in state Si ͑i = 1 , . . . , M͒ at time t is equal to

P͑n , . . . ,n ;t͒ = x! ͑p͑1͒͑t͒͒n1 ͑p͑M͒͑t͒͒nM,

1

M n1! ¯ nM!

¯

͑12͒

where ⌺iM=1ni = x. Equation ͑11͒ is one of the marginals of ͑12͒.

When the number of molecules is large, i.e., when xp͑i͒͑t͒͑1 − p͑i͒͑t͒͒ ӷ 1, the binomial distribution ͑11͒ can be

well approximated by a Gaussian distribution, or

ͱ Pn͑i͒͑t͒ Ӎ

1 e−͓͑n − i͑t͒͒2/2i2͑t͔͒ . 2i2͑t͒

The purpose for employing the Gaussian approximation is to save computation time when the number of molecules is very large, as well as to make the analysis more tractable. Numerical examples of the comparison between analytic distributions and their Gaussian approximations can be found in Ref. 12.

The same purpose can be achieved if the binomial distribution ͑11͒ is approximated by a Poisson distribution. Conditions for the validity of such approximation are stricter. First, the number of molecules should be large. Second, the probability for a single molecule in state i should be small

J. Chem. Phys. 122, 104101 ͑2005͒

͓e.g., p͑i͒͑t͒ Ͻ 0.1͔. Finally, the product of x and p͑i͒͑t͒ should

be a sizable number. Under such conditions, the variance of the binomial distribution roughly equals the mean. Then ͑11͒ can be approximated by the Poisson distribution with parameter i = xp͑i͒͑t͒, or

P͑i͒͑t͒ Ӎ in e−i.

n

n!

B. More sources

If there are several sources which are independent from

each other, the probability mass function for state Si is the convolution of binomial distributions. Suppose there are J

independent sources. Denote the number of molecules at state Si ͑i = 1 , . . . , M͒ at time t by N͑i͒͑t͒, and denote the number of molecules at state Si that are injected from source j by K͑ji͒͑t͒ ͑j = 1 , . . . , J͒. Then, N͑i͒͑t͒ = ⌺Jj=1K͑ji͒͑t͒. Since K͑ji͒͑t͒ ͑j = 1 , . . . , J͒ are independent random variables, the probability mass function of the random variable N͑i͒͑t͒ is the con-

volution of the probability mass functions of the random variables K͑ji͒͑t͒.13

For example, let x Sk molecules be injected by the source functions f ͑k͒͑t͒ and y Sl molecules by f ͑l͒͑t͒. Each of the x Sk molecules gives rise to the state probabilities pk͑i͒͑t͒ ͑i = 1 , . . . , M͒, whereas a molecule introduced by f ͑l͒͑t͒, to p͑i͒͑t͒ ͑i = 1 , . . . , M͒. Then, the number of molecules at state Sil that are injected from source k is Kk͑i͒͑t͒, which has a binomial distribution ͑ nx ͒͑pk͑i͒͑t͒͒n͑1 − pk͑i͒͑t͒͒x−n. Similarly, the number of molecules at state Si that are injected from source l is Kl͑i͒͑t͒, which has a binomial distribution ͑ ny ͒͑pl͑i͒͑t͒͒n͑1 − pl͑i͒͑t͒͒y−n. The probability mass function of N͑i͒͑t͒, the total number of molecules at state Si is equal to the convolution of two binomial distributions, i.e.,

ͩ ͪ Pn͑i͒͑t͒ = nx ͑pk͑i͒͑t͒͒n͑1 − pk͑i͒͑t͒͒x−n

ͩ ͪ ny ͑pl͑i͒͑t͒͒n͑1 − pl͑i͒͑t͒͒y−n, ͑13͒

where the convolution is along the n axis, that is

ϱ

͚ y1͑n͒ y2͑n͒ = y1͑m͒y2͑n − m͒. m=0

The mean i͑t͒ and variance i2͑t͒ of N͑i͒͑t͒ are

i͑t͒ = xpk͑i͒͑t͒ + ypl͑i͒͑t͒,

͑14͒

i2͑t͒ = xpk͑i͒͑t͒͑1 − pk͑i͒͑t͒͒ + ypl͑i͒͑t͒͑1 − pl͑i͒͑t͒͒,

͑15͒

respectively. As we mentioned before, when the number of molecules

is large, the binomial distribution can be well approximated by a Gaussian distribution with the same mean and variance. One property of Gaussian distributions is that the convolution of two Gaussian distributions yields a Gaussian distribution. Therefore, the convolution result in ͑13͒ can be approximated by a Gaussian distribution with mean ͑14͒ and variance ͑15͒. The Poisson distribution has the same prop-

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-4 Zhang et al.

erty: the convolution of Poisson distributions with parameters ␣ and , respectively, yields a Poisson distribution with parameter = ␣ + .

IV. EXAMPLES

In this section several examples of systems of ﬁrst-order reactions are discussed. For some systems, an analytic solution for the population distribution is derived. This section also includes two numerical examples: the ﬁrst is a onedirectional reaction chain and the second an enzyme substrate reaction.

A. Example 1: The ﬁrst-order irreversible reaction

c

Consider the ﬁrst-order reaction S1→S2, where only

molecules S1 are injected into the system with probability density function f ͑1͒͑t͒ = ␦͑t͒. The ﬁrst step is to rewrite ͑5͒ with M = 2, c12= c, c21= 0, and f ͑2͒͑t͒ = 0

ͩ ͪ ͩ ͪ ͩ ͪ dp͑1͒͑t͒

dt

− c c T p͑1͒͑t͒ ␦͑t͒

dp͑2͒͑t͒ = 0 0 p͑2͒͑t͒ + 0 .

͑16͒

dt

Taking the Laplace transform leads to

ͩ ͪ ͩ ͪ ͩ ͪ ͩ ͪ p͑1͒͑s͒

− c c T p͑1͒͑s͒

1

s p͑2͒͑s͒ = 0 0 p͑2͒͑s͒ + 0 .

The matrix L from ͑8͒ is given by

1

0

s+c

L=

.

1 11

−

s s+c s

After taking the inverse Laplace transform of L, we obtain

ͩ ͪ e−ct

0

G = 1 − e−ct u͑t͒ ,

where u͑t͒ is the unit step function. The ﬁnal solution is

p͑1͒͑t͒ = e−ctu͑t͒,

͑17͒

p͑2͒͑t͒ = ͑1 − e−ct͒u͑t͒.

͑18͒

In the rest of the paper we will drop the use of u͑t͒, and so

we point out that the given solutions for p͑t͒ are valid only for t ജ 0.When f ͑1͒͑t͒ is some other density function

͵ p͑1͒͑t͒ =

ϱ

f ͑1͒͑͒e−c͑t−͒d,

0

͵ p͑2͒͑t͒ =

ϱ

f ͑1͒͑͒͑1 − e−c͑t−͒͒d.

0

Note that if an S2 instead of an S1 molecule was injected into the system with f ͑2͒͑t͒, p͑1͒͑t͒ = 0 and p͑2͒͑t͒ = ͐0t f ͑2͒͑͒d.

By substituting the probabilities ͑17͒ and ͑18͒ in ͑11͒, we obtain the well-known result for f ͑͑t1͒͒ = ␦͑t͒

J. Chem. Phys. 122, 104101 ͑2005͒

ͩ ͪ Pn͑1͒͑t͒ = nx e−nct͑1 − e−ct͒x−n, ͩ ͪ Pn͑2͒͑t͒ = nx ͑1 − e−ct͒ne−c͑x−n͒t.

B. Example 2. The ﬁrst-order reversible reaction

c1

Consider the ﬁrst-order reversible reaction S1 S2 which

c2

starts with x S1 molecules and y S2 molecules. Each of the x

S1 molecules gives rise to the following probabilities, which are obtained by solving ͑5͒ with M = 2, f ͑1͒͑t͒ = ␦͑t͒, and f ͑2͒͑t͒ = 0:

p͑1͒͑t͒ = 1 ͑c + c e−͑c1+c2͒t͒,

1

c1 + c2 2 1

p͑2͒͑t͒ = c1 ͑1 − e−͑c1+c2͒t͒.

1

c1 + c2

͑19͒

For the effect of each of the y S2 molecules, we take f ͑1͒͑t͒ = 0 and f͑2͒͑t͒ = ␦͑t͒ and obtain

p͑1͒͑t͒ = c2 ͑1 − e−͑c2+c1͒t͒,

2

c2 + c1

p͑2͒͑t͒ = 1 ͑c + c e−͑c2+c1͒t͒.

2

c2 + c1 1 2

͑20͒

By substituting the probabilities of ͑19͒ and ͑20͒ in ͑13͒, we

obtain the distribution for the S1 molecules.

ͩ ͪͩ ͪ Pn͑1͒͑t͒ = nx

1

n

͑c + c e−͑c1+c2͒t͒

c1 + c2 2 1

ͩ ͪ ϫ 1 −

1

x−n

͑c + c e−͑c1+c2͒t͒

c1 + c2 2 1

ͩ ͪͩ ͪ y

n

c2

n

͑1 − e−͑c1+c2͒t͒

c2 + c1

ͩ ͪ ϫ 1 −

c2

y−n

͑1 − e−͑c1+c2͒t͒ ,

c2 + c1

͑21͒

and similarly for the S2 molecules. Note that when only S1 or only S2 molecules are present

at t = 0, the distribution reduces to the ﬁrst and second convolution factor in ͑21͒, respectively. The population distribu-

tion for this special case has already been given by McQuarrie.7

Also note that deriving the mean and variance for the

number of S1 and S2 molecules is straightforward, as explained in Sec. III B, and leads to the same results as re-

ported in Refs. 14 and 15.

C. Example 3: A combination of the ﬁrst-order irreversible and ﬁrst-order reversible reaction

Consider the reaction system

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-5 Probabilities in systems of ﬁrst-order reactions

J. Chem. Phys. 122, 104101 ͑2005͒

FIG. 1. Direct calculation of the probability mass function of S7 at time t = 1.5 s.

c1 c3

S1 S2→S3,

c2

͑22͒

and assume there are x S1 molecules and y S2 molecules at t = 0. Again, we ﬁrst solve ͑5͒ with f ͑1͒͑t͒ = ␦͑t͒, f ͑2͒͑t͒ = f ͑3͒͑t͒ = 0 to ﬁnd the probabilities due to one S1 molecule

and obtain

ͱ ͱ p1͑1͒͑t͒ =

m − c1 + c2 + c3 e−͑1/2͒͑−ͱm+c1+c2+c3͒t 2m

ͱ + m + c1 − c2 − c3 e−͑1/2͒͑ͱm+c1+c2+c3͒t, ͱ2 m

p1͑2͒͑t͒

=

ͱc1 e−͑1/2͒͑−ͱm+c1+c2+c3͒t

m

−

ͱc1 e−͑1/2͒͑ͱm+c1+c2+c3͒t,

m

͑23͒

FIG. 3. The probabilities for a single molecule to be in state Si ͑i = 1 , 2 , . . . , 7͒ at time t.

ͱ ͱ p1͑3͒͑t͒ = 1 −

m + c1 + c2 + c3 e−͑1/2͒͑−ͱm+c1+c2+c3͒t 2m

ͱ − m − c1 − c2 − c3 e−͑1/2͒͑ͱm+c1+c2+c3͒t, ͱ2 m

where

m = c12 + 2c1c2 − 2c1c3 + c22 + 2c2c3 + c32.

Next, we solve ͑5͒ with f ͑1͒͑t͒ = f ͑3͒͑t͒ = 0, f ͑2͒ = ␦͑t͒ to ﬁnd the probabilities due to one S2 molecule and ﬁnd that

p2͑1͒͑t͒

=

ͱc2 e−͑1/2͒͑−ͱm+c1+c2+c3͒t

m

−

ͱc2 e−͑1/2͒͑ͱm+c1+c2+c3͒t,

m

ͱ ͱ p2͑2͒͑t͒ =

m + c1 − c2 − c3 e−͑1/2͒͑−ͱm+c1+c2+c3͒t 2m

ͱ + m − c1 + c2 + c3 e−͑1/2͒͑ͱm+c1+c2+c3͒t, ͱ2 m

͑24͒

ͱ ͱ p2͑3͒͑t͒ = 1 −

m + c1 + c2 − c3 e−͑1/2͒͑−ͱm+c1+c2−c3͒t 2m

ͱ − m − c1 − c2 + c3 e−͑1/2͒͑ͱm+c1+c2+c3͒t. ͱ2 m

The substitution of the probabilities in ͑13͒ by the probabilities in ͑23͒ and ͑24͒ yields the population distribution for the S1, S2, and S3 molecules.

FIG. 2. The probability mass of S7 at time t = 1.5 s calculated using the SSA, with 500 realizations.

D. Example 4: A ﬁrst-order reaction chain

Consider a ﬁrst-order reaction chain in its general form

c1 c2

ci−1 ci

cN−1

S1→S2→ ¯ → Si→ ¯ → SN.

͑25͒

Assume that ci cj when i j ͑i = 1 , . . . , N and j = 1 , . . . , N͒ and that x molecules enter the system in state S1 at time t = 0. The input probability distribution of state S1 is f ͑1͒͑t͒

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-6 Zhang et al.

J. Chem. Phys. 122, 104101 ͑2005͒

FIG. 4. ͑a͒ A cell on which enzyme molecules are immobilized. The substrate molecules ﬂoat around freely. ͑b͒ Illustration of the chemical reaction. The substrate molecule AB approaches the enzyme molecule R ͑left͒. The enzyme recognizes the A molecule of the AB structure ͑middle͒. The enzyme destroys the bond between A and B, A remains attached to R, and B is released into the solution ͑right͒.

= ␦͑t͒. Then, using the method introduced in the previous sections, the state probabilities for one molecule are equal to

Ά ͚͚ · p͑i͒͑t͒ =

i

1 M͑1,i͒c e−cjt, ci j=1 j j

ci Ͼ 0

i−1

1 − M͑j1,i−1͒e−cjt,

j=1

ci = 0,i ജ 2,

1, ci = 0,i = 1

where

Ά ͟ M͑pk,m͒ =

m

cj ,

j=k,j p cj − cp 1, m = k.

mϾk·

The probability that n of the x molecules are in state Si at time t is given by

FIG. 5. Mean for the population of molecule species R ͑full line͒, ARB ͑dot-dashed line͒, and B ͑dashed line͒.

FIG. 6. Zoomed plots of the mean for the population of molecule species R ͑full line͒, ARB ͑dot-dashed line͒, and B ͑dashed line͒.

ͩ ͪ Pn͑i͒͑t͒ = nx ͑p͑i͒͑t͒͒n͑1 − p͑i͒͑t͒͒x−n. ͑26͒

If there are ci’s equal to cj’s, a closed-form solution is also possible.

E. Example 5: A numerical example for the ﬁrst-order reaction chain

Consider the ﬁrst-order reaction chain

c1 c2 c3 c4 c5 c6

S1→ S2→ S3→ S4→ S5→ S6→ S7 .

͑27͒

Suppose there are 1000 S1 molecules at time t = 0. The values of the reaction parameters c1 , . . . , c6 are 4.3, 16.6, 1.2, 2.8, 11.4, 11.9, respectively, with unit s−1. We used two methods to obtain the probability mass of S7 at time t = 1.5 s

͑1͒ by direct calculation ͑DC͒ using ͑26͒, and ͑2͒ by running 500 realizations of the SSA ͑stochastic simu-

lation algorithm͒, which is a Monte Carlo method for numerically computing the time evolution of a chemical system.9

The resulting probability mass functions are shown in

Figs. 1 and 2. From Fig. 2 it is clear that the 500 SSA

realizations are not enough for an accurate evaluation of the

probability distribution, even though the computation time it took is much longer than the direct calculation using ͑26͒. Obviously, methods based on analytic solutions are better

than those based on Monte Carlo simulations. However, in

cases of more complex reactions where neither analytic so-

lutions nor good approximations of population distributions

exist, the SSA is the method of choice. We have shown in Sec. III A that when x → ϱ and

p͑i͒͑t͒ → 0 ͓e.g., p͑i͒͑t͒ Ͻ 0.1͔, the binomial distribution ͑11͒

can be approximated by the Poisson distribution with mean = xp͑i͒͑t͒. The purpose of the approximation is to make the

analysis more tractable. In many situations, however, the condition p͑i͒͑t͒ → 0 is not valid. Using the reaction chain ͑25͒ and reaction parameter c1 , . . . , c6 of this example, the probabilities of a single molecule to be in state Si ͑i

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-7 Probabilities in systems of ﬁrst-order reactions

J. Chem. Phys. 122, 104101 ͑2005͒

FIG. 7. Variance of the population of molecule species R ͑full line͒, ARB ͑dot-dashed line͒, and B ͑dashed line͒.

= 1 , 2 , . . . , 7͒ at time t are drawn in Fig. 3, from which it can be seen that during the ﬁrst 1.5 s, many p͑i͒͑t͒’s have values

greater than 0.1, which implies that the Poisson approxima-

tion in that period is not suitable.

F. Example 6: An enzyme–substrate reaction

Let us consider an enzyme reaction where a soluble substrate with an A – B structure reacts with immobilized enzyme molecules located on the surface of cells ͓see Fig. 4͑a͔͒. The formulation of this chemical reaction is given by

c1͑2͒

c3

AB + R ARB→AR + B,

͑28͒

c2

where AB is the substrate molecule, R is the enzyme mol-

ecule, and ARB is the enzyme–substrate complex, known as

the intermediate product. This reaction is illustrated in Fig. 4͑b͒. The superscript of c1͑2͒ emphasizes that this speciﬁc probability rate constant is for a second-order reaction.

Several simpliﬁcations can be applied. First, AR can be

eliminated because the number of AR molecules has the

FIG. 8. Zoomed plots of the variance for the population of molecule species R ͑full line͒, ARB ͑dot-dashed line͒, and B ͑dashed line͒.

same probability distribution as the number of B molecules.

Such elimination obviously has no effect on the accuracy of

the obtained results. Furthermore, if the number of AB mol-

ecules is XAB, the product of XAB and the reaction parameter c͑2͒ represents the probability rate for an R molecule to transfe1r into state ARB. If c1 = c1͑2͒XAB, the reaction ͑28͒ can be simpliﬁed as

c1

c3

R ARB→B.

c2

This equation has the same structure as Eq. ͑22͒. The only difference is that the value of c1 varies with time.

A property for many enzyme reactions is that the number of substrate molecules is sufﬁciently large,16 and that the

number of substrate molecules consumed during the course

of the reaction is negligible in comparison to the total num-

ber of substrate molecules. Therefore, we can assume that

XAB is approximately equal to its initial value, and c1 can be seen as constant. Thus, after simpliﬁcation and approximation, the results in ͑23͒ can be used here.

As an example, we have calculated the mean and vari-

ance of the population of molecule species R, ARB, and B with the following initial values: XR͑0͒ = 10 000, XAB͑0͒ = 4.8176ϫ 1019, XAR͑0͒ = XARB͑0͒ = XB͑0͒ = 0, and reaction parameters c2 = 1 s−1, c3 = 100 s−1, and c1͑2͒ = 4.1930 ϫ 10−16 s−1. The volume in which the reactions are observed is equal to 0.01 l.

The obtained results are shown in Figs. 5–8. The ex-

pected population of molecule species are given in Figs. 5 and 6 ͑zoomed plot͒. The variance is shown in Figs. 7 and 8.

V. CONCLUSIONS

In this paper we have presented a general approach for studying systems of ﬁrst-order reactions. By exploiting the independence of the molecules in such systems, we can use this method for studying chemical reactions with any number of molecules and obtain the probability distributions of molecule populations.

With this approach, we have analyzed an enzyme catalyst reaction and have obtained results for both the transient and steady states of the reaction. In addition, we have also derived an analytic solution for the probability distribution of ﬁrst-order chain reactions. We have found that the Gaussian approximation is more appropriate than the Poisson approximation when the molecule population is large.

1A. M. Kierzek, J. Zaim, and P. Zielenkiewicz, J. Biol. Chem. 276, 8165 ͑2001͒. 2T. G. Liou and E. J. Campbell, J. Immunol. 157, 2624 ͑1996͒. 3B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter, Molecular Biology of the Cell ͑Garland, New York, 2002͒. 4S. Newlands, L. K. Levitt, C. S. Robinson, A. B. C. Karpf, V. R. M. Hodgson, R. P. Wade, and E. C. Hardeman, Genes Dev. 12, 2748 ͑1998͒. 5D. T. Gillespie, J. Phys. Chem. 81, 2340 ͑1977͒. 6N. G. Van Kampen, Stochastic Processes in Physics and Chemistry ͑North-Holland, Amsterdam, 1981͒. 7D. A. McQuarrie, J. Appl. Probab. 4, 413 ͑1967͒. 8I. J. Laurenzi, J. Chem. Phys. 113, 3315 ͑2000͒. 9D. T. Gillespie, J. Comput. Phys. 22, 403 ͑1976͒. 10J. J. Lukkien, J. P. L. Segers, P. A. J. Hilbers, R. J. Gelten, and A. P. J. Jansen, Phys. Rev. E 58, 2598 ͑1998͒.

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

104101-8 Zhang et al.

11A. B. Stundzia and C. J. Lumsden, J. Comput. Phys. 127, 196 ͑1996͒. 12K. De Cock, X. Zhang, M. F. Bugallo, and P. M. Djurić, Proceedings of

the 12th European Signal Processing Conference ͑2004͒. 13G. Grimmett and D. Stirzaker, Probability and Random Processes ͑Oxford

University Press, Oxford, 2001͒. 14M. Rathinam, L. R. Petzold, Y. Cao, and D. T. Gillespie, J. Chem. Phys.

J. Chem. Phys. 122, 104101 ͑2005͒

119, 12784 ͑2003͒. 15K. De Cock, X. Zhang, M. F. Bugallo, and P. M. Djurić, J. Chem. Phys.

121, 3347 ͑2004͒. 16X. Zhang, K. De Cock, M. F. Bugallo, and P. M. Djurić, Proceedings of

the 2003 IEEE Workshop on Statistical Signal Processing ͑2003͒.

Downloaded 30 Apr 2010 to 129.49.69.111. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp