Molecular Synchronization Waves in Arrays of Allosterically Regulated Enzymes
Vanessa Casagrande,1 Yuichi Togashi,2, ∗ and Alexander S. Mikhailov2, †
1Hahn-Meitner-Institut, Glienicker Straße 100, 14109 Berlin, Germany
2Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, 14195 Berlin, Germany
Spatiotemporal pattern formation in a product-activated enzymic reaction at high enzyme con-
centrations is investigated. Stochastic simulations show that catalytic turnover cycles of individual
enzymes can become coherent and that complex wave patterns of molecular synchronization can
develop. The analysis based on the mean-field approximation indicates that the observed patterns
result from the presence of Hopf and wave bifurcations in the considered system.
PACS numbers: 82.40.Ck, 87.18.Pj, 82.39.Fk, 05.45.Xt
Molecular machines, such as molecular motors, ion
pumps and some enzymes, play a fundamental role in
biological cells and can be also used in the emerging soft-
matter nanotechnology [1]. A protein machine is a cyclic
device, where each cycle consists of conformational mo-
tions initiated by binding of an energy-bringing ligand
[2, 3]. In motors, such internal motions generate me-
chanical work [4], while in enzymes they enable or facil-
itate chemical reaction events (see, e.g., [5, 6]). Much
attention has been attracted to studies of biomembranes
with ion pumps and molecular motors, where membrane
instabilities and synchronization effects have been ana-
lyzed [7, 8, 9]. Here, a different class of distributed ac-
tive molecular systems — formed by enzymes — is con-
sidered. The catalytic activity of an allosteric enzyme
protein is activated or inhibited by binding of small reg-
ulatory molecules; the role of such regulatory molecules
can be played by products of the same reaction [10]. Pre-
vious investigations of simple product-regulated enzymic
systems [11, 12] and enzymic networks [13] in small spa-
tial volume with full diffusional mixing have shown that
spontaneous synchronization of molecular turnover cycles
can take place there. External molecular synchroniza-
tion of enzymes of the photosensitive P-450 dependent
monooxygenase system by periodic optical forcing has
been experimentally demonstrated [14].
In this Letter, spatiotemporal pattern formation in en-
zymic arrays is investigated. In such systems, immobile
enzymes are attached to a solid planar support immersed
into a solution through which fresh substrate is supplied
and product molecules are continuously removed. Prod-
uct molecules released by an enzyme diffuse through the
solution and activate catalytic turnover cycles of neigh-
bouring enzymes in the array.
A simple stochastic model [12] of an enzyme as a cyclic
machine (a stochastic phase oscillator), shown in Fig. 1,
is used. Binding of a substrate molecule to an enzyme i
initiates an ordered internal conformational motion, de-
scribed by the conformational phase coordinate φi. The
initial state corresponds to the phase φi = 0. The cat-
alytic conversion event takes place and the product is
released at the state φp inside the cycle. After that, the
substrate
enzyme
regulatory
molecule
feedback
product
FIG. 1: (Color online) A sketch of the model.
conformational motion continues until the equilibrium
state of the enzyme (φi = 1) is finally reached. Initi-
ation of a turnover cycle is a random event, occurring
at a certain probability rate. We assume that substrate
is present in abundance, and its concentration is not af-
fected by the reactions. Conformational motion inside
the cycle is modeled as a stochastic diffusional drift pro-
cess, described by equation
φi = v+ ηi(t), where v is the
mean drift velocity and ηi(t) is an internal white noise
with 〈ηi(t)ηj(t
′)〉 = 2σδijδ(t− t
′) where σ specifies inten-
sity of intramolecular fluctuations.
Allosterically activated enzymes possess a site on their
surface where regulatory molecules can become bound.
Binding of a regulatory molecule leads to conformational
change that enhances catalytic activity of the enzyme. A
regulatory molecule binds to an enzyme with rate con-
stant β and dissociate from it with rate constant κ. Bind-
ing of a regulatory molecule at an enzyme raises its prob-
ability to start a cycle from α0 to α1. We assume that
a regulatory molecule can bind to an enzyme only in its
rest state and this molecule is released when the cycle
is started. The role of regulatory molecules is played
by product molecules of the same reaction. Immobile
enzymes are randomly distributed in space with concen-
tration c. Product diffuses at diffusion constant D and
undergoes decay at rate constant γ. The characteristic
diffusion length of product molecules is ldiff =
In our stochastic 2D simulations, the medium was dis-
cretized into spatial cells (up to 256 × 256), each con-
http://arxiv.org/abs/0704.0021v2
FIG. 2: Stochastic (a,b) and mean-field (c,d) simulations of
2D wave patterns; (a) τp = 0.14, c = 1, and β = 300, (b)
τp = 0.25, c = 10, and β = 10, (c) τp = 0.14, c = 1, and
β = 300, (d) τp = 0.34, c = 100, and β = 1.42. Other
parameters are α0 = 1, α1 = 1000, κ = 10, γ = 10, σ = 0,
D = 100. The linear size of the shown area is L = 40 ldiff in
all panels.
taining a number of enzyme molecules. The cells were
so small that diffusional mixing of product molecules in
a cell within the shortest characteristic time of the reac-
tion could always take place. Each enzyme was described
by the stochastic model given above; diffusion of product
molecules was modeled as a random walk over a discrete
cell lattice. The mean cycle time τ = 1/v was chosen as
the time unit (τ = 1). Systems including up to 655 360
enzymes were used in the simulations.
Figure 2a,b (see also Videos 1 and 2 in ref. [15]) shows
two typical examples of stochastic 2D simulations. Here,
spatial distributions of product molecules are displayed.
Waves of product concentration are propagating through
the medium. In a peak of a wave, many locally present
enzymes are simultaneously releasing product molecules.
Since product release can take place only at a certain
stage inside the cycle, this means that the cycles of en-
zymes are locally synchronized. Not only regular wave
structures, such as rotating spiral waves or target pat-
terns (Fig. 2a), but also complex regimes of wave turbu-
lence (Fig. 2b) have been observed.
To understand and interpret stochastic simulation re-
sults, an analytical study of the system in the mean-field
approximation, which holds in the limit of high enzyme
concentrations, has been performed. In this approxima-
tion, the system is characterised by three continuous vari-
ables n0(r, t), n1(r, t) and m(r, t) which represent local
concentrations of enzymes in the rest state without or
with regulatory molecules attached (n0 and n1) and local
concentration of the product (m). For simplicity, internal
fluctuations in enzymes are neglected (σ = 0). Thus, all
enzymes which have started their cycles at some time t
would release their products at a definite time t+τp (with
τp = φp/v) and finish their cycles, returning to the rest
state, at time t + τ . Therefore, the system is described
by a set of three reaction-diffusion equations with time
delays,
= βmn0 − κn1 − α1n1 (1a)
= −βmn0 + κn1 − α0n0 + α0n0(t− τ)
+α1n1(t− τ) (1b)
= −βmn0 + κn1 + α1n1 − γm+ α0n0(t− τp)
+α1n1(t− τp) +D∇
2m. (1c)
The system always has a uniform stationary state
with certain concentrations n0, n1 and m, which can be
found as solutions of the respective algebraic equations.
This state corresponds to the absence of synchronization.
However, it may become unstable if allosteric activation
is strong enough. To analyze stability, small perturba-
tions δn0, δn1 and δm are added to the stationary state,
equations (1) are linearized and their solutions are sought
as δn0 ∼ δn1 ∼ δm ∼ exp (λqt− iqx) with λq = µq+iωq.
Thus, each spatial mode with wavevector q is character-
ized by its frequency ωq and its rate of growth µq. The
properties µq and ωq are given by the roots of a charac-
teristic equation which is determined by the linearization
matrix of equations (1). The steady state becomes unsta-
ble when at least one spatial mode with some wavenum-
ber q0 starts to grow (µq0 > 0).
As the bifurcation parameter, coefficient β can be cho-
sen. If regulatory molecules cannot bind to enzymes
(β = 0), feedback is absent and instabilities are not pos-
sible. On the other hand, allosteric activation becomes
strong if regulatory molecules can easily bind and, in this
case, emergence of oscillations and wave patterns can be
expected. Our bifurcation analysis reveals that, depend-
ing on the parameters of the system, it can exhibit ei-
ther a Hopf or a wave bifurcation [16]. As a result of
the Hopf bifurcation, uniform oscillations with q = 0 de-
velop. Because of the presence of delays in equations (1),
the characteristic equation is nonpolynomial in terms of
λ and, generally, a number of oscillatory solutions with
different frequencies ω are possible. Physically, such so-
lutions correspond to formation of several synchronous
enzymic groups. This effect has been previously exten-
sively investigated for similar systems in small spatial
volumes with full diffusional mixing [11] and we shall
not further discuss it here. The most robust uniform
oscillations, which we consider, are characterized by the
frequency ω ≈ 2π/τ and correspond to the single-group
synchronization. As the result of a wave bifurcation (also
known as the Hopf bifurcation with a finite wave number
[17]), the first unstable modes are traveling waves with
a certain wavenumber q0. Figure 3 shows the bifurca-
tion diagram in the parameter plane (τp, β). Note the
presence of a codimension-2 bifurcation point where the
boundaries of the Hopf and the wave bifurcations join.
To investigate nonlinear dynamics of the system, nu-
merical simulations of equations (1) have been performed
0 0.1 0.2 0.3 0.4
oscillations
codimension−2
wave−Hopf bifurcation 
uniform
ripples
pacemakers/waves
higher frequency/
mixed modes
waves
standing−
traveling
standing waves
FIG. 3: Phase diagram (α0 = 1, α1 = 1000, κ = 10, γ = 10,
c = 100, D = 1000). The Hopf bifurcation (solid line) and the
wave bifurcation (dash-dotted line) boundaries are displayed.
Gray lines show instability of the stationary state with respect
to development of uniform oscillations with two (dashed) and
three (dotted) groups in the well-mixed case. Lines separating
parameter domains with different kinds of patterns are hand-
drawn, based on numerical simulations.
[16]. The explicit Euler integration method has been
used; no-flux boundary conditions were applied. Results
of 1D simulations are summarized in Fig. 3 and examples
of typical observed patterns are shown in Fig. 4. Stand-
ing waves (Fig. 4a) develop when the boundary of the
wave bifurcation (dash-dotted curve) is crossed and uni-
form oscillations are observed above the boundary of the
Hopf bifurcation. Near the codimension-2 point, more
complex behavior was found. This included rippled os-
cillations (Fig. 4b), self-organized pacemakers (Fig. 4c)
and modulated traveling waves (Fig. 4d). The observed
patterns are similar to those previously found in reaction-
diffusion systems with the wave bifurcation [18]. In the
right upper corner of the diagram in Fig. 3, higher fre-
quency oscillations with several synchronous groups take
place.
Two-dimensional simulations of reaction-diffusion
equations (1) with time delay have been performed for
selected parameter values. In 2D simulations, sponta-
neously developing concentric waves (target patterns)
and spiral waves have been observed; target patterns
were however unstable and evolved into pairs of rotat-
ing spiral waves (Fig. 2c and Video 3 [15]). Complex
wave regimes, which can be qualitatively characterized
as turbulence of standing waves, have also been observed
(Fig. 2d and Video 4 [15]).
The mean-field approximation is based on neglect-
ing statistical fluctuations in concentrations of reacting
species [11] and, therefore, it should hold in the high
concentration limit. In Fig. 4, two upper panel rows
display spatiotemporal patterns which are observed in
FIG. 4: Spatiotemporal patterns in a 1D system (in each
panel, the vertical axis is time, running down, and the hor-
izontal axis is the coordinate). The upper two rows are
stochastic simulations (σ = 0) with concentrations c = 1 and
c = 10, the bottom row shows mean-field simulations with
c = 100. (a) τp = 0.3, β = 95/c, (b) τp = 0.14, β = 260/c, (c)
τp = 0.22, β = 600/c, and (d) τp = 0.16, β = 300/c. Other
parameters as in Fig. 3; the system size shown is L = 51
ldiff .
stochastic simulations with parameter values correspond-
ing to the respective mean-field simulations. To compare
mean-field simulations with different enzyme densities,
the following property of equations (1) can be used: in-
troducing relative concentrations ñ0 = n0/c, ñ1 = n1/c
and m̃ = m/c, it can be noticed that they obey the same
equations, but with a rescaled coefficient β̃ = βc. Thus,
essentially the same patterns are observed as long as
the parameter combination βc remains constant. In the
stochastic simulations in Fig. 4, the coefficient β has been
increased to compensate for a decrease in the enzyme
concentration. For larger enzyme concentrations, good
agreement between mean-field predictions and stochas-
tic simulations has been found. In the mean-field equa-
tions (1), intramolecular fluctuations are not taken into
account (σ = 0 and therefore each turnover cycle has
the same fixed duration τ). Stochastic simulations have
been, however, also performed when such fluctuations
were present. Synchronization waves could still be found
even at internal noise levels which corresponded to the
mean relative dispersion ξ of turnover times of about 10%
(with ξ =
/τ ≃ (2στ)
Although the emphasis in this Letter is on the phenom-
ena in two-dimensional enzymic arrays, analogous effects
should be expected for three-dimensional systems repre-
senting aqueous enzymic solutions. The linear stability
analysis, yielding Hopf and wave bifurcation boundaries
(see Fig. 3), is valid also for the 3D geometry. We have
performed preliminary stochastic simulations for thin so-
lution layers with high enzyme concentrations and could
observe synchronization patterns similar to those found
for the enzymic arrays.
A product molecule, released by an enzyme, diffuses
in the solution until it either binds, as a regulatory
molecule, to another enzyme or undergoes a decay. Here,
it should be taken into account that a regulatory molecule
can bind to an allosteric enzyme only at a certain bind-
ing site of characteristic radius R. Using the theory of
diffusion-controlled reactions, the average time ttransit
after which a regulatory product would find a binding
site of one of the enzymes can be roughly estimated
[11] as ttransit = 1/cDR, if enzymes are uniformly dis-
tributed inside the reaction volume with concentration
c. Therefore, binding typically occurs within the dis-
tance Lcorr = (Dttransit)
= (cR)
from the point
where a molecule is released. Obviously, it can only
take place if the product molecule has not undergone
decay until that moment, i.e. if γttransit < 1. This
condition puts a restriction on the enzyme concentration
c, which must be higher than the critical concentration
c∗ = γ/DR. Choosing γ = 103 s−1, D = 10−5 cm2s−1
and R = 10−7 cm, the critical enzyme concentration is
c∗ = 1015 cm−3 = 10−6 M. A similar estimate can be
obtained when enzymes are immobilized on a plane im-
mersed into a reactive solution; in this case the mean dis-
tance between the enzymes on the plane should be less
than lc = (Rldiff)
[22]. Although the required en-
zyme concentrations are relatively large, they are within
the range characteristic for biological cells (glycolytic en-
zymes are present [19] in a cell at even higher concentra-
tion of more than 10−5 M). The characteristic temporal
period of developing patterns is determined by the en-
zyme turnover time τ , which typically varies from mil-
liseconds to seconds. The characteristic length scale of
developing wave patterns is determined by the diffusion
length ldiff , which can vary under these conditions from
a fraction of a micrometer to tens of micrometers.
Our analysis shows that spontaneous molecular syn-
chronization of allosteric product-activated enzymes can
be observed in enzymic arrays. Artificial arrays formed
by immobilized protein machines (molecular motors) are
already used in experiments on active nanoscale trans-
port (see [20]). Many enzymes in biological cells are
membrane-bound, thus forming natural enzymic arrays.
Similar phenomena are possible in dense enzyme solu-
tions. In the study by Petty et al. [21], traveling waves
of NAD(P)H and proton concentrations with the wave-
length of about a micrometer were observed inside neu-
trophil cells. These metabolic waves had the temporal
period of about 300 ms, which is by two orders of magni-
tude shorter than the characteristic period of glycolytic
oscillations in the cells and lies closer to the time scales
of turnover cycles of individual enzymes. An intriguing
question, requiring further detailed analysis, is whether
molecular synchronization waves may have already been
seen in these experiments.
Molecular synchronization waves are principally dif-
ferent from classical concentration waves in reaction-
diffusion systems. Under synchronization conditions,
internal conformational states of individual enzyme
molecules in their turnover cycles become strongly cor-
related. In optics, a similar situation is found when a
transition to coherent laser generation has taken place.
Our theoretical analysis may open a way to the investiga-
tions of a new class of spatio-temporal pattern formation
in chemically active molecular systems.
The authors are grateful to M. Falcke and P. Stange for
valuable discussions. Financial support of Japan Society
for the Promotion of Science through a fellowship for
research abroad (Y. T.) is acknowledged.
∗ Present address: Nanobiology Laboratories, Graduate
School of Frontier Biosciences, Osaka University, 1-3 Ya-
madaoka, Suita, Osaka 565-0871, Japan; Electronic ad-
dress: togashi@phys1.med.osaka-u.ac.jp
† Electronic address: mikhailov@fhi-berlin.mpg.de
[1] K. Kinbara, T. Aida, Chem. Rev. 105, 1377 (2005).
[2] L. A. Blumenfeld, A. N. Tikhonov, Biophysical Thermo-
dynamics of Intracellular Processes: Molecular Machines
of the Living Cell (Springer, Berlin 1994).
[3] M. Gerstein, A. M. Lesk, C. Chothia, Biochemistry 33,
6739 (1994).
[4] F. Jülicher, A. Ajdari, J. Prost, Rev. Mod. Phys. 69, 1269
(1997).
[5] H.-Ph. Lerch, A. S. Mikhailov, B. Hess, Proc. Natl. Acad.
Sci. (USA) 99, 15410 (2002).
[6] H.-Ph. Lerch, R. Rigler, A. S. Mikhailov, Proc. Natl.
Acad. Sci. (USA) 102, 10807 (2005).
[7] S. Ramaswamy, J. Toner, J. Prost, Phys. Rev. Lett. 84,
3494 (2000).
[8] P. Lenz, J.-F. Joanny, F. Jülicher, J. Prost, Phys. Rev.
Lett. 91, 108104 (2003).
[9] H.-Y. Chen, Phys. Rev. Lett. 92, 168101 (2004).
[10] A. Goldbeter, Biochemical Oscillations and Cellular
Rhythms (Cambridge University Press, Cambridge 1996).
[11] P. Stange, A. S. Mikhailov, B. Hess, J. Phys. Chem. B
102, 6273 (1998).
[12] P. Stange, A. S. Mikhailov, B. Hess, J. Phys. Chem. B
103, 6111 (1999).
[13] K. Sun, Q. Ouyang, Phys. Rev. E 64, 026111 (2001).
[14] M. Schienbein, H. Gruler, Phys. Rev. E 56, 7116 (1997).
[15] See EPAPS Document No. E-PRLTAO-99-041730
for dynamical evolutions in the 2D simula-
tions. For more information on EPAPS, see
http://www.aip.org/pubservs/epaps.html .
[16] V. Casagrande, Doctoral thesis, Technical University,
Berlin (2006),
http://opus.kobv.de/tuberlin/volltexte/2006/1273/ .
[17] D. Walgraef, Spatio-Temporal Pattern Formation
(Springer, Berlin 1997).
[18] A. M. Zhabotinsky, M. Dolnik, I. R. Epstein, J. Chem.
Phys. 103, 10306 (1995).
[19] B. Hess, A. Boiteux, J. Krüger, Adv. Enzyme Regul. 7,
mailto:togashi@phys1.med.osaka-u.ac.jp
mailto:mikhailov@fhi-berlin.mpg.de
http://www.aip.org/pubservs/epaps.html
http://opus.kobv.de/tuberlin/volltexte/2006/1273/
149 (1969).
[20] H. Hess, G. D. Bachand, Materials Today 8 (12, Suppl.
1), 22 (2005).
[21] H. R. Petty, R. G. Worth, A. L. Kindzelskii, Phys. Rev.
Lett. 84, 2754 (2000).
[22] Diffusion perpendicular to the plane is considered as di-
lution within a layer of effective thickness ≃ ldiff .

  Spatiotemporal pattern formation in a product-activated enzymic reaction at
high enzyme concentrations is investigated. Stochastic simulations show that
catalytic turnover cycles of individual enzymes can become coherent and that
complex wave patterns of molecular synchronization can develop. The analysis
based on the mean-field approximation indicates that the observed patterns
result from the presence of Hopf and wave bifurcations in the considered
system.

Molecular Synchronization Waves in Arrays of Allosterically Regulated Enzymes
Vanessa Casagrande,1 Yuichi Togashi,2, ∗ and Alexander S. Mikhailov2, †
1Hahn-Meitner-Institut, Glienicker Straße 100, 14109 Berlin, Germany
2Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, 14195 Berlin, Germany
Spatiotemporal pattern formation in a product-activated enzymic reaction at high enzyme con-
centrations is investigated. Stochastic simulations show that catalytic turnover cycles of individual
enzymes can become coherent and that complex wave patterns of molecular synchronization can
develop. The analysis based on the mean-field approximation indicates that the observed patterns
result from the presence of Hopf and wave bifurcations in the considered system.
PACS numbers: 82.40.Ck, 87.18.Pj, 82.39.Fk, 05.45.Xt
Molecular machines, such as molecular motors, ion
pumps and some enzymes, play a fundamental role in
biological cells and can be also used in the emerging soft-
matter nanotechnology [1]. A protein machine is a cyclic
device, where each cycle consists of conformational mo-
tions initiated by binding of an energy-bringing ligand
[2, 3]. In motors, such internal motions generate me-
chanical work [4], while in enzymes they enable or facil-
itate chemical reaction events (see, e.g., [5, 6]). Much
attention has been attracted to studies of biomembranes
with ion pumps and molecular motors, where membrane
instabilities and synchronization effects have been ana-
lyzed [7, 8, 9]. Here, a different class of distributed ac-
tive molecular systems — formed by enzymes — is con-
sidered. The catalytic activity of an allosteric enzyme
protein is activated or inhibited by binding of small reg-
ulatory molecules; the role of such regulatory molecules
can be played by products of the same reaction [10]. Pre-
vious investigations of simple product-regulated enzymic
systems [11, 12] and enzymic networks [13] in small spa-
tial volume with full diffusional mixing have shown that
spontaneous synchronization of molecular turnover cycles
can take place there. External molecular synchroniza-
tion of enzymes of the photosensitive P-450 dependent
monooxygenase system by periodic optical forcing has
been experimentally demonstrated [14].
In this Letter, spatiotemporal pattern formation in en-
zymic arrays is investigated. In such systems, immobile
enzymes are attached to a solid planar support immersed
into a solution through which fresh substrate is supplied
and product molecules are continuously removed. Prod-
uct molecules released by an enzyme diffuse through the
solution and activate catalytic turnover cycles of neigh-
bouring enzymes in the array.
A simple stochastic model [12] of an enzyme as a cyclic
machine (a stochastic phase oscillator), shown in Fig. 1,
is used. Binding of a substrate molecule to an enzyme i
initiates an ordered internal conformational motion, de-
scribed by the conformational phase coordinate φi. The
initial state corresponds to the phase φi = 0. The cat-
alytic conversion event takes place and the product is
released at the state φp inside the cycle. After that, the
substrate
enzyme
regulatory
molecule
feedback
product
FIG. 1: (Color online) A sketch of the model.
conformational motion continues until the equilibrium
state of the enzyme (φi = 1) is finally reached. Initi-
ation of a turnover cycle is a random event, occurring
at a certain probability rate. We assume that substrate
is present in abundance, and its concentration is not af-
fected by the reactions. Conformational motion inside
the cycle is modeled as a stochastic diffusional drift pro-
cess, described by equation
φi = v+ ηi(t), where v is the
mean drift velocity and ηi(t) is an internal white noise
with 〈ηi(t)ηj(t
′)〉 = 2σδijδ(t− t
′) where σ specifies inten-
sity of intramolecular fluctuations.
Allosterically activated enzymes possess a site on their
surface where regulatory molecules can become bound.
Binding of a regulatory molecule leads to conformational
change that enhances catalytic activity of the enzyme. A
regulatory molecule binds to an enzyme with rate con-
stant β and dissociate from it with rate constant κ. Bind-
ing of a regulatory molecule at an enzyme raises its prob-
ability to start a cycle from α0 to α1. We assume that
a regulatory molecule can bind to an enzyme only in its
rest state and this molecule is released when the cycle
is started. The role of regulatory molecules is played
by product molecules of the same reaction. Immobile
enzymes are randomly distributed in space with concen-
tration c. Product diffuses at diffusion constant D and
undergoes decay at rate constant γ. The characteristic
diffusion length of product molecules is ldiff =
In our stochastic 2D simulations, the medium was dis-
cretized into spatial cells (up to 256 × 256), each con-
http://arxiv.org/abs/0704.0021v2
FIG. 2: Stochastic (a,b) and mean-field (c,d) simulations of
2D wave patterns; (a) τp = 0.14, c = 1, and β = 300, (b)
τp = 0.25, c = 10, and β = 10, (c) τp = 0.14, c = 1, and
β = 300, (d) τp = 0.34, c = 100, and β = 1.42. Other
parameters are α0 = 1, α1 = 1000, κ = 10, γ = 10, σ = 0,
D = 100. The linear size of the shown area is L = 40 ldiff in
all panels.
taining a number of enzyme molecules. The cells were
so small that diffusional mixing of product molecules in
a cell within the shortest characteristic time of the reac-
tion could always take place. Each enzyme was described
by the stochastic model given above; diffusion of product
molecules was modeled as a random walk over a discrete
cell lattice. The mean cycle time τ = 1/v was chosen as
the time unit (τ = 1). Systems including up to 655 360
enzymes were used in the simulations.
Figure 2a,b (see also Videos 1 and 2 in ref. [15]) shows
two typical examples of stochastic 2D simulations. Here,
spatial distributions of product molecules are displayed.
Waves of product concentration are propagating through
the medium. In a peak of a wave, many locally present
enzymes are simultaneously releasing product molecules.
Since product release can take place only at a certain
stage inside the cycle, this means that the cycles of en-
zymes are locally synchronized. Not only regular wave
structures, such as rotating spiral waves or target pat-
terns (Fig. 2a), but also complex regimes of wave turbu-
lence (Fig. 2b) have been observed.
To understand and interpret stochastic simulation re-
sults, an analytical study of the system in the mean-field
approximation, which holds in the limit of high enzyme
concentrations, has been performed. In this approxima-
tion, the system is characterised by three continuous vari-
ables n0(r, t), n1(r, t) and m(r, t) which represent local
concentrations of enzymes in the rest state without or
with regulatory molecules attached (n0 and n1) and local
concentration of the product (m). For simplicity, internal
fluctuations in enzymes are neglected (σ = 0). Thus, all
enzymes which have started their cycles at some time t
would release their products at a definite time t+τp (with
τp = φp/v) and finish their cycles, returning to the rest
state, at time t + τ . Therefore, the system is described
by a set of three reaction-diffusion equations with time
delays,
= βmn0 − κn1 − α1n1 (1a)
= −βmn0 + κn1 − α0n0 + α0n0(t− τ)
+α1n1(t− τ) (1b)
= −βmn0 + κn1 + α1n1 − γm+ α0n0(t− τp)
+α1n1(t− τp) +D∇
2m. (1c)
The system always has a uniform stationary state
with certain concentrations n0, n1 and m, which can be
found as solutions of the respective algebraic equations.
This state corresponds to the absence of synchronization.
However, it may become unstable if allosteric activation
is strong enough. To analyze stability, small perturba-
tions δn0, δn1 and δm are added to the stationary state,
equations (1) are linearized and their solutions are sought
as δn0 ∼ δn1 ∼ δm ∼ exp (λqt− iqx) with λq = µq+iωq.
Thus, each spatial mode with wavevector q is character-
ized by its frequency ωq and its rate of growth µq. The
properties µq and ωq are given by the roots of a charac-
teristic equation which is determined by the linearization
matrix of equations (1). The steady state becomes unsta-
ble when at least one spatial mode with some wavenum-
ber q0 starts to grow (µq0 > 0).
As the bifurcation parameter, coefficient β can be cho-
sen. If regulatory molecules cannot bind to enzymes
(β = 0), feedback is absent and instabilities are not pos-
sible. On the other hand, allosteric activation becomes
strong if regulatory molecules can easily bind and, in this
case, emergence of oscillations and wave patterns can be
expected. Our bifurcation analysis reveals that, depend-
ing on the parameters of the system, it can exhibit ei-
ther a Hopf or a wave bifurcation [16]. As a result of
the Hopf bifurcation, uniform oscillations with q = 0 de-
velop. Because of the presence of delays in equations (1),
the characteristic equation is nonpolynomial in terms of
λ and, generally, a number of oscillatory solutions with
different frequencies ω are possible. Physically, such so-
lutions correspond to formation of several synchronous
enzymic groups. This effect has been previously exten-
sively investigated for similar systems in small spatial
volumes with full diffusional mixing [11] and we shall
not further discuss it here. The most robust uniform
oscillations, which we consider, are characterized by the
frequency ω ≈ 2π/τ and correspond to the single-group
synchronization. As the result of a wave bifurcation (also
known as the Hopf bifurcation with a finite wave number
[17]), the first unstable modes are traveling waves with
a certain wavenumber q0. Figure 3 shows the bifurca-
tion diagram in the parameter plane (τp, β). Note the
presence of a codimension-2 bifurcation point where the
boundaries of the Hopf and the wave bifurcations join.
To investigate nonlinear dynamics of the system, nu-
merical simulations of equations (1) have been performed
0 0.1 0.2 0.3 0.4
oscillations
codimension−2
wave−Hopf bifurcation 
uniform
ripples
pacemakers/waves
higher frequency/
mixed modes
waves
standing−
traveling
standing waves
FIG. 3: Phase diagram (α0 = 1, α1 = 1000, κ = 10, γ = 10,
c = 100, D = 1000). The Hopf bifurcation (solid line) and the
wave bifurcation (dash-dotted line) boundaries are displayed.
Gray lines show instability of the stationary state with respect
to development of uniform oscillations with two (dashed) and
three (dotted) groups in the well-mixed case. Lines separating
parameter domains with different kinds of patterns are hand-
drawn, based on numerical simulations.
[16]. The explicit Euler integration method has been
used; no-flux boundary conditions were applied. Results
of 1D simulations are summarized in Fig. 3 and examples
of typical observed patterns are shown in Fig. 4. Stand-
ing waves (Fig. 4a) develop when the boundary of the
wave bifurcation (dash-dotted curve) is crossed and uni-
form oscillations are observed above the boundary of the
Hopf bifurcation. Near the codimension-2 point, more
complex behavior was found. This included rippled os-
cillations (Fig. 4b), self-organized pacemakers (Fig. 4c)
and modulated traveling waves (Fig. 4d). The observed
patterns are similar to those previously found in reaction-
diffusion systems with the wave bifurcation [18]. In the
right upper corner of the diagram in Fig. 3, higher fre-
quency oscillations with several synchronous groups take
place.
Two-dimensional simulations of reaction-diffusion
equations (1) with time delay have been performed for
selected parameter values. In 2D simulations, sponta-
neously developing concentric waves (target patterns)
and spiral waves have been observed; target patterns
were however unstable and evolved into pairs of rotat-
ing spiral waves (Fig. 2c and Video 3 [15]). Complex
wave regimes, which can be qualitatively characterized
as turbulence of standing waves, have also been observed
(Fig. 2d and Video 4 [15]).
The mean-field approximation is based on neglect-
ing statistical fluctuations in concentrations of reacting
species [11] and, therefore, it should hold in the high
concentration limit. In Fig. 4, two upper panel rows
display spatiotemporal patterns which are observed in
FIG. 4: Spatiotemporal patterns in a 1D system (in each
panel, the vertical axis is time, running down, and the hor-
izontal axis is the coordinate). The upper two rows are
stochastic simulations (σ = 0) with concentrations c = 1 and
c = 10, the bottom row shows mean-field simulations with
c = 100. (a) τp = 0.3, β = 95/c, (b) τp = 0.14, β = 260/c, (c)
τp = 0.22, β = 600/c, and (d) τp = 0.16, β = 300/c. Other
parameters as in Fig. 3; the system size shown is L = 51
ldiff .
stochastic simulations with parameter values correspond-
ing to the respective mean-field simulations. To compare
mean-field simulations with different enzyme densities,
the following property of equations (1) can be used: in-
troducing relative concentrations ñ0 = n0/c, ñ1 = n1/c
and m̃ = m/c, it can be noticed that they obey the same
equations, but with a rescaled coefficient β̃ = βc. Thus,
essentially the same patterns are observed as long as
the parameter combination βc remains constant. In the
stochastic simulations in Fig. 4, the coefficient β has been
increased to compensate for a decrease in the enzyme
concentration. For larger enzyme concentrations, good
agreement between mean-field predictions and stochas-
tic simulations has been found. In the mean-field equa-
tions (1), intramolecular fluctuations are not taken into
account (σ = 0 and therefore each turnover cycle has
the same fixed duration τ). Stochastic simulations have
been, however, also performed when such fluctuations
were present. Synchronization waves could still be found
even at internal noise levels which corresponded to the
mean relative dispersion ξ of turnover times of about 10%
(with ξ =
/τ ≃ (2στ)
Although the emphasis in this Letter is on the phenom-
ena in two-dimensional enzymic arrays, analogous effects
should be expected for three-dimensional systems repre-
senting aqueous enzymic solutions. The linear stability
analysis, yielding Hopf and wave bifurcation boundaries
(see Fig. 3), is valid also for the 3D geometry. We have
performed preliminary stochastic simulations for thin so-
lution layers with high enzyme concentrations and could
observe synchronization patterns similar to those found
for the enzymic arrays.
A product molecule, released by an enzyme, diffuses
in the solution until it either binds, as a regulatory
molecule, to another enzyme or undergoes a decay. Here,
it should be taken into account that a regulatory molecule
can bind to an allosteric enzyme only at a certain bind-
ing site of characteristic radius R. Using the theory of
diffusion-controlled reactions, the average time ttransit
after which a regulatory product would find a binding
site of one of the enzymes can be roughly estimated
[11] as ttransit = 1/cDR, if enzymes are uniformly dis-
tributed inside the reaction volume with concentration
c. Therefore, binding typically occurs within the dis-
tance Lcorr = (Dttransit)
= (cR)
from the point
where a molecule is released. Obviously, it can only
take place if the product molecule has not undergone
decay until that moment, i.e. if γttransit < 1. This
condition puts a restriction on the enzyme concentration
c, which must be higher than the critical concentration
c∗ = γ/DR. Choosing γ = 103 s−1, D = 10−5 cm2s−1
and R = 10−7 cm, the critical enzyme concentration is
c∗ = 1015 cm−3 = 10−6 M. A similar estimate can be
obtained when enzymes are immobilized on a plane im-
mersed into a reactive solution; in this case the mean dis-
tance between the enzymes on the plane should be less
than lc = (Rldiff)
[22]. Although the required en-
zyme concentrations are relatively large, they are within
the range characteristic for biological cells (glycolytic en-
zymes are present [19] in a cell at even higher concentra-
tion of more than 10−5 M). The characteristic temporal
period of developing patterns is determined by the en-
zyme turnover time τ , which typically varies from mil-
liseconds to seconds. The characteristic length scale of
developing wave patterns is determined by the diffusion
length ldiff , which can vary under these conditions from
a fraction of a micrometer to tens of micrometers.
Our analysis shows that spontaneous molecular syn-
chronization of allosteric product-activated enzymes can
be observed in enzymic arrays. Artificial arrays formed
by immobilized protein machines (molecular motors) are
already used in experiments on active nanoscale trans-
port (see [20]). Many enzymes in biological cells are
membrane-bound, thus forming natural enzymic arrays.
Similar phenomena are possible in dense enzyme solu-
tions. In the study by Petty et al. [21], traveling waves
of NAD(P)H and proton concentrations with the wave-
length of about a micrometer were observed inside neu-
trophil cells. These metabolic waves had the temporal
period of about 300 ms, which is by two orders of magni-
tude shorter than the characteristic period of glycolytic
oscillations in the cells and lies closer to the time scales
of turnover cycles of individual enzymes. An intriguing
question, requiring further detailed analysis, is whether
molecular synchronization waves may have already been
seen in these experiments.
Molecular synchronization waves are principally dif-
ferent from classical concentration waves in reaction-
diffusion systems. Under synchronization conditions,
internal conformational states of individual enzyme
molecules in their turnover cycles become strongly cor-
related. In optics, a similar situation is found when a
transition to coherent laser generation has taken place.
Our theoretical analysis may open a way to the investiga-
tions of a new class of spatio-temporal pattern formation
in chemically active molecular systems.
The authors are grateful to M. Falcke and P. Stange for
valuable discussions. Financial support of Japan Society
for the Promotion of Science through a fellowship for
research abroad (Y. T.) is acknowledged.
∗ Present address: Nanobiology Laboratories, Graduate
School of Frontier Biosciences, Osaka University, 1-3 Ya-
madaoka, Suita, Osaka 565-0871, Japan; Electronic ad-
dress: togashi@phys1.med.osaka-u.ac.jp
† Electronic address: mikhailov@fhi-berlin.mpg.de
[1] K. Kinbara, T. Aida, Chem. Rev. 105, 1377 (2005).
[2] L. A. Blumenfeld, A. N. Tikhonov, Biophysical Thermo-
dynamics of Intracellular Processes: Molecular Machines
of the Living Cell (Springer, Berlin 1994).
[3] M. Gerstein, A. M. Lesk, C. Chothia, Biochemistry 33,
6739 (1994).
[4] F. Jülicher, A. Ajdari, J. Prost, Rev. Mod. Phys. 69, 1269
(1997).
[5] H.-Ph. Lerch, A. S. Mikhailov, B. Hess, Proc. Natl. Acad.
Sci. (USA) 99, 15410 (2002).
[6] H.-Ph. Lerch, R. Rigler, A. S. Mikhailov, Proc. Natl.
Acad. Sci. (USA) 102, 10807 (2005).
[7] S. Ramaswamy, J. Toner, J. Prost, Phys. Rev. Lett. 84,
3494 (2000).
[8] P. Lenz, J.-F. Joanny, F. Jülicher, J. Prost, Phys. Rev.
Lett. 91, 108104 (2003).
[9] H.-Y. Chen, Phys. Rev. Lett. 92, 168101 (2004).
[10] A. Goldbeter, Biochemical Oscillations and Cellular
Rhythms (Cambridge University Press, Cambridge 1996).
[11] P. Stange, A. S. Mikhailov, B. Hess, J. Phys. Chem. B
102, 6273 (1998).
[12] P. Stange, A. S. Mikhailov, B. Hess, J. Phys. Chem. B
103, 6111 (1999).
[13] K. Sun, Q. Ouyang, Phys. Rev. E 64, 026111 (2001).
[14] M. Schienbein, H. Gruler, Phys. Rev. E 56, 7116 (1997).
[15] See EPAPS Document No. E-PRLTAO-99-041730
for dynamical evolutions in the 2D simula-
tions. For more information on EPAPS, see
http://www.aip.org/pubservs/epaps.html .
[16] V. Casagrande, Doctoral thesis, Technical University,
Berlin (2006),
http://opus.kobv.de/tuberlin/volltexte/2006/1273/ .
[17] D. Walgraef, Spatio-Temporal Pattern Formation
(Springer, Berlin 1997).
[18] A. M. Zhabotinsky, M. Dolnik, I. R. Epstein, J. Chem.
Phys. 103, 10306 (1995).
[19] B. Hess, A. Boiteux, J. Krüger, Adv. Enzyme Regul. 7,
mailto:togashi@phys1.med.osaka-u.ac.jp
mailto:mikhailov@fhi-berlin.mpg.de
http://www.aip.org/pubservs/epaps.html
http://opus.kobv.de/tuberlin/volltexte/2006/1273/
149 (1969).
[20] H. Hess, G. D. Bachand, Materials Today 8 (12, Suppl.
1), 22 (2005).
[21] H. R. Petty, R. G. Worth, A. L. Kindzelskii, Phys. Rev.
Lett. 84, 2754 (2000).
[22] Diffusion perpendicular to the plane is considered as di-
lution within a layer of effective thickness ≃ ldiff .
