PHYSICAL REVIEW E 69, 051904 (2004)

Phase transition of traveling waves in bacterial colony pattern
Joe Yuichiro Wakano, Atsushi Komoto, and Yukio Yamaguchi
Department of Chemical System Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-Ku, Tokyo 113-8656, Japan(Received 9 February 2003; published 12 May 2004)
Depending on the growth condition, bacterial colonies can exhibit different morphologies. Many previous
studies have used reaction diffusion equations to reproduce spatial patterns. They have revealed that nonlinear
reaction term can produce diverse patterns as well as nonlinear diffusion coefficient. Typical reaction term
consists of nutrient consumption, bacterial reproduction, and sporulation. Among them, the functional form of
sporulation rate has not been biologically investigated. Here we report experimentally measured sporulation
rate. Then, based on the result, a reaction diffusion model is proposed. One-dimensional simulation showed theexistence of traveling wave solution. We study the wave form as a function of the initial nutrient concentrationand find two distinct types of solution. Moreover, transition between them is very sharp, which is analogous tophase transition. The velocity of traveling wave also shows sharp transition in nonlinear diffusion model,
which is consistent with the previous experimental result. The phenomenon can be explained by separatrix inreaction term dynamics. Results of two-dimensional simulation are also shown and discussed.
DOI: 10.1103/PhysRevE.69.051904

PACS number(s): 87.18.Hf, 47.20.Hw, 87.23.Cc

I. INTRODUCTION

Cooperative self-organization of bacterial colonies has
been intensively studied both experimentally and theoretically [1–8]. The series of these studies should be considered
as a part of the large stream where we pursue the universal
comprehension of pattern formation. One of the most impressive and illustrating way of study is to perform numerical calculation of nonlinear partial differential equations in
multiple-dimensional space to show various beautiful patterns [9–14]. However, they are so complicated and it is very
difficult to understand why a certain pattern is achieved only
for appropriate parameters. Another effective way is to study
a simple element extracted from the complicated system. In
this study we focus on traveling wave solution (TWS) which
is a very simple pattern formation in which a fix-shaped
wave propagates at a constant velocity. TWS often exists
even in the nonlinear reaction diffusion model which produces complicated pattern under appropriate condition. TWS
is clearly one of the most important element of pattern formation and some analytical works have already been done on
the study of bacterial colony pattern formation [15,16].
Cohen et al. reported that colony expanding velocity increased as the initial nutrient concentration increased in theirexperiment [3]. Moreover, they reported the sharp change inthe functional form of velocity as a function of nutrient level.
The function was discontinuous at the point where twodimensional pattern changed. Although the observed twodimensional colony growth was not simple traveling wave,
this study implies that the discontinuous change of TWS
might be the source of the pattern formation. If a travelingwave propagates at a constant velocity everywhere, the result
should be simple disklike (two-dimensional) or spherelike(three-dimensional) pattern. Experimentally observed diverse
spatial patterns imply that the propagating velocity or theform of interface varies depending on the local properties,e.g., curvature or nutrient concentration. If the system has
two distinct TWSs both of which are quasistable but notglobally stable, the spatial pattern in which the two TWSs
1539-3755/2004/69(5)/051904(9)/$22.50

are mixed might evolve. Discontinuous change of TWS is
one of the keys which describes the formation of spatial patterns.
It is experimentally known that some bacterial strains belonging to the genus Bacillus or Paenibacillus produce very
complicated and clear patterns, i.e., branching. One characteristic of the genus is the ability to sporulate. Spores areinactive and dormant state of cells, which can survive starvation or dry environment. Bacteria species with sporulation
ability produce the complex patterns of spore distribution.
Previous theoretical studies revealed that colony pattern isbest understood when we consider the pattern as the historyof active bacteria density. Spores are the history of bacteria
activity and sporulation is an important process in which
active bacteria density is recorded in the history. Many models of bacterial colony pattern formation assume sporulation,
however, there are few experimental studies that measured
the sporulation rate.
Mimura et al. suggested a model with linear diffusion andnonlinear reaction terms [14]. The model has the ability toreproduce four out of five observed patterns of bacterialcolony. Mathematically, the ability comes from the separatrix
in reaction term dynamics (a dynamical system where diffusion terms are neglected). As Mimura et al. describe in theirpaper, the model is constructed so that reaction term dynamics is an excitable system. When the global behavior of the
trajectory discontinuously changes depending on the initialpoint (i.e., separatrix), this may produce very strong nonlinear effect. Biologically, however, separatrix in Mimuramodel is not very realistic. It assumes that sporulation rate isa decreasing function of nutrient and bacteria density. As
sporulation is considered as adaptation to environmental deterioration, the rate is biologically expected to be high when
nutrient density per bacterium is low, if it depends on bacteria density. Anyway, without the experimental data of sporulation rate, it is impossible to build a biologically reasonablemodel.
In this study, we first present the experimentally measured
sporulation rate of bacteria as a function of the initial nutri-

69 051904-1©2004 The American Physical Society

PHYSICAL REVIEW E 69, 051904 (2004)

WAKANO, KOMOTO, AND YAMAGUCHI

FIG. 1. Optical micrographs of B. circulans after 48 h incubation; (a) N = 0 g / L and (b) N = 4.5 g / L, respectively. Each scale bar
indicates 5 ␮m. The arrows indicate a vegetative cell and a spore.

ent. Due to the experimental difficulty, the quantitative relationship between sporulation rate and reproduction rate isunclear. The relationship can be qualitatively classified into
three cases. Based on the experimental result, we propose
three types of models. One of them has separatrix in its reaction term dynamics while the others do not. We perform
one-dimensional computer simulations to study the existenceof TWS and the dependence of the wave form and velocity
on the initial nutrient concentration. Linear diffusion andnonlinear degenerate diffusion are studied. Results of twodimensional simulation are also shown.
II. EXPERIMENTS AND RESULTS

In order to study sporulation rate, Bacillus bacteria wereincubated in the liquid culture medium and the number of
total cells and the ratio of spores were measured (Appendix
A). Total density of both active (vegetative) cells and spores
increased exponentially in the first 5 h (data not shown).
This result is consistent with the previous work that the cell
division of active bacteria is assumed to be a first-order reaction; db / dt = kb. Reaction rate coefficient k is estimatedfrom the data of exponential growth stage. The dependenceof k on the initial nutrient concentration N, well fits a wellknown Michaelis-Menten-type functionk=kmaxCN,
K + CN

where constant parameters are estimated as kmax = 0.6 h−1 andK = 2.2 g / L, respectively.
We investigated the relationship between sporulation andthe nutrient concentration. As sporulation is an adaptation tobad environmental conditions such as nutrient starvation or
drying, sporulation rate has been considered as a decreasing
function of nutrient concentration [14]. However, as shownin Fig. 1, few spores are formed under nutrient-free conditioneven after 48 h incubation. Thus we suppose that bacteria
require some nutrient in order to sporulate. Biologicalknowledge of sporulation process supports the hypothesis[17]. In sporulation, DNA in a bacterial cell is replicated asin regular cell division. Inside the parent cell, one of them isFIG. 2. The ratio of spores to total cells plotted vs the initialnutrient concentration. When nutrient concentration is very small,bacteria do not sporulate. N = 3.0 g / L is best suited for sporulation.

coated with four-layered hard protein. The coated DNA iscalled a “forespore.” After a forespore matures, the parent
cell releases the forespore, which is called an “endospore”
and can survive long periods of time in the bad environment.
In short, as bacteria cannot replicate their DNA undernutrient-free condition, they cannot sporulate.
The ratio of spores to total cells was observed under thevarious initial nutrient concentration (Fig. 2). Predictably,
few spores are formed at high nutrient level. The number ofspores is highest when N = 3.0 g / L. Although the ratio ofspores is not equal to sporulation rate, the result strongly
suggests that sporulation rate as a function of nutrient concentration has a single peak at a certain nutrient level.
Spores are observed about 23 h after the inoculation. By
this duration sporulation rate is estimated to be much smaller
than reproduction rate whose increasing time constant (doubling time) is estimated 1.0– 1.5 h from the analysis of thegrowth curve in early stage. However, as long as the pattern
dynamics of bacterial colony is concerned, spores are defined
as inactive cells which cannot reproduce. As we could distinguish only endospores, it is presumable that time whenvegetative cells transform into forespores is much earlier
than the time of our observation [17]. So we measured theconcentration of ␣-amylase in liquid culture, which is extracted by bacteria only when they sporulate into a forespore[18]. The drastic increase of ␣-amylase concentration (datanot shown) suggests that formation of forespores occurs
about 6 – 8 h after inoculation. The result implies that sporulation rate might be as large as reproduction rate for a certain
interval of nutrient concentration. The result also implies thatsporulation and germination process’s take a long time andthat no spore returned to the vegetative phase in our experiment. However, to evaluate the quantitative relationship between reproduction rate and sporulation rate, further experimental analysis is necessary.

051904-2

III. MODELS AND RESULTS
A. Reaction term dynamics

First, we consider a model without spatial structure,

PHYSICAL REVIEW E 69, 051904 (2004)

PHASE TRANSITION OF TRAVELING WAVES IN…

FIG. 3. Three possible models for the qualitative relationship ofreproduction rate (dotted) and sporulation rate (solid).

db= ␣共n兲b − ␤共n兲b − ␥b,dtdn= ␣共n兲b,
dt

共1兲

ds
= ␤共n兲b,dt
where b , n, and s represent the density of active bacteria(vegetative cells), nutrient concentration, and the density of
spores, respectively. ␣共n兲, ␤共n兲, and ␥ correspond to reproduction rate, sporulation rate, and death rate, respectively.
For simplicity, we assume that these rates do no depend on
bacteria density. As for reproduction rate, Michaelis-Mententype equation,

␣共n兲 = c␣n,n+1

is known to be appropriate for various bacteria species. The
natural death of the bacteria is not commonly observed.
However, as spores are the special dormant state of bacteria,
it is reasonable to assume that active bacteria cannot survive
starvation. For simplicity, we assume the constant death rateof active bacteria, ␥ which is set to very low values compared to sporulation rate. The unit of bacteria concentrationis rescaled so that nutrient conversion factor is 1. Our experiment suggests that sporulation rate is zero at n = 0 , ⬁ andtakes the maximum value at a certain nutrient concentration.
Due to experimental difficulty, the quantitative relationshipbetween sporulation rate and reproduction rate is unclear.
Therefore, a large degree of freedom exists to determine theactual form of ␤共n兲. Qualitatively, there are three possible
relationships 共Fig. 3兲: 共a兲 sporulation is always slower than
reproduction, 共b兲 sporulation is faster than reproduction only
when nutrient is poor, and 共c兲 sporulation is faster than reproduction only for a certain interval of nutrient concentration.
In order to investigate the global behavior of these threemodels, we perform phase plane analysis. Vector field on
共b , n兲 space is qualitatively determined, regardless of actual
function form (Fig. 4). In models 1 and 2, it is obvious thatthe infinitesimal perturbation of the initial state only affects
the terminal state infinitesimally. On the other hand, numerical calculation suggests that separatrix exists in model 3.
When the initial nutrient concentration is lower than a critical value, all bacteria sporulate or die before nutrient is ex-

FIG. 4. The phase plane analysis of the models (see text). Dotted lines represent isoclines. (a) Models 1 and 2 have only one
isocline and their vector fields are qualitatively the same. (b) Model
3 has three isoclines. Numerical calculation suggests the existenceof separatrix around the critical initial value nⴱ. Two orbits startingfrom the initial values which are lower and higher than the criticalvalue are shown.

hausted. When the initial nutrient concentration is higher
than the critical value, almost all nutrient is consumed before
active bacteria disappear. We present the mathematical proofof the existence of the critical value in Appendix B.
We focus on the separatrix of model 3. It is natural toquestion whether the separatrix still exists when spatial diffusion is introduced. Furthermore, such a model may reproduce the observed sharp change of traveling wave velocity.
The two-dimensional spatial pattern is also of great interest.
These motivated us to construct a model with spatial structure based on model 3.
B. Linear diffusion modelFirst, we introduce linear diffusion⳵b= Dⵜ2b + ␣共n兲b − ␤共n兲,⳵t⳵n= ⵜ2n − ␣共n兲b,⳵t

共2兲where D represents diffusion coefficient of active bacteria.
Some experimental knowledge is neglected for simplicity,such as nutrient chemotaxis or the dependence of bacterial
activity on local nutrient concentration. Due to the simplification, D is constant and the model becomes easy to deal
with. Under appropriate rescaling of time and space, diffusion coefficient of nutrient can be chosen as one without lossof generality. The density of spores s is dependent variable
and is not analyzed here. We assume the sporulation rate to
be a normal function:

␣共n兲 = c␣n,n+1

␤共n兲 = c␤ exp兵− ␴共n − n␤兲2 + ␥其.

共3兲

For the sake of simple notification, death rate is included in␤共n兲. Constant parameters 共c␣ , c␤ , n␤ , ␥兲 are chosen so that a

051904-3

PHYSICAL REVIEW E 69, 051904 (2004)

WAKANO, KOMOTO, AND YAMAGUCHI

FIG. 5. Results of the numerical calculation of linear diffusion
model 共D = 1兲. The spatiotemporal patterns of bacteria (upper row)and nutrient (lower row) are shown. Horizontal and vertical axes
correspond to spatial position and time, respectively. Brightness
represents density (concentration) of bacteria (nutrient), scales ofwhich are kept constant. (a) N = 0.485, (b) N = 0.490, and (c) N= 0.495.

model corresponds to model 3 in the preceding section.
In all numerical calculations, we put 共c␣ , c␤ , n␤ , ␥兲
= 共0.3, 0.1, 0.3, 0.01兲.Crank-Nicholson scheme is applied to
calculate diffusion term with zero flux boundary condition.
Initially, nutrient is uniformly distributed at a level N and thesmall amount of active bacteria are put at the origin while no
spore is present.
As a result of numerical simulation of the model in onedimensional space, traveling wave is formed for most initialconditions. However, we found interesting behavior for acertain small parameter region (Fig. 5). When N is smaller
than the critical value Nⴱ, a bacterial wave has a single peak[Fig. 5(a)]. When N is larger than Nⴱ, a new wave appears
behind the original wave [Fig. 5(b)]. A wave in the backpropagates initially slower, but when the distance between
two waves reaches a certain value, a wave in the back begins
to propagate at the same velocity as the fore wave. In other

FIG. 7. Results of numerical calculation of nonlinear diffusion
model 共D = 1兲. Scales are different for each N. (a) N = 0.64, (b) N= 0.65, and (c) N = 0.95.

words, a wave becomes twofold [Fig. 5(c)]. Moreover, at thesame critical value Nⴱ, nutrient concentration left behind thewave sharply changes (from 0.24 to 0.0026). This critical
initial value Nⴱ is estimated as 0.485⬍ Nⴱ ⬍ 0.49. The result
suggests that according to the initial condition this systemhave two phases, a singlefold traveling wave leaving someamount of nutrient behind and a twofold traveling wave leaving little nutrient behind. We conclude that the separatrix inreaction term dynamics is maintained against introduction of
linear diffusion.
Surprisingly, however, propagating velocity does notshow sharp change (Fig. 6). To guess the velocity analytically, we derive the linear spreading velocity of b in thelinearized equation about the unstable state 共b , n兲 = 共0 , N兲
[19],

⳵b
= Dⵜ2b + Sb,⳵t⳵n= ⵜ2n − ␣共N兲b,⳵t

共4兲S = ␣共N兲 − ␤共N兲,in which the first equation is the linearized FisherKolmogoroff equation whose linear spreading velocity c is
c = 冑SD,

FIG. 6. Velocity of traveling wave (circle), velocity of that in thelinearized system (bold line), and nutrient concentration left behind
the wave (triangle) are shown against the initial nutrient concentration. The phase transition of the wave form occurs at the critical
initial nutrient concentration while the wave velocity changes continuously and smoothly, D = 1.

共5兲when the initial condition has a compact support. The numerical calculation shows that the velocity of TWS in the
nonlinear model 关Eq. 共2兲兴 is equal to the linear spreading
velocity 共Fig. 6兲. This means that the wave has a “pulled
front” in the sense that the global nonlinear dynamics of thetraveling wave is governed by the dynamics around the leading edge where the linearized system is a good approximation 关19兴. This might be the reason why the velocity is acontinuous function of the initial nutrient concentration despite the existence of the separatrix in reaction term dynamics.

051904-4

PHYSICAL REVIEW E 69, 051904 (2004)

PHASE TRANSITION OF TRAVELING WAVES IN…
C. Nonlinear diffusion model

Nonlinear diffusion coefficient has been first applied tobacterial colony pattern problem by Kawasaki et al. [10].
Biologically, such diffusion is considered as a result of cooperative cell motion. There are several evidences of cellcell signal exchange and it is natural to consider that bacterial diffusion is not simple linear diffusion of bacteria.
Mathematically, the following Fisher-Kolmogoroff equation
with nonlinear diffusion is well studied:⳵u= ⵱ „D共u兲 ⵱ u… + u共1 − u兲.
⳵t共6兲When D共u兲 = u, a complete analysis has been carried out
关20–22兴 while more general problems are also intensively
studied 关15,23–25兴. When u = 0, the diffusion term vanishes
and the equation degenerates into an ordinary differential
equation. So this type of equation is often referred to as
nonlinear degenerate parabolic equation. In relation to singularity of u = 0, a characteristic solution is known in this system. The TWS u共x , t兲 = U共x − ct兲 = U共z兲 satisfying
U共− ⬁兲 = 1,

U共z兲 = 0∀ z 艌 z ⴱ,U⬘共zⴱ兲 ⫽ 0,zⴱ ⬍ + ⬁is called as the TWS of sharp type. In this solution, u = 0 in
front of the wave front 共z 艌 zⴱ兲 while u ⬎ 0 behind the front
共z ⬍ zⴱ兲. In the linear diffusion equation, on the other hand, u
only approaches zero as z → ⬁. The TWS of sharp type isknown to exist in models with various reaction terms andnonlinear degenerate diffusion term.
We focus on the property of the sharp wave front. The
preceding discussion of linear conjecture of wave velocity
does not hold in such TWS. Therefore, the wave velocity as
a function of the initial nutrient concentration may change
sharply as the wave form does. In order to keep the model as
simple as possible, we introduce nonlinear diffusion coefficient of bacteria, which is a function of only bacteria density.
We will study the following equations:⳵b= D ⵱ 共b ⵱ b兲 + ␣共n兲b − ␤共n兲,⳵t⳵n= ⵜ2n − ␣共n兲b,⳵t

共7兲and perform numerical calculation of the model in onedimensional space. We will use the same ␣共n兲 and ␤共n兲 [Eq.
(3)].
As a result of computer simulation, traveling wave is always formed (Fig. 7). The separatrix of the wave form remains. Nutrient concentration left behind the wave front
sharply changes at the critical value Nⴱ1 ⯝ 0.64. However, awave form remains single wave until N exceeds another critical value Nⴱ2 ⯝ 0.93. Above the critical value, a wave form
becomes twofold. As is discussed, we are interested in theasymptotic behavior of the wave front. Global wave form
and the detail of the wave front of typical twofold waves inlinear and nonlinear diffusion models are shown in Fig. 8.
The wave front of the linear diffusion model is smooth whilethat of the nonlinear diffusion model is sharp, convex, and

FIG. 8. Wave form of bacteria density in the (a) linear and (b)nonlinear diffusion models 共D = 1兲. Enlarged figures show that thefront is smooth in the linear diffusion model and sharp in the nonlinear diffusion model. (a) N = 0.495 and (b) N = 0.95.

intersects b = 0 line. It is confirmed that all TWSs in thelinear diffusion model are of smooth type and that all TWSs
in the nonlinear diffusion model are of sharp type. As is
expected from these facts, wave velocity showed sharpchange at both initial conditions where wave form sharply
changes [Fig. 9(a)]. Notice the remarkable similarity between this and the experimental result of the colony growthvelocity of Paenibacillus dendritiformis bacteria [Fig. 9(b)],although in the experiment the growth velocity of complicated two-dimensional pattern was measured. Another interesting point is that two types of phase transition, which occur
at the same time in the linear diffusion model, occur independently. At the present state, however, we have no analytical explanation why two different critical values appear.
D. Two-dimensional simulation

In order to study the effect of the phase transition of wave
form and velocity in a two-dimensional version of the model,
we performed numerical calculation. Preliminary simulation
showed that the linear diffusion model [Eqs. (2) and (3)]does not show any structured pattern other than simple disklike pattern. So we focus on the nonlinear diffusion model[Eqs. (3) and (7)]. The density of spores s is calculated asds= c␤ exp兵− ␴共n − n␤兲2其.
dt

共8兲Owing to limited computer time, calculation was done in thepositive quadrant divided into 2000⫻ 2000 square lattices.

051904-5

PHYSICAL REVIEW E 69, 051904 (2004)

WAKANO, KOMOTO, AND YAMAGUCHI

FIG. 9. (a) Velocity of traveling wave (circle) and nutrient concentration left behind the wave (triangle) are shown against theinitial nutrient concentration. The phase transition of both the wave
form and the velocity occurs, D = 1. (b) Experimentally observedcolony growth velocity. Reprinted from Ref. [3], Copyright (2004),
with permission from Elsevier.

The response of the pattern to the initial nutrient concentration N is studied (Fig. 10). As a result, four typical patterns are obtained; namely, fine branching, ring, branching,and dense finger. When N is very low, fine branching patternis observed. This pattern is also similar to dense-branchingmorphology (DBM) pattern. In experiments, both finebranching and DBM are observed when the initial nutrient islow, which is consistent with the present numerical result[3,26].
When N exceeds the first critical value Nⴱ1 ⯝ 0.66, thewhole pattern changes very sharply. Nutrient concentration
left behind changes from 0.35 to 0.002 and the pattern itself
also changes to ring pattern. The corresponding transition in
one-dimensional simulation has been already described in
the preceding section. The result means that the separatrix
originating from the reaction term still remains in twodimensional model with nonlinear diffusion. In the ring pattern, no spore is left behind the bacterial front. There is noexperimental observation which supports this. However, thephase transition of two-dimensional pattern is mathematically interesting.
When N exceeds the second critical value Nⴱ2 ⯝ 0.78, theFIG. 10. Results of two-dimensional simulation of nonlinear
diffusion model 共D = 1兲. Spatial distribution of b + s is shown. The
brightness of b is expressed five times stronger than that of s. Except N = 0.7, a branching pattern is formed. A small amount of nutrient is left behind the front when N = 0.65.

pattern changes into branching pattern. This change is also
very sharp, suggesting that it is another separatrix. It seems
that there are two independent traveling waves of bacteria.
Only forward wave produces branching pattern and back051904-6

PHYSICAL REVIEW E 69, 051904 (2004)

PHASE TRANSITION OF TRAVELING WAVES IN…

FIG. 11. Enlarged figure of two-dimensional branching patternwhen N = 0.9.

wave, propagating at slower speed, contributes little to theglobal pattern.
When N exceeds the third critical value Nⴱ3 ⯝ 0.89, the twobacterial waves are combined to form the twofold wave,
while the global pattern does not change very much. Enlarged pattern is shown in Fig. 11. Singlefold and twofoldwaves are mixed along the front line. Twofold wave in bacterial pattern corresponds to the tip of the growing branch inspore pattern. Twofold wave consists of a fore wave corresponding to the first positive region of net bacterial reproduction rate where nutrient is high 共n ⬎ n3兲 and a back wave
corresponding to the second region where nutrient is low
共n1 ⬍ n ⬍ n2兲 [Fig. 4(b)]. Spores are produced in the regionn2 ⬍ n ⬍ n3 while bacteria simply die out in the region n ⬍ n1.
Therefore, most spores are left behind the fore wave, which
only the twofold wave has.
When N is increased further, the pattern graduallychanges into dense fingering. In experiments, dense fingering
pattern is observed when the initial nutrient is high, too. In
our model, the change is caused by increase in twofold wave
region in the front line.
There is no unique method to measure the growth velocity
of such complex patterns. When we roughly estimate thevelocity as a function of the initial nutrient, it is similar to theresult of one-dimensional model. The velocity is very slow
when N ⬍ Nⴱ1, almost the same when N1 ⬍ N ⬍ Nⴱ2, and thenincreases rapidly when N ⬎ Nⴱ3.
Generally, in branching pattern formation, active spot exists only at the tip of growing branch. The possible mechanism for such spatiotemporal pattern formation is as follows.
First, simple disklike pattern (radial traveling wave solution)
becomes unstable due to some nonlinearity in the model.
Such interface instability has been well studied, e.g.,
Mullins-Sekerka dispersion relation in crystal growth theory.
In the particular field of bacterial colony growth, the interface instability of a simple nonlinear model is intensively
studied [16,27]. Dispersion relation, however, only provides
the most unstable wave number and its stability. In order toform clear pattern, active spot should be independent distinct

aggregate. It has been unclear why wave front evolves toactive and inactive parts which are so clearly distinguishable.
Based on the existence of separatrix, we here propose ahypothesis why clear patterns are observed in our model.
From our calculation on the one-dimensional space, it isnaturally expected that the model potentially has two typesof stable states, i.e., a singlefold TWS and a twofold TWS.
The initial condition determines which is realized. Transition
between these states is discontinuous, i.e., separatrix. Whenthe model is extended to the two-dimensional version, global
TWS may become unstable under some condition. In suchcase, spatial heterogeneity may evolve into the mixture of
singlefold and twofold waves, both of which are locallystable in so small region that the system is approximately
equal to the one-dimensional model in radial direction. The
length of the region along the front line might be determinedby dispersion relation analysis. Separatrix originating fromthe reaction term produces two distinct areas in the wave
front, which moves outward to draw clear branching pattern.
IV. CONCLUSION

We performed an experiment to determine the sporulation
rate as a function of the nutrient concentration. Based on the
result, a reaction diffusion model of bacterial colony growthis proposed. Especially, we focus on separatrix in reaction
term dynamics. We study the effect of separatrix on travelingwave solution to find a different behavior; transition between
a singlefold and a twofold wave. The phase transition of thevelocity of the traveling wave shows great similarity with theprevious experiment. Results of a two-dimensional simulation of our model are generally consistent with the observed
patterns.
In our model with N ⬎ Nⴱ2, the spatial pattern of sporesshows clear branching while that of bacteria does not. This
seems inconsistent with the observed patterns, which should
be studied in future. The effect of the diffusion coefficient of
nutrient, the elimination of the effect of square lattice innumerical calculation, and the mathematical proof or the explanation of separatrix in reaction diffusion equation are also
open problems.
As we mentioned at the beginning of the paper, the studyof bacterial colony growth is not only a special biological
topic but also it must be an attempt to explore universal lawof pattern formation. Good theoretical study suggests new
biological experiments while a new biological result inspires
new models. Although bacteria are one of the simplest formof life, they produce so complex patterns. We hope thepresent study contribute to our comprehension of selforganized pattern formation.
APPENDIX A: EXPERIMENTAL METHOD

There are some Bacillus species each of which produces
the original colony patterns. We used Bacillus circulans (B.
circulans) as a spore forming species, but the qualitative features of sporulation learned in one species can be applied toanother. B. circulans was incubated in the liquid culture medium which was prepared by the following methods [28].

051904-7

PHYSICAL REVIEW E 69, 051904 (2004)

WAKANO, KOMOTO, AND YAMAGUCHI

10 g / L sodium chloride and nutrients with various concentration were dissolved in distilled water. We used Bacto
Tryptone Peptone (TP) and Bacto Yeast Extract (YE) as nutrients, which contain vitamin, amino acid, minerals, and carbon sources. They were purchased from Difco, Detroit, M1,
USA. The weight ratio was kept TP/ YE= 2 and the initial
total concentration N was varied from 0 g / L to 15 g / L. Thissolution was autoclaved for 15 min at 120° C and 2 atm andthen cooled down to room temperature. Bacteria, which had
been stored at −20° C to prevent them from mutation, werepreincubated in this liquid culture with N = 15 g / L. After
15– 18 h preincubation, bacteria density became about 108− 109 number/mL while no spore was formed. The bacteria
density was adjusted to about 107 number/mL and then constant quantity of the solution was incubated in a plastic tube
filled with fresh liquid culture. The tube was let in the bathincubator whose temperature was controlled at 37° C andshaken during incubation. We adopted batch-type incubation,in which fresh nutrient was not fed during incubation.
We took 5 ␮L of test portion from each liquid culture
after 30 h incubation. The total bacteria number in each test
portion was counted by using Thoma blood cell counter andan optical microscope 共⫻400兲 and then the bacteria density
was calculated. The ratio of spores to total cells was analyzed by the following method. Bacteria in every test portion
were stained by Möller’s spore staining method for the easy
optical microscopic observation. Optical micrographs
共⫻1000兲 of every sample were taken by a digital camera. In
the micrographs, spores and vegetative cells are easily distinguishable. The ratio of spores were calculated by imaging
software.

␭ = 0,S

共B4兲

S ⬅ ␣共n兲 − ␤共n兲,

共B5兲

and

with eigenvectors

冉冊
0
1

and

冉 冊

S
,
− ␣共n兲

respectively. The zero eigenvalue means that the equilibrium
points are neutrally stable against perturbation in a n axis
direction, which is trivial because whole n axis is equilibrium. The sign of S determines whether an equilibrium is arepeller or an attractor.
Here we assume the qualitative relationship between ␣共n兲and ␤共n兲 is the same as in model 3 in text. According to thestability, n axis is divided into four subsets; stable line segment S1共0 ⬍ n ⬍ n1兲, unstable line segment U1共n1 ⬍ n ⬍ n2兲,
stable line segment S2共n2 ⬍ n ⬍ n3兲, and unstable half line
U2共n3 ⬍ n兲 [Fig. 4(b)]. As n is monotonically decreasing, any
orbit starting from U1 must enter S1. Let f共p兲 as a terminal
point of an orbit starting from a point p. According to f共p兲,U2 consists of the following two subsets:

APPENDIX B: EXISTENCE OF SEPARATRIXU21 ⬅ 兵p 苸 U2:f共p兲 苸 S1其,共B6兲

U22 ⬅ 兵p 苸 U2:f共p兲 苸 S2其.

共B7兲We consider a dynamical system
As we proved, any orbit converges to the equilibrium anddb= ␣共n兲b − ␤共n兲b,
dt
共B1兲
dn
= − ␣共n兲b,dt
defined in the positive quadrant R+2 = 兵共b , n兲 : b 艌 0 , n 艌 0其. We
assume ␣共0兲 = 0 which guarantees any orbit remains in R+2.
The set of equilibrium points of Eq. (B1) is n axis 共b = 0兲. We
assume continuous functions ␣共n兲 and ␤共n兲 satisfy
0 ⬍ ␣共n兲 ⬍ A ∀ n ⬎ 0,

共B2兲

0 ⬍ ␤1 ⬍ ␤共n兲 ⬍ ␤2 ,共B3兲and thus n is monotonically decreasing function of t. Therefore, chaos or limit cycle never occurs. Apparently, n cannot
diverge so we focus on whether b diverges or not. As long asb ⬎ 0, n decreases. When n = 0, db / dt ⬍ −␤1b holds, whichmeans b converges to zero. Therefore, all orbits converge tothe equilibrium points on n axis.
Eigenvalues of linear analysis around an equilibrium
point 共0 , n兲 on n axis are

U2 = U21 艛 U22 ,共B8兲

U21 艚 U22 = ␾ ,共B9兲hold. Hereafter, we prove the following theorem.
Theorem 1. (existence of separatrix). There exists uniqueNⴱ ⬎ n3 such thatU21 傶 兵共0,n兲 苸 U2:n ⬎ Nⴱ其,U22 傶 兵共0,n兲 苸 U2:n3 ⬍ n ⬍ Nⴱ其.
Proof. It is geometrically apparent that both U21 and U22 areconnected because phase plane is two dimensional andany two orbits never intersect each other. Apparently, U21
must be above U22 on n axis; otherwise there should beintersecting orbits. The theorem is proved if U21⫽ ␾ and U22 ⫽ ␾ are proved.

051904-8

PHYSICAL REVIEW E 69, 051904 (2004)

PHASE TRANSITION OF TRAVELING WAVES IN…

is bounded for n ⬎ n2 because ␣ and ␤ are bounded andbecause ␣ ⬎ 0 for n ⬎ n2. Thus, b can change infinitesimallywith infinitesimal change of n. Moreover, db / dn ⫽ 0 for n⫽ n3. This means that 兩f共p兲 − p3兩 ⬍ O共⑀兲 if 兩p − p3兩 ⬍ ⑀, wherep3 ⬅ 共0 , n3兲. In other words, an orbit starting from a point
infinitesimally above p3 must enter a point infinitesimallybelow p3, which is an element of S2. Thus, U22 ⫽ ␾.
When p → 共0 , ⬁兲, the orbit starting from p intersects n= n3 line at 共⬁ , n3兲 as db / dn ⬎ 0 for n ⬎ n3. However, decrease in b is finite in finite interval 关n2 , n3兴 as inclination of
the orbit is finite. Thus, the orbit starting from a point withvery large n must enter S1, which means U21 ⫽ ␾.䊐

[1] E. O. Budrene and H. C. Berg, Nature (London) 349, 630(1991).
[2] E. O. Budrene and H. C. Berg, Nature (London) 376, 49(1995).
[3] I. Cohen, A. Czirok, and E. Ben-Jacob, Physica A 233, 678(1996).
[4] I. Cohen, I. G. Ron, and E. Ben-Jacob, Physica A 286, 321(2000).
[5] J. Y. Wakano, S. Maenosono, A. Komoto, N. Eiha, and Y.
Yamaguchi, Phys. Rev. Lett. 90, 258102 (2003).
[6] N. H. Mendelson and J. Lega, J. Bacteriol. 180, 3285 (1998).
[7] N. H. Mendelson, A. Bourque, K. Wilkening, K. R. Anderson,
and J. C. Watkins, J. Bacteriol. 181, 600 (1999).
[8] A. Komoto, K. Hanaki, S. Maenosono, J. Y. Wakano, Y.
Yamaguchi, and K. Yamamoto, J. Theor. Biol. 225, 91 (2003).
[9] E. Ben-Jacob, O. Schochet, A. Tenenbaum, I. Cohen, A.
Czirok, and T. Vicsek, Nature (London) 368, 46 (1994).
[10] K. Kawasaki, A. Mochizuki, M. Matsushita, T. Umeda, and N.
Shigesada, J. Theor. Biol. 188, 177 (1997).
[11] M. Matsushita, J. Wakita, H. Itoh, I. Rafols, T. Matsuyama, H.
Sakaguchi, and M. Mimura, Physica A 249, 517 (1998).
[12] I. Golding, Y. Kozlovsky, I. Cohen, and E. Ben-Jacob, Physica
A 260, 510 (1998).
[13] Y. Kozlovsky, I. Cohen, I. Golding, and E. Ben-Jacob, Phys.
Rev. E 59, 7025 (1999).

[14] M. Mimura, H. Sakaguchi, and M. Matsushita, Physica A 282,
283 (2000).
[15] R. A. Satnoianu, P. K. Maini, F. S. Garduno, and J. P. Armitage, Discrete Contin. Dyn. Syst., Ser. B 1, 339 (2001).
[16] J. Muller and W. van Saarloos, Phys. Rev. E 65, 061111(2002).
[17] C. Woese, Microbiology (Saunders, Orland, 1997).
[18] W. L. Nicholson and G. H. Chambliss, J. Bacteriol. 161, 875(1985).
[19] W. van Sarloos, Phys. Rep. 386, 29 (2003).
[20] D. G. Aronson, Density-Dependent Iinteraction-Systems (Academic Press, New York, 1980).
[21] W. I. Newman, J. Theor. Biol. 85, 325 (1980).
[22] J. D. Murray, Mathematical Biology (Springer, New York,1989).
[23] W. I. Newman, J. Theor. Biol. 104, 473 (1983).
[24] F. Sanchez-Garduno and P. K. Maini, J. Math. Biol. 33, 163(1994).
[25] Y. Hosono, Jpn. J. Ind. Appl. Math. 3, 163 (1986).
[26] U. W. Wakita, K. Komatsu, A. Nakahara, T. Matsuyama, and
M. Matsushita, J. Phys. Soc. Jpn. 63, 1205 (1994).
[27] S. Kitsunezaki, J. Phys. Soc. Jpn. 66, 1544 (1997).
[28] N. Eiha, A. Komoto, S. Maenosono, J. Y. Wakano, K. Yamamoto, and Y. Yamaguchi, Physica A 313, 609 (2002).

Inclination of an orbit
db ␣共n兲 − ␤共n兲
=
dn
␣共n兲

051904-9