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 first 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 first-order chemical reactions. The method is based on the observation that the molecules in first-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 first-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 first 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 first approach.
In 1967, McQuarrie summarized the stochastic study using the master equation approach.7 For a reaction, it is assumed that in an infinitesimal 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 first 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 first-order reactions in yet another way.
In this paper we study systems of first-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 first-order reactions are dominant and for systems where second-order reactions can be approximated by firstorder reactions ͑see Sec. IV for an example͒.
The independence of the molecules in a first-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 first solve the master equation for only one molecule and use its solution to find 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 first-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 specific 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 first-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 first-order reaction the specific 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 firstorder 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 first-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 first-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 first applying the Laplace transform, followed by solving the resulting algebraic equation, and finally using the inverse Laplace transform. More specifically, 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 first-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 finding 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 first-order reactions are discussed. For some systems, an analytic solution for the population distribution is derived. This section also includes two numerical examples: the first is a onedirectional reaction chain and the second an enzyme substrate reaction.
A. Example 1: The first-order irreversible reaction
c
Consider the first-order reaction S1→S2, where only
molecules S1 are injected into the system with probability density function f ͑1͒͑t͒ = ␦͑t͒. The first 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 final 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 first-order reversible reaction
c1
Consider the first-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 first 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 first-order irreversible and first-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 first-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 first solve ͑5͒ with f ͑1͒͑t͒ = ␦͑t͒, f ͑2͒͑t͒ = f ͑3͒͑t͒ = 0 to find 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 find the probabilities due to one S2 molecule and find 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 first-order reaction chain
Consider a first-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 float 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 first-order reaction chain
Consider the first-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 first-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 first 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 specific probability rate constant is for a second-order reaction.
Several simplifications 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 simplified 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 sufficiently 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 simplification 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 first-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 first-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 first 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 first-order chemical reactions. The method is based on the observation that the molecules in first-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 first-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 first 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 first approach.
In 1967, McQuarrie summarized the stochastic study using the master equation approach.7 For a reaction, it is assumed that in an infinitesimal 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 first 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 first-order reactions in yet another way.
In this paper we study systems of first-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 first-order reactions are dominant and for systems where second-order reactions can be approximated by firstorder reactions ͑see Sec. IV for an example͒.
The independence of the molecules in a first-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 first solve the master equation for only one molecule and use its solution to find 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 first-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 specific 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 first-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 first-order reaction the specific 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 firstorder 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 first-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 first-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 first applying the Laplace transform, followed by solving the resulting algebraic equation, and finally using the inverse Laplace transform. More specifically, 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 first-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 finding 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 first-order reactions are discussed. For some systems, an analytic solution for the population distribution is derived. This section also includes two numerical examples: the first is a onedirectional reaction chain and the second an enzyme substrate reaction.
A. Example 1: The first-order irreversible reaction
c
Consider the first-order reaction S1→S2, where only
molecules S1 are injected into the system with probability density function f ͑1͒͑t͒ = ␦͑t͒. The first 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 final 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 first-order reversible reaction
c1
Consider the first-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 first 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 first-order irreversible and first-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 first-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 first solve ͑5͒ with f ͑1͒͑t͒ = ␦͑t͒, f ͑2͒͑t͒ = f ͑3͒͑t͒ = 0 to find 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 find the probabilities due to one S2 molecule and find 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 first-order reaction chain
Consider a first-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 float 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 first-order reaction chain
Consider the first-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 first-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 first 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 specific probability rate constant is for a second-order reaction.
Several simplifications 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 simplified 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 sufficiently 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 simplification 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 first-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 first-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