PHYSICAL REVIEW E 81, 051920 共2010兲

Discreteness-induced slow relaxation in reversible catalytic reaction networks
Akinori Awazu1 and Kunihiko Kaneko2,3

1

Department of Mathematical and Life Sciences, Hiroshima University, Kagami-yama 1-3-1, Higashi-Hiroshima 739-8526, Japan2
Department of Basic Science, University of Tokyo, Komaba, Meguro-ku, Tokyo 153-8902, Japan
3
ERATO Complex Systems Biology, JST, Komaba, Meguro-ku, Tokyo 153-8902, Japan
共Received 1 March 2010; revised manuscript received 15 April 2010; published 21 May 2010兲
Slowing down of the relaxation of the fluctuations around equilibrium is investigated both by stochastic
simulations and by analysis of master equation of reversible reaction networks consisting of reactions between
a pair of resource and the corresponding high-energy product that works as a catalyst for another resourceproduct reaction. As the number of molecules N is decreased, the relaxation time to equilibrium is prolonged
due to the deficiency of catalysts, as demonstrated by the amplification compared to that by the continuum
limit. This amplification ratio of the relaxation time is represented by a scaling function as h = N exp共−␤V兲, andit becomes prominent as N becomes less than a critical value h ⬃ 1, where ␤ is the inverse temperature and V
is the energy required to the transformation from resources to the corresponding products.
DOI: 10.1103/PhysRevE.81.051920

PACS number共s兲: 87.16.Yc, 05.40.⫺a, 82.39.Rt

I. INTRODUCTION

The study of reaction processes in catalytic reaction networks is generally important to understand the dynamics andfluctuations in biochemical systems and their functionality.
Obviously, understanding the generic features of equilibrium
characteristics and relaxation to equilibrium is the first step
toward gaining such an understanding. Indeed, such reaction
systems often exhibit anomalous slow relaxation to equilibrium due to some kinetic constraints such as diffusioninfluenced 共limited兲 reaction 关1兴 and formations of transient
Turing patterns 关2兴. In this paper, we consider a mechanism
to realize such a slow relaxation in catalytic reaction networks, where the discreteness in molecule number that may
reach zero induces drastic slowing down.
Most intracellular reactions progress with the aid of catalysts 共proteins兲, whereas catalysts have to be synthesized as aresult of such catalytic reactions. Indeed, reaction dynamics
in catalytic networks have been extensively investigated. In
most such studies, a limiting case with a strong nonequilibrium condition was assumed by adopting a unidirectional
reaction process 共i.e., by neglecting backward reactions兲. To
understand the basic properties of biochemical reactions,however, it is important to study both equilibrium and nonequilibrium characteristics by including forward and backward reactions that satisfy the detailed balance condition.
Such a study is not only important for statistical thermodynamics, but it also provides some insight on the regulation ofsynthesis or degradation reactions for homeostasis in cells.
Recently, we discovered a slow relaxation process to equilibrium, which generally appears in such catalytic reaction
networks, and we proposed “chemical-net glass” as a novel
class of nonequilibrium phenomena. In this case, relaxationin the vicinity of equilibrium is exponential, whereas far
from it much slower logarithmic relaxation with some bottlenecks appears due to kinetic constraints in catalytic relationships 关3兴. In this study, we adopted continuous rate equations
and assumed that the molecule number is sufficiently large.
In biochemical reaction processes, however, some chemical species can play an important role at extremely low con1539-3755/2010/81共5兲/051920共8兲

centrations of even only a few molecules per cell 关4–6兴. In
such systems, fluctuations and discreteness in the molecule
number are important. Indeed, recent studies by using a stochastic simulation of catalytic reaction networks have demonstrated that the smallness in the molecule number induces
a drastic change with regard to statistical and spatiotemporal
behaviors of molecule abundances from those obtained by
the rate equation, i.e., at the limit of large molecule numbers
关7–20兴. In these studies, the strong nonequilibrium conditionis assumed by taking a unidirectional reaction.
Now, it is important to study the relaxation process toequilibrium by considering the smallness in the molecule
number. Does the discreteness in molecule number influence
the equilibrium and relaxation behaviors? Is the relaxation
process slowed down by the smallness in the molecule number? To address these questions, we have carried out several
simulations of the relaxation dynamics of random catalytic
reaction networks by using stochastic simulations. Numericalresults from several networks 关21,22兴 suggest that the relaxation time is prolonged drastically when the number of molecules is smaller. Usually, the temperature dependence of therelaxation time follows the Kramers form exp共␤E兲, where ␤is the inverse temperature and E is a certain energy, depending on the system. When the discreteness in molecule number is relevant to the relaxation, this energy deviates from thevalue at the continuum limit, as the molecule number is decreased. The deviation of the relaxation time from the continuum limit is expressed by the factor exp共␤␦E兲. If ␦E ispositive, it is considered as the effective additional energy
required to pass through the bottleneck due to the discreteness in molecule number.
In this paper, we analyze such slowing down of a reaction
process to equilibrium that is induced by the smallness inmolecule numbers. Instead of taking complex reaction networks, we choose simple networks or network motives toestimate the relaxation time analytically. In fact, complex
networks are often constructed by combining a variety ofsimple network motives with simple branch or loop structures. We focus on the relaxation dynamics of reversible
catalytic reaction systems with such simple network motives
as a first step toward understanding the general relaxation

051920-1©2010 The American Physical Society

PHYSICAL REVIEW E 81, 051920 共2010兲

AKINORI AWAZU AND KUNIHIKO KANEKO(a) Cascade
R1

R2

RS-1

RS

P3

PS-1

PS

R3

RS-1

RS

R3

label a pair of resource and the corresponding product Ri and
Pi by the following manner:
First, we choose and label a resource and the corresponding product as R1 and P1. Then there is one product thatcatalyzes this R1-P1 reaction. We label this product as P2 andthe resource corresponding to this product as R2. Next, theproduct which catalyzes the P2-R2 reaction and the corresponding resource is labeled as R3 and P3, and so forth.
Here, we define the rate of the reaction from Pi to Ri and thatfrom Ri to Pi proceeded by the catalyst Pi+1 as k Pi,Ri andkRi,Pi. If all chemicals are catalyzed by one of them, we canwrite the reaction as

...
P1

P2

(b) Loop
R1

R2

...

k P ,R

P1

P2

P3

PS-1

i

i

kR ,P
i

i

k P ,R
S

S

kR ,P

S

Pi + Pi+1  Ri + Pi+1

PS

FIG. 1. Illustration of 共a兲 cascade and 共b兲 loop systems. Solidarrows indicate reaction paths 共their width indicates the transition
tendency兲 and dashed arrows indicate catalyzation.

for i = 1 , 2 , . . . , S − 1, where
PS + P1  RS + P1 ,S

properties in complex catalytic reaction networks.
In Sec. II, we introduce two network motives, where thesynthesis of a product from resource molecules 共and its reverse reaction兲 is catalyzed by one of the other products.
Here, we note that some specific network motives may exhibit incomplete equilibration when the molecule number decreases, and the average chemical concentration in the steadystate deviates from the equilibrium concentration derived by
the continuous rate equations.
In Sec. III, we show relaxation characteristics from thestochastic simulations. The relaxation of the fluctuation
around the steady state slows down as the molecule number
is decreased below a critical value. The increase in the relaxation time is represented by a scaling function by using
h = N exp共− ␤V兲,
where N is the molecule number and V is the energy gap
between a product and a resource. In Sec. IV, we present ananalytical estimate for this relaxation suppression due to thesmallness in molecule number by using a suitable approximation for master equation. In Sec. V, we present a summaryand discuss the generality of our results.
II. MODELS

Here, we consider reversible catalytic reaction systems
with two simple network structures, cascade and loop systems, as shown in Fig. 1, which may function as network
motives for complex reaction networks. These systems consist of 2S chemical species, S product chemicals, and S resource chemicals; and each product is transformed to thecorresponding resource, and vice versa, by the catalyzationof one of the other products. Here, we assume that eachproduct chemical can catalyze at most one of the otherresource-product reactions. 共Instead, we can interpret that
there exist S chemical species with excited and nonexcitedstates, and chemicals in an excited state can catalyze an excitation reaction of one of the other molecules.兲 Here, we

which leads to the loop system 关Fig. 1共b兲兴. When there exists
a product which does not catalyze any reactions and we label
such a product and the corresponding resource as P1 and R1,
the cascade system in Fig. 1共a兲 is obtained, wherek P ,R
S

S

kR ,P

S

PS  RS .
S

共By neglecting cases in which some pair of resource andproduct is totally disconnected from others, the loop and cascade systems are the only possibilities.兲
The reaction rates k Pi,Ri and kRi,Pi are set so that they satisfy the detailed balance condition. We assume that the energy of the chemical Pi is larger than that of Ri, and set
k Pi,Ri = 1 and kRi,Pi = exp共−␤Vi兲, where Vi is the energy gap
between Pi and Ri and ␤ is the inverse temperature. We
define pi and ri as the numbers of molecules of the chemical
species Pi and Ri, respectively. We fix the total number ofmolecules as SN, and pi + ri = N holds for each i. The state of
the system is represented by a set of numbers 共p1 , p2 , . . . , pS兲.
In both the systems, it is noted that for N → ⬁ 共i.e., the−␤Vi/ 共1 + e−␤Vi兲 andcontinuous limit兲, 具pi典 / N → peq
i /N=eeq−␤Vi
具ri典 / N → ri / N = 1 / 共1 + e 兲 hold at the equilibrium distribution, which is reached at t → ⬁. The derivations of theseequilibrium concentrations from the rate equation will begiven in Sec. IV.
For finite N, however, there is a difference between the
distributions of the cascade and the loop systems. In the cascade system, the average of the equilibrium chemical concentrations is identical to the continuum limit and is given by
具pi典 / N = e−␤V / 共1 + e−␤V兲, that is, they are independent of Nand ␤. This is because all the states 共p1 , p2 , . . . , pS兲 共0 ⱕ pi
ⱕ N兲 are connected by reactions and the above equilibrium
distribution is the only stationary solution for the masterequation.
On the other hand, in the loop system, there is a deviationin the steady chemical concentration from the continuum
limit, which becomes more prominent as N becomes smaller.

051920-2

PHYSICAL REVIEW E 81, 051920 共2010兲

DISCRETENESS-INDUCED SLOW RELAXATION IN…
(a)

0.6

0.4
C(t)

(d)

N=1
N=4

Cascade S = 2
β=3

0.6
C(t)
0.4

N = 16N = 1024

0.2
0
1

10

(b)

0

1

10

100

0.6
C(t)
0.4

0.6
C(t)
0.4
0

t

100

N=1
N=4
N = 16N = 1024
1

10

100

100

1000

10000 tS = 2, β = 3
S = 3, β = 1

C( )S = 3, β = 3
S = 4, β = 1

t

10000

1000
S = 2, β = 1

0.4

0.2
10

t

Loop S = 3
β = 3 CODE
C*

(f)

N=1
N=4
N = 16N = 1024

CODE
1

10

0.8

1000

Cascade S = 2
β=3

0.2

1

0.2

(c) 0.8

0

C*
CODE

N=1
N=4
N = 16N = 1024

0.4
C(t)
0.2

N=4
N = 16N = 1024(e)

Cascade S = 3
β=3

0.6

0

t

100

0.8

C*

N=1

0.2

CODE

CODE

Loop S = 2
β=3

0.8

S = 4, β = 3
0
0.1

1

h

10

FIG. 2. 共Color online兲 C共t兲 of cascade systems with 共a兲 S = 2, 共b兲 S = 3, and 共c兲 S = 4, and loop systems with 共d兲 S = 2 and 共e兲 S = 3 forseveral N with ␤ = 3. 共f兲 C共⬁兲 as a function of h in loop systems for several ␤ and S. CODE indicates the autocorrelations given by Eq. 共4兲
in 共a兲–共c兲 and Eq. 共3兲 in 共d兲 and 共e兲. Cⴱ = exp共−e−␤Vt兲 in 共d兲, and Cⴱ = exp关−共e−2␤V / 2兲t兴 in 共b兲 and 共e兲 with ␤ = 3 共V = 1兲.

This is because the state 共p1 , p2 , . . . , pS兲 = 共0 , 0 , . . . , 0兲 cannot
be reached from other states, whereas the state cannot moveto any other states. Hence, the steady distribution from theinitial conditions without 共p1 , p2 , . . . , pS兲 = 共0 , 0 , . . . , 0兲 deviates from the continuum limit. This deviation becomesprominent as N becomes smaller. For example, for N = 1 and
Vi = V, the distribution from the initial condition without
共p1 , p2 , . . . , pS兲 = 共0 , 0 , . . . , 0兲 is given by 具pi典 = e−␤V
共1 + e−␤V兲S−1 / 关共1 + e−␤V兲S − 1兴. Note that 具pi典 tends to 1 / S withan increase in ␤.
III. SIMULATION RESULTS

In this section, we present the results of stochastic simulations and show the dependence of the relaxation process on
the number of molecules N and the inverse temperature ␤.
For simplicity, we consider Vi to be uniform for all species
共=V兲; however, this assumption can be relaxed.
Numerical simulations are carried out by iterating the following stochastic processes. 共i兲 We randomly pick up a pairof molecules, say, molecules 1 and 2. 共ii兲 Molecule 1 is transformed with its reaction rate 共if it is P, it is transformed to R,
and vice versa兲 if molecule 2 can catalyze the reaction ofmolecule 1. In the cascade case, there is a reaction that
progresses without a catalyst and, in this case, if molecule 1
is the one that reacts without a catalyst, then it is transformed
with the reaction rate independently of molecule 2. Here, aunit time is defined as the time span in which the above
processes for catalytic reactions are repeated SN times. In
each unit time, each molecule is picked up on average tocheck if the transformation occurs. In the following, we focus on the behavior of the system after a sufficiently longtime from the initial time where the numbers of each molecule pi and ri are set randomly from 关0 , N兴 under the con-straint pi + ri = N and 共p1 , p2 , . . . , pS兲 ⫽ 共0 , 0 , . . . , 0兲.
Now, we define the autocorrelation function c共t兲 as c共t兲
eq
eq
eq
by
= 具兺i兵关pi共t兲 − peqi 兴关pi共0兲 − pi 兴 + 关ri共t兲 − ri 兴关ri共0兲 − ri 兴其典,scaling so that c共t兲 = 0 when the equilibrium is reached. Byfurther normalizing the function, so that it takes unity in theinitial distribution, we define C共t兲 ⬅ c共t兲 / c共0兲. Figures
2共a兲–2共e兲 show C共t兲 of the cascade system 关Figs. 2共a兲–2共c兲兴
and the loop system 关Figs. 2共d兲 and 2共e兲兴 for some S and Nwith ␤ = 3. As already discussed, C共⬁兲 → 0 in the cascadesystem, whereas C共⬁兲 → 0 for large N but C共⬁兲 ⬎ 0 for small
N in the loop system. Here, the value C共⬁兲 in the loop system starts to deviate when h = Ne−␤V becomes less than 1.
Hence, we have plotted C共⬁兲 of the loop system as a function of h in Fig. 2共f兲 for ␤ = 1 and 3. As shown, C共⬁兲 ⬎ 0
holds for h ⬍ 1 independently of ␤. On the other hand, inboth systems, the relaxation to the equilibrium value C共⬁兲 isdrastically slowed down for small N, as compared to that forlarge N when S ⬎ 2, whereas the relaxation for small N isfaster when S = 2.
To observe the dependence of the relaxation time on N,
we measured the integrated relaxation time defined as ␶
= 兰⬁0 兵关C共t兲 − C共⬁兲兴 / 关1 − C共⬁兲兴其dt. Figures 3共a兲 and 3共b兲 show
␶ as a function of N for ␤ = 3 with S = 2 , 3 , 4 for the cascade关Fig. 3共a兲兴 and loop 关Fig. 3共b兲兴 systems. For S ⱖ 3, the relaxation time ␶ increases by several orders of magnitude with adecrease in N in both systems. On the other hand, ␶ for S= 2 does not exhibit any drastic change with the decrease in
N in both systems.
This prolongation of ␶ for S ⬎ 2 becomes more prominentas ␤ is increased. From several data, ␶ is suggested to increase as a function of exp共␤V兲. Combining N and ␤ dependencies, we introduce a parameter h = N exp共−␤V兲. The discreteness effect is dominant when h = N exp共−␤V兲 is lessthan unity. Figures 3共c兲 and 3共d兲 show ␳ = ␶ / ␶N→⬁ as a function of h for the cascade 关Fig. 3共c兲兴 and loop 关Fig. 3共d兲兴

051920-3

PHYSICAL REVIEW E 81, 051920 共2010兲

AKINORI AWAZU AND KUNIHIKO KANEKO
(a)

(c)

Cascade β = 3

3

10

τ

CascadeρτC (2)
τC (3)
τC (4)

S=2
S=3
S=4

For (c) and (d)S = 2, β = 1
S = 2, β = 2
S = 2, β = 3
S = 2, β = 4

10

100
1

10

1

(b) 4
10

10

100
N

1000

(d)
100

Loop β = 3

τ

τ (2)
τL (3)
τL (4)

1

10

100

h

S = 4, β = 1
S = 4, β = 2
S = 4, β = 3
S = 4, β = 4

Loop

ρ

L

S=2
S=3
S=4

3

10

0.1

10000

S = 3, β = 1
S = 3, β = 2
S = 3, β = 3
S = 3, β = 4

10

1

100
10
1

0.1

10

100
N

1000

1

10

100

10000

h

FIG. 3. 共Color online兲 ␶ as a function of N in 共a兲 cascade and 共b兲 loop systems, and ␶C共S兲 and ␶L共S兲 for ␤ = 3 and S = 2 , 3 , 4. ␳ as afunction of h in 共c兲 cascade 共d兲 loop systems with S = 2 , 3 , 4 for several ␤.

systems for several values of ␤ and S = 2 , 3 , 4. For S ⬎ 2, thedeviation of ␳ from the continuum limit 共␳ = 1兲 becomesprominent when h is below unity in both systems. The increase in ␳ appears to become steeper with an increase in S.
On the other hand, ␳ for S = 2 does not exhibit a drastic
increase with a decrease in h.
IV. ORIGIN OF SLOW RELAXATIONS AND CROSSOVER
A. Relaxation processes for N \ ⴥ and N = 1

Now, we analytically estimate the enhancement in relaxation time and explain its representation in the form h= N exp共−␤V兲. For this purpose, we compare the estimate bymaster equation analysis for small N and compare it with thatfrom the continuum limit N → ⬁.
In the continuum limit, the reaction dynamics are represented by the following rate equation:冋 冉 冊 册ẋi = xc e−␤V1− xi − xi ,S

共1兲with xi = pi / SN. Here, xc = 1 for i = S in the cascade system,
xc = x1 for i = S in the loop system, and xc = xi+1 for i ⬍ S in−␤V/ 关S共1 + e−␤V兲兴 holds for t → ⬁.
both systems. xi → xeq
i =eWhen the deviation from equilibrium ␦xi = xi − xeq
i is small, its
evolution for the loop systems obeys the following linearized
equation:

␦ẋi = −

e −␤V
␦xi .
S

共2兲

For the cascade system, this equation is also valid for theelements i ⫽ S, whereas ␦ẋS = −␦xS. Then, the autocorrelation
function of a small fluctuation of pi around peq
i is obtained as冉

C共t兲 = exp −for the loop system, and

e −␤V
t
S

冊
冉

1
S−1
e −␤V
C共t兲 = exp共− t兲 +
exp −
t
S
S
S

共3兲

冊

共4兲

for the cascade system. Indeed, these agree quite well withthe simulation results for a sufficiently large N 共e.g., N= 1024 in Fig. 2.兲. Thus, the characteristic time of the relaxation is estimated as ␶L共S兲 ⬃ Se␤V for the loop system and␶C共S兲 ⬃ S1 + 共S − 1兲e␤V for the cascade system, which are consistent with the simulation results shown in Fig. 3.
As the other extreme limit, consider the case with N = 1.
In this case, the relaxation dynamics are dominated by acompletely different process induced by the absence of catalysts whose number can often go to zero. In such cases,
states are trapped at some local energy minimum that appears due to the deficiency of catalysts. Then, the hopping
processes among them play an important role in the relaxation dynamics, as shown below. In the following, we focus
on the cases with S = 2 and 3 to clarify that such an effect is
induced by discreteness in the molecule number. Note that,as shown in the last section, the behavior for S ⱖ 3 is distinct
from that for S = 2; in the former case, the relaxation time isenhanced by the decrease in N, in contrast to the latter case.
First, we study the loop system. When S = 2, the system
realizes three states from the initial conditions—共p1 , p2兲
= 共1 , 0兲, 共0,1兲, and 共1,1兲—as shown in Fig. 4共a兲. In this case,the transition rate from the state 共1,0兲 to 共1,1兲 is estimated asfollows: for this transition, a pair of molecules from theproduct of the first species and the resource of the second
1,
species has to be chosen. This probability is given by 21 2−1−␤V
while the reaction rate is given by e . In unit time, this051920-4

PHYSICAL REVIEW E 81, 051920 共2010兲

DISCRETENESS-INDUCED SLOW RELAXATION IN…
(a) Loop, S = 2( p1 , p 2 )(c) Loop, S = 3( p1 , p2 , p 3 )(0, 1)

e-βV
1

1
e βV

(1, 1)

the characteristic time of the correlation of each pi is givenby ⬃e␤V.
On the other hand, for S = 3, the system realizes eightstates—共p1 , p2 , p3兲 = 共0 , 0 , 0兲, 共1,0,0兲, 共0,1,0兲, 共0,0,1兲, 共1,1,0兲,
共1,0,1兲, 共0,1,1兲, and 共1,1,1兲—as shown in Fig. 4共d兲. The
slowest characteristic time of the relaxation is given by thetransition time from the branch 共1,1,0兲-共0,1,0兲 since the system must go through the highest-energy state 共1,1,1兲, which
is a limiting process for this case. Then, in a manner similar
to the loop system with S = 3, the characteristic time is obtained as ⬃2e2␤V. This gives the characteristic time of theslowest motions of the system. This estimation fits well with
the numerical result shown in Fig. 2共b兲.

(1, 0)

1/2(1, 0, 1)(1, 0, 0)βV
e βV/2
1/2 e /2
1/2
1/2
(1, 1, 1)(0, 1, 0)(1, 1, 0)βVe /2e βV/2e-βV/2 1/2
1/2(0, 0, 1)(0, 1, 1)e-βV/2(b) Cascade, S = 2 ( p1 , p 2 )(0, 1)e-βV
1(1, 1)1e βV(1, 0)e-βV(0, 0)

1

B. N and ␤ dependencies of C(ⴥ) and relaxation time(d) Cascade, S = 3 ( p1 , p2 , p 3 )(1, 0, 1)
1/2 e-βV/2(1, 1, 1)e-βV/2 1/2

e-βV/21/2
e-βV/2(1, 0, 0)1/2e-βV/2(0, 0, 0)

1/2

1/2
1/2
(0, 0, 1)(0, 1, 1)βVe /2e βV/2(1, 1, 0)

1/2(0, 1, 0)e βV/2

FIG. 4. 共a兲 Illustration of transition diagrams of 共a兲 loop systemwith S = 2, 共b兲 loop system with S = s, 共c兲 cascade system with S= 2, and 共d兲 cascade system with S = 3, where arrows indicate possible transitions and the values next to them specify the transition
ratios.

process is iterated two times. Hence, the rate is given by 2
1 −␤V⫻ 21 2−1e = e−␤V. Then, the characteristic time of the correlation of each pi is given by ⬃e␤V, which is consistent withthe results shown in Fig. 2共d兲.
On the other hand, for S = 3, the system realizes sevenstates—共p1 , p2 , p3兲 = 共1 , 0 , 0兲, 共0,1,0兲, 共0,0,1兲, 共1,1,0兲, 共1,0,1兲,
共0,1,1兲, and 共1,1,1兲—as shown in Fig. 4共b兲. The characteristic time of the correlation of each pi is given by the transition
time among the three branches including lowest-energy
states, 共1,0,1兲-共1,0,0兲, 共1,1,0兲-共0,1,0兲, and 共0,1,1兲-共0,0,1兲.
Here, in order to hop from one branch to another, the system
must go through the highest-energy state 共1,1,1兲 due to therestriction by the catalytic relation. Thus, the escape ratefrom each branch is estimated by ⬃ 21 e−2␤V, and the characteristic time of the correlation of each pi is estimated as⬃2e2␤V 共see Appendix A兲. Because the relaxation time in thecontinuum limit is proportional to exp共␤V兲, the deviation ␳from it increases with exp共␤V兲, which is consistent with theresults shown in Fig. 2共e兲. Thus, the enhancement of therelaxation time from the continuous case is explained.
Essentially the same argument is also valid for the cascade systems. When S = 2, the system can realize transitions
among four states—共0,1兲-共1,1兲-共1,0兲-共0,0兲—as shown in Fig.
4共c兲. Here, 共0,1兲 is a metastable state and 共0,0兲 is the lowestenergy state. The relaxation is characterized by the escaperate from a metastable state, which is given by ⬃e−␤V. Thus,

Next, we extend the argument of the last section to analyze the N and ␤ dependencies of C共⬁兲 and the relaxation
time in greater detail. In particular, we explain why h= N exp共␤V兲 ⬃ 1 gives a critical value and how the amplification of relaxation time depends on h for h ⬍ 1. Because of
the simplicity due to the symmetry in the catalytic relationship, we only study loop systems; however, the argument
presented below can be extended to cascade systems.
Figure 5共a兲 shows the transition diagram of the loop system with S = 2, where each circle indicates each state 共p1 , p2兲
and the arrows indicate possible transitions. Generally, forany values of S, the transition rate from a state
共p1 , p2 , . . . , pi = n , pi+1 , . . . , pS兲 to a state 共p1 , p2 , . . . , pi = n+ 1 , pi+1 , . . . , pS兲 per unit time is estimated as follows. For
this transition, a pair of molecules from the resource of theith species 共Ri兲 and the product of the 共i + 1兲th species 共Pi+1兲
has to be chosen. This probability is given by 关共N− pi兲 / SN兴关pi+1 / 共SN − 1兲兴, and the reaction rate is given bye−␤V. Hence, the transition rate per unit time is given byi
Wn→n+1=

共N − n兲pi+1 −␤Ve .
SN − 1

Similarly, the transition rate in the opposite direction is given
asi=
Wn+1→n

共n + 1兲pi+1.
SN − 1

If the molecule number is so large or ␤ is so small thatii
⬎ Wn+1→nholds for small n and
h = Ne−␤V Ⰷ 1, Wn→n+1
ii
Wn→n+1 ⬍ Wn+1→n holds for large n. Then, the dominantstates of the system are located in an intermediate region in
the phase space 关0 , N兴. For example, the intermediate gray
共blue online兲 region in Fig. 5共a兲 indicates such dominantstates for S = 2. Under such conditions, xi = 具pi典 / SN obeys Eq.
共1兲 for a sufficiently large value of N 共see Appendix B兲.
On the other hand, if N is so small or ␤ is so large thatii
Ⰶ Wn+1→nholds for all i and n. Thus, pi for allh Ⰶ 1, Wn→n+1
i tends to decrease to zero. Then, there exist SN metastablestates—共n , 0 , 0 , . . . , 0兲,
共0 , n , 0 , . . . , 0兲,. . .,
共0 , 0 , . . . , 0 , n , 0 , . . , 0兲, . . ., and 共0 , 0 , 0 , . . . , n兲 共1 ⱕ n ⱕ N兲.
Among them, the following S states, 共1 , 0 , 0 , . . . , 0兲,
共0 , 1 , 0 , . . . , 0兲, . . ., and 共0 , 0 , 0 , . . . , 1兲, have the lowest en-

051920-5

PHYSICAL REVIEW E 81, 051920 共2010兲

AKINORI AWAZU AND KUNIHIKO KANEKO
(a)

N

β=2

Loop S = 2

β=3

β=4

Simulation(a) S = 2

τLd (2)
100

τ

p1

10

exp(2) exp(3) exp(4)

2
1

1

10

100

Loop S = 3
Simulation

β=2

1000

N

0
0

1

2

(b) S = 2 (p1, p2)(1,0)

b = 1/(SN-1)
c = exp(-βV)(SN-1)Np2(b)2b (N-1)c(N-1)c
(1,1)
b
2b
b Nc

Nc

β=3

β=4

τLd (3)

1000

τ

(0,1)

100
(c) S = 3(N-1)c 2b(p1, p2, p3)(N-1)c(N-1)c2b2b10(1,1,1)(N-1)c

Nc

2b

b Nc

(1,0,1)Ncb(0,1,1)Ncb(1,0,0)b

Nc (N-1)c

Nc

2b
b

(0,0,1)2b(1,1,0)(N-1)cb(0,1,0)

FIG. 5. 共Color online兲 共a兲 Illustration of the transition diagramsfor S = 2, and effective transition diagrams for 共b兲 S = 2 and 共c兲 S= 3, where bold arrows indicate the focused transitions in the text.

ergy. For example, in the cases with S = 2, the states 共0 , p兲
and 共q , 0兲 共p , q ⫽ 0兲 are metastable states and 共1,0兲 and 共0,1兲
are the lowest-energy states.
It should be noted that the lowest-energy states are thedominant states for h Ⰶ 1. The probability to realize theselowest-energy states tends to 1 / S with an increase in ␤.
Thus, with the increase in ␤, 具pi典 approaches 1 / S for smallN, which indicates C共⬁兲 = const⬎ 0 for small N and large ␤.
Moreover, for h Ⰶ 1, the transitions among lowest-energy
states contribute dominantly to the relaxation process. Then,we estimate the characteristic time of the fluctuations of thesystem for h Ⰶ 1 by considering the transition processes fromone lowest-energy states such as 共0 , 0 , . . . , 0 , p j = 1 , 0 , . . . , 0兲
to the other lowest-energy states such as 共0 , 0 , . . . , 0 , p j= 0 , 0 , p j⬘ = 1 , 0 , . . . , 0兲. In the following, we consider only
the cases with S = 2 and 3. We only focus on the dynamics ofp j under the constraint that p j has only 0 or 1, because hⰆ 1.
First, consider the case with S = 2. Figure 5共b兲 shows adetailed transition diagram around the region where pi
共i = 1 , 2兲 is only 0 or 1. The escape rate from 共1,0兲 and 共0,1兲
N
is given by ⬃ 2N−1
e−␤V. Thus, the characteristic time of thecorrelation of each pi is given by

1

10
100
exp(2) exp(3) exp(4)

1000

N

FIG. 6. 共Color online兲 Relaxation time ␶ obtained from simulations 共points兲 and its approximate analytical expression ␶Ld 共S兲
共curves兲 estimated in the text. Plotted as a function of N for the loop
systems with 共a兲 S = 2 and 共b兲 S = 3 with ␤ = 2 , 3 , 4. The analyticalexpression agrees with the simulation data both for small N, and forlarge N, where ␶ approaches a constant value expected from the rateequation. The crossover occurs at around h = Ne−␤V ⬃ 1
关N ⬃ exp共2兲, exp共3兲, and exp共4兲 for ␤ = 2, 3, and 4兴.

␶Ld 共2兲 ⬃2N − 1 ␤V
e ,
N

共5兲

which is consistent with the results shown in Fig. 6共a兲.
Next, we study the case with S = 3. The transition diagramof the states 共p1 , p2 , p3兲 is shown in Fig. 5共c兲 when pi
共i = 1 , 2 , 3兲 takes only 0 or 1. Similar to the N = 1 case, thecharacteristic time of the transition among the three branches
including lowest-energy states, 共1,0,1兲-共1,0,0兲, 共1,1,0兲共0,1,0兲, and 共0,1,1兲-共0,0,1兲 through the state 共1,1,1兲, is considered. In a manner similar to the N = 1 case, the transition
rate from each branch is estimated by ⬃关Ne−␤V / 共3N− 1兲兴关Ne−␤V / 共1 + Ne−␤V兲兴 = N2e−␤2V / 关共3N − 1兲共1 + Ne−␤V兲兴.
Thus, the relaxation time of the fluctuation of p1 is estimatedas the decrease with N as

␶Ld 共3兲 ⬃共3N − 1兲共1 + Ne−␤V兲 2␤V
e .
N2

共6兲

Considering the e␤V dependence of ␶N→⬁, the above estimate
is consistent with Fig. 6共b兲.

051920-6

PHYSICAL REVIEW E 81, 051920 共2010兲

DISCRETENESS-INDUCED SLOW RELAXATION IN…

For S larger than 3, the transition diagram becomes rather
complicated. However, a similar analysis should be possible
to estimate the prolongation in the relaxation time.
V. SUMMARY AND DISCUSSIONSIn the present paper, the slowing down of the relaxation inreversible catalytic reaction networks induced by the smallness of molecule number is investigated as a general property of catalytic reaction networks. This prolongation of relaxation is a result of bottlenecks in reactions; these appear
due to the deficiency of the catalyst required for a reaction.
The number of molecules can be so small that the number ofcatalysts becomes zero. In this case, a pair of a substrate andthe corresponding catalyst molecule species can hardly exist
simultaneously. Such a constraint makes it difficult to realize
a specific configuration necessary for the relaxation. The
probability for realization is given by exp共−␤Eb兲, with Eb as
the corresponding effective energy barrier to realize suchrare conditions, or the sum of such energy barriers. Thisbottleneck energy is generally different from the energy gap
in the continuum limit that is obtained from the relaxation
time of the rate equation 共ordinary differential equation兲.
Thus, the relaxation time at a small molecule number deviates from the continuum case by the factor exp共␤␦E兲 with anappropriate effective energy difference ␦E.
By considering the models of simple catalytic reaction
networks consisting of resource chemicals of S species andthe corresponding products, we have demonstrated this deviation of relaxation time from both direct simulations andanalysis by using master equation. From the numerical andanalytical estimates, V ⱕ Eb ⱕ 2V and 0 ⱕ ␦E ⱕ V for S = 3,
where V is the energy gap between the resource and the
product chemicals. For S ⬎ 2, in general, the prolongation of
the relaxation time becomes prominent when h = N exp
共−␤V兲 is less than unity, and its amplification ratio from thecontinuum limit is represented as a function of S and h. Notethat the cascade system in the N = 1 case is equivalent to the“asymmetrically constrained Ising chain” 共ACIC兲, hierarchically constrained Ising model, or east model, which are studied as simple abstract models for glassy states 关23–25兴. Following the interpretation therein, the increase in relaxation
time at h ⬍ 1 as a result of the decrease in N or temperature
may be regarded as a type of glass transition. According to
the recent studies on ACIC, the correlation time of the motion of p1 共not the relaxation time of the total system兲 isestimated as ␶1 ⬃ 共1 + 2e␤V兲k, where the integer k obeys
2共k−1兲 ⬍ S ⱕ 2k 关24,25兴. In cases with S = 2 , 3 , 4, this fact isconsistent with our estimate of the relaxation time of thecascade system with N = 1. The estimation of ␦E as a function of S and h for general cases both for cascade and loop
systems is an important issue that should be studied in the
future.
In addition to the slowdown in relaxation, the equilibrium
distribution deviates in a network called a loop system,
where all the reactions are catalyzed by one of the products.
The constraint that the numbers of a certain pair of chemical
species cannot simultaneously be zero leads to the deviationof the average distribution of molecule numbers from thecontinuum limit. Again, this deviation becomes prominent
when h is less than unity.
Although we have adopted simple network motives toanalyze the relaxation, the prolongation of relaxation time isquite general in catalytic reaction networks. Catalytic bottlenecks often appear as the number of molecules is decreased
in a large variety of reaction networks in which catalysts aresynthesized within 关21,22兴. The present study can provide abasis for the general case with complex networks, as the
motives here are sufficiently small and can exist within suchcomplex networks.
Biochemical reactions generally progress in the presenceof catalysts that are themselves synthesized as products ofsuch reactions. These reactions form a network of a varietyof chemical species. Here, the molecule number of each species is generally not very large. Hence, the slow relaxation
process and deviation from equilibrium discussed in thisstudy may underlie intracellular reaction processes. Moreover, the present network motives are so simple that they aresuggested to exist in biochemical networks. We also note thatthe resource and product in our model can be interpreted as
nonexcited and excited states of enzymatic molecules. Indeed, many molecules are known to exhibit catalytic activity
only when they are in an excited state, which can help otherchemicals to switch to an excited state. In fact, such networks with mutual excitation are known in signaltransduction networks 关26–28兴, where the present slow relaxation mechanism may be relevant to sustain the excitabilityof a specific enzyme type over a long time span. It is important to pursue the relevance of the present mechanism incell-biological problems by considering more realistic models in the future.
We also note that not only the discreteness in the molecule number but also the negative correlation between asubstrate and the corresponding catalyst within a reaction
network or in a spatial concentration pattern suppresses therelaxation process 关2,3,21,22兴. The present mechanism due
to discreteness may work synergetically with the earlier
mechanism to further suppress the relaxation to equilibrium.
The construction of reaction networks to achieve slower relaxation together with the network analysis will be an important issue in the future.
ACKNOWLEDGMENTS

The authors would like to thank Shinji Sano for informing
us of his finding on the prolongation of relaxation in reaction
networks due to the discreteness in molecule numbers, which
triggers the present study. A.A. was supported in part by aGrant-in-Aid for Young Scientists 共B兲 共Grant No. 19740260兲
from MEXT of Japan.
APPENDIX A

We define the probability that the states in the branch
共1,0,1兲-共1,0,0兲 are realized as Q1,0,0. Then, the probability torealize the state 共1,0,1兲 is given by 关e−␤V / 共1 + e−␤V兲兴Q1,0,0.
Here, the transition rate from 共1,0,1兲 to 共1,1,1兲 is given bye−␤V / 2. Then, the probability current from the branch051920-7

PHYSICAL REVIEW E 81, 051920 共2010兲

AKINORI AWAZU AND KUNIHIKO KANEKO

共1,0,1兲-共1,0,0兲
isestimatedby⬃共e−␤V / 2兲关e−␤V /
−␤V
−2␤V
−␤V
Q1,0,0 共e Ⰶ 1兲. Because of the
共1 + e 兲兴Q1,0,0 ⬃ 1 / 2esymmetry among the catalytic reactions, the probability currents from the other branches are obtained in the same way
to get the same form. Thus, the escape rate from each branch
is estimated by ⬃ 21 e−2␤V.
APPENDIX B

We define the probability that pi = n as Qin, and the jointi. Here, Qinprobability to realize pi = n and pi+1 = m as Qn,m
N
i
i
i i+1
= 兺m=0Qn,m and Qn,m = QnQm . Then the time evolution ofifollows
Qn,m
i
=
Q̇n,m

Q̇in =i+ 共n + 1兲Qn+1− nQin − 共N − n兲e−␤VQin其,

共B1兲

共B2兲

i+1
N
N
nQin 共具pi+1典 = 兺m=0
mQm
兲. Using this equawhere 具pi典 = 兺n=0tion, we obtain the time evolution of 具pi典 as

具ṗi典 =mii
兵关N − 共n − 1兲兴e−␤VQn−1,m+ 共n + 1兲Qn+1,m
SN − 1ii
其.
− nQn,m− 共N − n兲e−␤VQn,m

具pi+1典i
兵关N − 共n − 1兲兴e−␤VQn−1
SN − 1

具pi+1典关− 具pi典 + 共N − 具pi典兲e−␤V兴.
SN − 1

共B3兲

Then, we obtain

This implies that xi = 具pi典 / SN obeys Eq. 共1兲 for a sufficiently
large value of N.

关1兴 K. Kang and S. Redner, Phys. Rev. Lett. 52, 955 共1984兲; M.
Yang, S. Lee, and K. J. Shin, ibid. 79, 3783 共1997兲; I. V.
Gopich, A. A. Ovchinnikov, and A. Szabo, ibid. 86, 922
共2001兲; D. Pines and E. Pines, J. Chem. Phys. 115, 951 共2001兲.
关2兴 A. Awazu and K. Kaneko, Phys. Rev. Lett. 92, 258302 共2004兲.
关3兴 A. Awazu and K. Kaneko, Phys. Rev. E 80, 041931 共2009兲.
关4兴 N. Olsson, E. Piek, P. ten Dijke, and G. Nilsson, J. Leukoc
Biol. 67, 350 共2000兲.
关5兴 P. Guptasarma, BioEssays 17, 987 共1995兲.
关6兴 H. H. McAdams and A. Arkin, Trends Genet. 15, 65 共1999兲.
关7兴 Y. Togashi and K. Kaneko, Phys. Rev. Lett. 86, 2459 共2001兲;J. Phys. Soc. Jpn. 72, 62 共2003兲; J. Phys.: Condens. Matter 19,
065150 共2007兲.
关8兴 J. Ohkubo, N. Shnerb, and D. A. Kessler, J. Phys. Soc. Jpn.
77, 044002 共2008兲.
关9兴 A. Awazu and K. Kaneko, Phys. Rev. E 76, 041915 共2007兲.
关10兴 A. Awazu and K. Kaneko, Phys. Rev. E 80, 010902共R兲 共2009兲.
关11兴 B. Hess and A. S. Mikhailov, Science 264, 223 共1994兲; J.
Theor. Biol. 176, 181 共1995兲.
关12兴 P. Stange, A. S. Mikhailov, and B. Hess, J. Phys. Chem. B
104, 1844 共2000兲.
关13兴 N. M. Shnerb, Y. Louzoun, E. Bettelheim, and S. Solomon,
Proc. Natl. Acad. Sci. U.S.A. 97, 10322 共2000兲.

关14兴 Y. Togashi and K. Kaneko, Physica D 205, 87 共2005兲.
关15兴 G. Marion, X. Mao, E. Renshaw, and J. Liu, Phys. Rev. E 66,
051915 共2002兲.
关16兴 V. P. Zhdanov, Eur. Phys. J. B 29, 485 共2002兲.
关17兴 T. Dauxois, F. Di Patti, D. Fanelli, and A. J. McKane, Phys.
Rev. E 79, 036112 共2009兲.
关18兴 K. Kaneko, Adv. Chem. Phys. 130, 543 共2005兲.
关19兴 C. Furusawa and K. Kaneko, Phys. Rev. Lett. 90, 088102
共2003兲.
关20兴 C. Furusawa et al., Biophysics 1, 25 共2005兲.
关21兴 S. Sano, Master thesis, University of Tokyo, 2009 共in Japanese兲.
关22兴 S. Sano, A. Awazu, and K. Kaneko 共unpublished兲.
关23兴 J. Jäckle and S. Eisinger, Z. Phys. B: Condens. Matter 84, 115
共1991兲.
关24兴 F. Mauch and J. Jäckle, Physica A 262, 98 共1999兲.
关25兴 P. Sollich and M. R. Evans, Phys. Rev. Lett. 83, 3238 共1999兲;Phys. Rev. E 68, 031504 共2003兲.
关26兴 A. Goldbeter and D. K. Koshland, Jr., Proc. Natl. Acad. Sci.
U.S.A. 78, 6840 共1981兲; J. Biol. Chem. 259, 4441 共1984兲.
关27兴 A. Levchenko and P. A. Iglesias, Biophys. J. 82, 50 共2002兲.
关28兴 W. Ma, A. Trusina, H. El-Samad, W. A. Lim, and C. Tang, Cell
138, 760 共2009兲.

051920-8