PHYSICAL REVIEW E 84, 011927 (2011)

Growth states of catalytic reaction networks exhibiting energy metabolism
Yohei Kondo and Kunihiko Kaneko
Department of Basic Science, University of Tokyo, Komaba, Meguro-ku, Tokyo 153-8902, Japan(Received 31 March 2011; published 28 July 2011)
All cells derive nutrition by absorbing some chemical and energy resources from the environment; theseresources are used by the cells to reproduce the chemicals within them, which in turn leads to an increase in theirvolume. In this study we introduce a protocell model exhibiting catalytic reaction dynamics, energy metabolism,
and cell growth. Results of extensive simulations of this model show the existence of four phases with regard to
the rates of both the influx of resources and cell growth. These phases include an active phase with high influxand high growth rates, an inefficient phase with high influx but low growth rates, a quasistatic phase with low
influx and low growth rates, and a death phase with negative growth rate. A mean field model well explains thetransition among these phases as bifurcations. The statistical distribution of the active phase is characterized by apower law, and that of the inefficient phase is characterized by a nearly equilibrium distribution. We also discuss
the relevance of the results of this study to distinct states in the existing cells.
DOI: 10.1103/PhysRevE.84.011927

PACS number(s): 87.17.Aa, 89.75.Fb, 82.39.Rt

I. INTRODUCTION

To unveil the minimal requirement of life, we need to firstunderstand how an ensemble of a variety of catalytic reactions
generates a stable cell that grows to reproduce itself. A cell in
any form—right from its prototype form at its inception to its
currently evolved form—receives an influx of some nutrientchemicals; the cell then effectively uses the energy obtained
from the decomposition of these nutrient chemicals to develop
useful catalysts and membranes. It should be noted that thechemicals and energy obtained from the nutrients throughcatalytic reactions are used for cell growth and cell replication.
Given the fact that cell growth requires a continuous flowof energy and nutrients, it can be said that cell growth will not
be sustained if the cell is in thermal equilibrium. Indeed, therelevance of nonequilibrium open systems to life processes has
been discussed over several decades [1,2]. On the other hand,in a Carnot cycle, the maximal thermodynamic efficiency isachieved under an adiabatic condition, where a weak influx
flow is generated under a condition close to the equilibrium
condition. In this case a cellular state far from equilibrium
may not be appropriate for efficient thermodynamic processes.
Therefore it is pertinent to determine the condition under whichhighly efficient cellular growth will be achieved [3].
This problem regarding efficiency in cell-volume growth isa general problem, so it is essential to understand the primitive
form of a cell, or a protocell. In this case, one needs tounderstand a universal property of a growing cell following
catalytic chemical reaction processes upon external nutrients.
On the other hand, this issue is also important to the present
cells. In this case, it is likely that the rate of conversion ofchemical resources to increase the cell volume depends on thecurrent cellular state, or one may distinguish distinct cellularstates on the basis of the efficiency of conversion.
With regard to such cellular states, it is a known factthat a cell undergoes several “phases” depending on culture
conditions, where each phase is characterized by a distinct
growth rate [4–6]. The appearance of these phases depends
on the nutrient condition; moreover, each phase represents
a distinct cellular state whose chemical composition differsfrom that of any other state. For example, in the case ofbacteria, a high growth rate is achieved in the log phase, which
1539-3755/2011/84(1)/011927(8)is characterized by exponential cell growth under a nutrientrich condition; having said that, bacteria also undergo somesecondary phases, such as a stationary phase characterized
by suppressed growth and a dormant phase characterized
by suppressed metabolism. These observations imply theexistence of qualitatively distinguishable cellular states. With
this background, is it now possible to understand the existenceof such distinct phases, characterized by distinct growth andinflux rates, in a cell? The purpose of the present study is tointroduce a simple “protocell” model having distinct cellularstates and analyze efficiency for cellular growth depending oneach state.
To investigate the origin of multiple cellular states andunderlying energetics theoretically, a cell model including
metabolism and growth is needed. There are extensive studies
conducted on protocell models that are based on catalytic
reactions [7–17] to determine possible conditions that facilitate
the reproduction of cells. Previously, Furusawa and one of
the authors (K.K.) studied a protocell model consisting of alarge number of chemicals that mutually catalyze each other,in order to understand the manner in which a cellular stateis realized [12–15]. Random catalytic networks were used to
identify the condition for steady cell growth, where we founda universal statistical law on the abundance distribution ofchemicals. Interestingly, in spite of the simplicity of the modelwith just catalytic reaction networks, this statistical law was
also confirmed over most of the existing cells [12].
However, in the model, energetics was not discussed. It was
assumed that the energy required for catalytic reactions was
supplied automatically, and also that the growth rate of thecell was equal to the inflow rate of the resources. Furthermore,
multiple distinct cellular states were not observed with regard
to this model. Therefore, the growth rate was uniquelydetermined by the nutrition condition, and the efficiency ofcell growth depending on cellular states could not be studied
at all.
To study the cell-state dependence of growth efficiency, we
extend here the above-mentioned protocell model to include
energy metabolism. For this purpose, we first need to consider
reversible reactions, with the reaction rate constants given by
the chemical potentials of molecules to be consistent with

011927-1©2011 American Physical Society

YOHEI KONDO AND KUNIHIKO KANEKOPHYSICAL REVIEW E 84, 011927 (2011)

thermodynamics [18,19]. This consideration is in contrast tothat pertaining to most models of catalytic reaction networks[7,11–15,20–22], where only irreversible(unidirectional) reactions are considered. Second, we introduce specific chemicals
for “energy currency molecules” such as adenosine triphosphate (ATP) to derive energy from the reactions or to supply
energy to the reactions. This energy currency component isassumed to be used for cell growth and is synthesized by othercatalytic reactions. In fact, the cell growth usually requires
some high-energy chemicals, the synthesis of which requires
energy currency molecules. For example, the synthesis of DNA
requires ATP consumption.
In the case of our model, first we consider that there isan influx of nutrient resources. These nutrients are used in
the synthesis of a variety of catalysts through reactions given
by catalytic reaction networks. An energy currency molecule(e.g., ATP) is synthesized through such reactions, while the cellvolume increases by consuming the energy generated by thismolecule. Further, the metabolic efficiency of this model, i.e.,
the rate of cell growth for consumption of nutrients, depends
on the internal cellular state, i.e., the concentrations of the catalysts. Hence we can study the relationship between the influx
rate and the cell growth rate, depending on the cellular state,
and examine the possibility of existence of distinct cellularstates with unique metabolic efficiencies and cell growth rates.
The present paper is organized as follows. In the next
section we introduce a model consisting of a catalytic reaction
network with energetics, and cellular growth. In Sec. III,
direct simulation of the model and the mean field analysis arepresented, where we report the existence of different phaseswith distinct cellular growth speed and efficiency, as well
as different statistical characteristics of chemical abundances.
Section IV is devoted to summary and discussion.
II. MODEL

Here we consider a very primitive “minimal” protocell
compartmentalized by a membrane, within which encapsulated catalytic reactions progress (see Fig. 1 for schematic
representation). This catalytic reaction is necessary for thesynthesis of all the chemicals within the cell and also brings
about cell growth. For this synthesis, some supply of matterand energy is generally required. For the former, we consider
some nutrient chemical species, from which catalysts aresynthesized. For the latter, a common molecular species thatworks as an “energy currency,” say ATP, is assumed; it isrequired for many endothermal reactions to progress. There
are influxes and effluxes of these and other chemicals. As thenumber of molecules increases, the cell volume increases (as aresult of synthesis of membrane molecules) and, accordingly,
the concentrations of all chemicals within the cell decrease.
For system variables we consider the following: (i) A
set of concentrations of catalytic components whose concentrations are given by (x0 ,x1 , . . . xN−1 ). These chemicals aresynthesized by catalytic reaction(s), and other metabolites are
assumed to exist constantly; hence they appear as constantparameters in the reaction. (ii) An energy currency molecule.
With ATP and adenosine diphosphate (ADP) in mind, we
assume that the conversion of A to B supplies energy (Acorresponds to ATP, and B to ADP). The concentration

j

i

k

i

k

A

B

material

membrane

FIG. 1. (Color online) Schematic representation of our “protocell” model. It consists of an encapsulated chemical reaction network
within a membrane. Green circles indicate catalytic components,
and i,j,k are the component indices; purple circles indicate energy
currency molecules. Arrows across the membrane represent material
exchange by diffusion through the membrane, whereas arrows within
the membrane represent chemical reactions. Conversion from energy
currency component A to B is used for the synthesis of membrane
molecules, resulting in an increase in the cell volume.

variables of these are (yA ,yB ). (iii) The volume of the protocell
is set at V and increases or decreases with time.
Catalytic reaction coupled to the conversion of energy
currency molecules: We consider a reaction process from ito j catalyzed by l, coupled with the conversion of energy
currency molecules. This reaction is given byi + l + A(B) ↔ j + l + B(A),(1)

where i,j,l are reactant, product, and catalyst, respectively. By denoting the standard chemical potentials forthe catalytic components and energy currency components
as (μ0i ,μ0j ,μ0A ,μ0B ), the reaction rate constant satisfying thedetailed balance condition is determined as


0
0
0
0
kij = Aij min 1,e−β(μj +μB(A) −μi −μA(B) ) ,(2)
where β is the inverse temperature and Aij = Aj i [19]. The
reaction rate for i to j is then represented as kij xi xl yA(B) . Notethat the reaction progresses to reduce the free energy of theassociated chemicals, not the energy. From this point of view,
the assumptions of the mass-action law and the rate constants
correspond to taking the form of the chemical potential asμi ≡ μ0i + (1/β) log xi . In fact, the net reaction flux from ito j , i.e., xl (kij xi yA(B) − kj i xj yB(A) ), always reduces the freeenergy per unit volume g = xi μi + xj μj + yA μA + yB μB .
We consider a variety of reactions. The whole set ofreactions constitutes a network given by the tensors CA (i,j,l)and CB (i,j,l), which represent the structure of a catalytic
reaction network, defined by
CA(B) (i,j,l)
1 [If reaction i + l + A(B) → j + l + B(A) exists]=
0 [Otherwise](3)
The nutrient chemical is assumed to diffuse through themembrane and its external concentration is set at X0 ; hence
the efflux rate is given by D0 (X0 − x0 ), with the diffusion
coefficient D0 . For simplicity, other components except thechemical N − 1 are assumed to be impermeable, i.e., Di = 0.
The external concentration XN−1 is set much smaller than X0 ;hence the cell takes 0 as nutrient inflow and diffuses out N − 1
as waste. However, this condition is not essential and can be

011927-2

GROWTH STATES OF CATALYTIC REACTION NETWORKS . . .

PHYSICAL REVIEW E 84, 011927 (2011)

relaxed. Several other Di s could be set at nonzero, and theresult would be qualitatively unchanged.
Cell growth progresses as a result of membrane increases(and the influx of water in proportion). Instead of explicitly
introducing a membrane molecule, here we simply note that
membrane growth needs some energy currency consumption.
Assuming that the precursor molecules for the membrane areabundant in the medium, we consider a phenomenological“growth reaction” that progresses accompanied by the reaction
A → B. We assume that the rate of the growth reaction alsosatisfies the detailed balance condition. Then, defining P and
μ0P as the concentration and standard chemical potential,
respectively, of the membrane precursor and μM as the(effective) chemical potential of the membrane, we write therate of the growth reaction as

0
0
0 
k P yA eβ(μP +μA ) − yB eβ(μM +μB ) ,(4)

First, all of the coefficients of the chemical reaction rate
(Aij s) are set at 1 for simplicity. Second, the standard chemical
potentials of the catalysts (μ0i s) are sampled from a uniform
distribution from 0 to 1. Third, if the component 0 is anutrient, the energy currency A is assumed to be synthesized when the component 0 is degraded by appropriately
setting the directions between the catalytic reactions andthe conversion of the energy currency component. For allother reactions, the coupling directions have been determined
randomly.
Finally, we make a cautious comment on the effect of theinflux of the energy currency molecules on the intracellularenergy currency composition. The compositions of the energy
currency molecules A and B outside the cell are expected to
be in chemical equilibrium. Then, when the energy currency
molecules in such a composition flow in with cell growth, thecompositions of the intracellular energy currency components(yA ,yB ) will be equilibrated. We can treat this effect byeqeq(yA − yA ) V1 dV
onto the time evolution of yA , where yA is
dt
the concentration of yA at equilibrium, and we carried out
simulations for such a model. However, we confirmed thatthe results discussed in the present paper are qualitatively thesame, independent of this effect; hence we neglect this effectfor simplicity.

where k is a rate constant. In addition, we assume that theconcentration of the precursor is fixed. Hence, the volume
growth is simply given by
1 dV= γ kg (yA − yB eβ )
V dt

(5)

where  [= μM + μ0B − μ0P − (1/β) log P − μ0A ] is the en00ergy required for the process, kg (=kP eβ(μP +μA ) ) is its rate
constant, and γ is the conversion rate from the membrane
molecules to the cell volume.
Here, P ,μ0P ,μM ,μ0A ,μ0B ,k are direct independent parameters controlling the reactions, and kg , are functions of theseparameters. However, the behaviors of the model as dynamicalsystems are well represented by choosing kg , (and X0 ) as
relevant control parameters, so we use these as the relevant
parameters in the present model.
As a result of the growth in cell volume V , given by Eq. (5),all the chemicals are diluted accordingly. This dilution leads
to a decrease in the concentrations of all the chemicals in
proportion to themselves, i.e., xi V1 dV. As a result of thisdt
dilution, it is expected that the total concentration of theenergy currency components (yA + yB ) would decrease and,without a supply, the system would cease to grow. In reality,
such molecules would also be diffused from the outside,
say phosphorate molecules, the material components for ATP
and ADP. Hence the energy currency molecules are assumed
to flow in with increasing the cell volume so that yA + yB
is constant. Consequently, the two concentrations are notindependent. We therefore need to consider only the time
evolution of yA .
Taken together, the time evolution of the concentrations of
the catalysts xi s and energy currency yA are given by



dxi
Cn (j,i,l)xl (kj i xj yn − kij xi ym )=dt(n,m)=(A,B),(B,A) i,j,l+ Di (Xi − xi ) − xi

1 dV
V dt

dyA
= −kg (yA − yB eβ )dt+CB (i,j,l)xl (kij xi yB − kj i xj yA ).
i,j,l

(6)

III. RESULTS
A. Emergence of states with different metabolic efficiencies

We simulated the present model using several parametervalues and taking a variety of catalytic reaction networks.
Each network was generated randomly with a given total
reaction number G under the constraint that it is percolated.
From these simulations we found common behaviors in thepresent protocell model, as long as the reaction path number
was sufficient. Indeed, as long as the average path number per
chemical, i.e., K = G/N, satisfies K > K ∗ ∼ 8, the common
behavior is observed, independent of the network. If K issmaller than that, the behavior often depends on the network,
which we do not discuss here.
For most of these cases, the chemical concentrations andgrowth rates converge to a fixed point after some relaxation
time. If the growth rate is positive, a fixed point represents
a cellular state with steady growth. If it is negative, the cell
shrinks and will ultimately collapse. We are interested in theformer case here. Since the approach to fixed points is quite
common, we focus on the properties of such fixed points,
without going into details of the transient relaxation process.
For most parameter values the stable fixed point is unique, but
at some parameter values there are two stable fixed points,
leading to bistability in the cellular state; this is discussed
below.
In cell state characterization, the growth rate λ ≡ V1 dV, anddt
the influx Jin ≡ D0 (X0 − x0 ) are important order parameters
for characterizing the protocell; the “metabolic efficiency” η,defined as η ≡ λ/Jin , is another basic order parameter.
With regard to the internal properties of a protocell, thedensity m is defined as the total concentration of catalysts, asm ≡ i xi . The density is an important order parameter andis related to influx, efflux, and growth rate. In a steady state of

011927-3

YOHEI KONDO AND KUNIHIKO KANEKOPHYSICAL REVIEW E 84, 011927 (2011)
3

0

10

-1

10

5

-2

10

X0 = 10
X0 = 0.5-310

density

growth rate

10

1

kg

10

0

2

10

1

10

0

1

(a)

kg

10

10

1

(b)

kg 10
(c)

FIG. 2. (Color online) Dependencies of the steady-state properties on growth reaction constant kg for X0 = 0.5 (blue) and X0 =10 (red/purple) using a catalytic network with N = 24,K = 10. The parameters were set as XN−1 = 0.01, D = 1.0, γ = 1.0, β =; (b) influx: Jin = D0 (X0 − x0 ); and (c) density:
1.0,  = 0.0001, μ0A = μ0B = 1.0, and yA + yB = 1.0. (a) Growth rate: λ = V1 dV
dt
m = i xi .

λm = Jin − Jout ,(7)which gives the relationship among the growth rate, density, influx, and efflux. Then the efficiency η is represented by η = (1 − Jout /Jin )/m. The three order parameters(λ,Jin ,m) are sufficient to specify a cellular state macroscopically.
Figure 2 shows the dependencies of the growth rate λ, influx
Jin , and density m on the growth reaction rate constant kg
for a certain catalytic reaction network with N = 24,K = 10.
(Indeed, this behavior is commonly observed independent of
the networks.) When X0 is large (X0 = 10), it is found that
there are two branches of fixed points with different growth
rates, as shown in Fig. 2(a). In this case, these two fixed points
are bistable when kg takes an intermediate value. The transition
between the two states with the increase or decrease in kg
occurs with hysteresis.
These two states with higher/lower growth rates correspondto lower/higher influx and density, respectively, as shown in
Figs. 2(b) and 2(c). Here, however, the dependence of the influx
on kg is much weaker and there is not a large gap between theinflux values of the two states. Thus the difference between
the growth rates of the two states directly represents the difference in metabolic efficiency η. This significant difference inthe efficiency reflects the internal state of the cell, as discussed
later.
Figure 2 also shows that only one state with the intermediate
growth rate exists when the external nutrient concentration issmall (X0 = 0.5). This state has a very low influx value, andaccordingly, a high efficiency and is therefore distinguishable
from the above two states. When X0 decreases, we observedthat the efficiency increases, which is a common behavior,independent of the values of kg .
Besides the states shown in Fig. 2, another state with anegative growth rate appears when kg is large enough. In thiscase cells shrink, and hence that the state is regarded as a“death” state. We note that the death state coexists with anothergrowing state if X0 is not small, but, for small X0 values (saysmaller than 0.1), it is the only attractor.
Summarizing these distinct states, we have plotted the phase
diagram against kg and X0 as shown in Fig. 3. Each “phase”

here is defined and characterized as in Table I. In particular,
as expected from the plot in Fig. 2, the range of kg values
in which phases I and II coexist narrows continuously as X0
decreases and finally disappears.
B. Mean field analysis

The existence of each phase and the structures of thetransitions between phases are explained by carrying out asimple mean field approximation for our model. To construct
the mean field equations of our model, we consider the
average concentration x of chemicals over all chemicalcomponents (i = 0,1, . . . ,N − 1), i.e., x = m/N, and replace
the rate equation by this concentration. To get this mean
field approximation, the catalytic network is approximatedby a fully connected network with weak paths so that it has
the the same average path density ρ in the original model(=K/N). Also, the standard chemical potentials μ0i of all
the chemicals are set to be identical (to their mean value),and the standard chemical potentials of the energy currencycomponents μ0A ,μ0B are the same. In addition, by assumingthat diffusion through the membrane is fast enough, theinfluent component concentration x0 and effluent component
concentration xN−1 are regarded to be equilibrated with theirexternal values, and they are replaced by X0 and XN−1 ,respectively.
Combining these approximations, we obtain the following
temporal evolution equations, which consist of the concentration variables of the catalyst x and the energy currency
2

10
kg

the cell, there has to be a balance between material influx anddilution by cell volume growth so that the density is constant.
Hence we get

0 IV

10

I or III& IV

Phase I
I & II

III
II

10 -1

X0

101

FIG. 3. Phase diagram of the model on the X0 − kg plane. Each
phase (I, II, III, IV) is defined as in Table I. Solid lines represent
discontinuous transitions between phases; broken lines represent
continuous transitions.

011927-4

GROWTH STATES OF CATALYTIC REACTION NETWORKS . . .

PHYSICAL REVIEW E 84, 011927 (2011)

TABLE I. Definition of phases and characterization by orderparameters, and the statistical distribution of abundances of chemicals
to be discussed in Sec. III C.

diagram agrees rather well with that in the original model(Fig. 3).
From the mean field equations one can see that themean catalyst concentration x (or the density m = N x)
has a negative effect on the growth rate λ; this introduces
bistability in phases I and II. The third and fourth terms
in the right-hand side of Eq. (9) indicate the consumptionof A by the growth reaction and the degradation of A by
catalytic reactions, respectively. The dependence of each termon catalyst concentration x shows that the degradation of
A increases faster than consumption, so that the allocationof energy currency to the growth processes decreases withincreasing x. Hence, for large x the abundance of A and thegrowth rate are suppressed. Consequently, the dilution effect,
represented by the third term in the right-hand side of Eq. (8),does not work effectively and x is sustained to take a large
value. On the other hand, when x is initially small it continues
to take a small value, as it is suppressed by dilution. As aresult, two states with different growth rates coexist. Notingthat m = xN in the mean field approximation and replacing
x by the density m, the above argument is also applied to theoriginal model. To sum up, there is mutual inhibition between
the growth rate and catalyst density. This mutual inhibition
leads to bistability. Note, however, that for this mechanism
to work, the density m should be sufficiently high so thatbistability is not present for small X0 . In this case, a single
phase, III, exists against the change in kg .
Since the mean field equation involves just two variables,
the transitions among each phase are analyzed as bifurcationsof fixed points. From the standard analysis of the flow of a
two-dimensional state space and the linear stability of the fixed
points, it is straightforward to show that there are saddle-nodebifurcations at the transitions from phase I to II and II toI, with kg as a bifurcation parameter. These two bifurcation
points correspond to the boundaries between the phase I andthe bistable region, and the phase II and the bistable region.
As X0 decreases, these bifurcation points approach and finallycollide and disappear. This disappearance corresponds to thetransition to phase III.
Transition to the death phase, IV, has the same bifurcation
structure. There are saddle-node bifurcations at the transition
between the death phase IV and the growing phases with X0 as
a bifurcation parameter. These bifurcation points approach andfinally collide with each other and disappear as kg decreases.

Phase

λ

Distributionof abundancesηJin

I: Active phase
High
High Medium Power law
II: Inefficient phase Low
High Low
Close to equilibrium
III: Quasistatic phase Low
Low High
Mixed
IV: Death phase
Negativemolecule yA , asdx
1
= ρx(X0 yB − xyA ) − ρx(yA + yB )dt
2× (x − XN−1 ) − xkg γ (yA − yB eβ ),dyA= ρ(N − 2)x(X0 yB − xyA )dt
1− ρ(N − 2)x(yA − yB )(x − XN−1 )2− kg (yA − yB eβ ) − ρ(N − 2)2 x 2 (yA − yB ).

(8)

(9)

As with the direct numerical results, the existence of two
branches of fixed points is confirmed for large X0 , withhysteresis as a change of the growth reaction rate constant kg .
In Fig. 4(a), the growth rates are plotted vs kg for X0 = 0.5 and10, which are obtained from the mean field equations, Eqs. (8)and (9). Comparison with the results in Fig. 2(a) indicates thatthe obtained dependence of the two fixed point values on theparameter kg and the existence of hysteresis for large X0 agree
well qualitatively with the results from the original model.
In addition, the death phase, IV, with a negative growthrate, is also obtained from the mean field equations. Again,
for a medium value of X0 , this fixed point coexists with thatfor phase I or III. In contrast, with decreasing kg the region ofcoexistence narrows and finally disappears.
By summing all the solutions, a phase diagram is obtained
from the mean field equations, as shown in Fig. 4(b). The
0

2

10

I or III
& IV

10 -1
0

10 -2
10

10
X0 = 10
X0 = 0.5

-3

1

kg
(a)

10

Phase I
I & II

IV

kg

growth rate

10

III
II

10-1

X0

C. Statistical natures of internal states

101

(b)

FIG. 4. (Color online) Behaviors of the mean field Eqs. (8) and(9)]. The other parameters are set to be identical with those adopted inFig. 2. (a) Dependence of growth rate on kg , obtained from the steadystates of the mean field equations, for X0 = 0.5,10. Colors represent
different X0 values. (b) Phase diagram from mean field equations
on the X0 − kg plane. Each phase is defined in the text as given in
the original model. Solid lines represent discontinuous transitions
between the phases with saddle-node bifurcations, and broken lines
represent continuous transitions distinguished by changes in thegrowth rate and influx without bifurcations.

So far, we have shown the distinct phases characterized
by “macroscopic” variables, such as the growth rate, influx,
and total density of chemicals, which are well represented
by the mean field equations. In a direct simulation involving
many chemicals, however, there are microscopic variables
such as the catalyst concentrations (x0 ,x1 , . . . xN−1 ). In ourmodel, there is no specific meaning to each chemical; hence
the distribution of abundances is a relevant characteristic of acellular state. For this statistical property of the concentrations,
we take a rank-ordered concentration distribution, i.e., weplot the abundances in order of magnitude following earlier
studies [12]. As a result, it is revealed that each phase has aquite different statistical property, as discussed below.

011927-5

10

X0 = 0.6
X0 = 0.3
X0 = 0.1-3∝ (rank)-1

10 -5

Phase I101

0.1
Phase II
0

103

rank

X0 = 0.6
X0 = 0.3
X0 = 0.1

1000
rank

(a)

10

2000

(b)

-2
0.005

-3

10 -4
0

0.001

1

growth rate

concentration

10

PHYSICAL REVIEW E 84, 011927 (2011)

1

10 -1

concentration

concentration

YOHEI KONDO AND KUNIHIKO KANEKO

rank 30

0.002
0.001

Phase II
0

Phase III
1000
rank

Phase I

0

2000

2

(d)

KL divergence

(c)

1
X0

10

1

100

10 -1
10

-2

10

-4

-3

10 -5

10

0.1

0.1 k g 10

0.3

0.5

0.7

X0
(e)

FIG. 5. (Color online) Statistical nature of the model. A networkwith N = 2000 and K = 20 was used. (a) Rank-ordered concentration distributions in phase I for X0 = 0.1,0.3,0.6. Colors represent
different values of X0 . Here kg = 10.0, and the other parameters aresame as those in Fig. 2. (b) Distributions in phase II. The networkand parameters are the same as those in (a), except kg = 0.1. (c) A
distribution in phase III. The network and parameters are the same as
those in (a), except for X0 = 0.02,kg = 5.0. The inset shows a log-log
plot of the distribution up to rank 30. (d) Dependencies of the growthrate on X0 in phases I and II. Colors and shapes of points represent
data corresponding to the plots in (a) and (b) sharing the same colorand shape. (e) Dependencies of Kullback-Leibler divergence from the
distributions to the equilibrium compositions on X0 . The inset shows
the dependencies on kg .

Figures 5(a)–5(c) show rank-ordered concentration distributions in phases I, II, and III, respectively. The ordinate
indicates the concentration of catalysts, and the abscissa shows
the rank determined by the concentration. For phase I weobserved that the slope of the rank-ordered concentration
distribution increases with increments of X0 and reaches a
power-law distribution with exponent –1. However, we notethat this power-law distribution with an exponent –1, in otherwords, Zipf’s law, can be achieved only when the catalytic
reaction network is sufficiently sparse (approximately K < 40
for N = 2000). For a dense network, although the slope of the
distribution increases with X0 , the distribution remains curved
in the log-log plot and thus does not fall on the power law.
The change in the slope is accompanied by a change in growthrate, which, for reference, is plotted in Fig. 5(d) with increasing
X0 for phase I. Each point in the figure corresponds to each
distribution in Figs. 5(a) and 5(b) of the same color and shape.

With increasing the growth rate, the slope of the distribution
also increases and approaches the power law with exponent–1, near the optimal growth rate. With further increases in
X0 , however, the growth rate decreases until it reaches zero.
In this region before the arrest of growth, the concentrations
approximately obey a power-law distribution with exponents–1 or slightly higher (data not shown). In this region, however,the steady state is often no longer stable and the growth rate
shows oscillations with time. The plot with the downward
arrow in Fig. 5(d) denotes a regime with such oscillations.
Figure 5(b) shows rank-ordered distributions in phase II for
several values of X0 by a semilog plot. In contrast to phase I,
the plot is linear on a semilog scale and the slope is independentof the values of X0 . In other words, the distribution is
exponential, which, as shown below, comes from the factthat the concentrations have near equilibrium compositionsin these states, defined by the standard chemical potentials
eqof the catalysts as xi ∝ exp(−βμ0i ). Because of the uniform
distribution of the potentials μ0i , the rank-ordered equilibrium
distributions take the form of exponential distributions.
To confirm that the distributions in phase II are actuallyclose to the equilibrium distribution, we measure the distancefrom equilibrium by calculating the Kullback-Leibler divergence from the equilibrium distributions of the concentrations
eq
(eq)
(eq)
(eq)
(= i pi log(pi /pi ), where pi ≡ xi / i xi ) [23], as
plotted in Fig. 5(e), vs the change in X0 ; the change of thedivergence vs kg is shown in the inset of Fig. 5(e), ranging
between phases I and II. It is shown that the divergence inphase II is quite low compared to that in phase I and decreases
with increments of X0 , to approach equilibrium. Although
the emergence of exponential distributions itself depends on
the assumption of uniform distribution of μ0i s, the decreases
in the divergence are independent of assumptions on the
distribution of μ0i s. Thus, phase II can be characterized by
the distribution close to the thermal equilibrium. As shownin the inset of Fig. 5(e), this branch connected continuously
from the equilibrium shows a discontinuous jump to phase I,
where the composition of the concentration distribution is far
from that at thermal equilibrium.
Thus phases I and II are characterized by distributions withpower-law distributions and close to equilibrium, respectively.
On the other hand, the states in phase III are mixed, as
shown in Fig. 5(c). For chemicals with smaller abundances therank-abundance plot follows exponential decrease and is nearequilibrium, as in phase II. The concentrations of abundantchemicals, however, deviate distinctly from the equilibrium
distribution. The rank-abundance plot for abundant chemicals
follows a power law, as in phase I, while the fraction of suchchemicals that deviate from equilibrium distribution increases
with kg . This partial deviation from equilibrium distribution
characterizes phase III.
IV. SUMMARY AND DISCUSSION

We have studied a simple “protocell” model in which
chemicals are encapsulated within a membrane. The protocell
synthesizes these chemicals by using the nutrients it derives
from the environment; the protocell also utilizes the energy
obtained from intracellular reactions to increase its volume.
These reactions are catalyzed by other chemicals that are also

011927-6

GROWTH STATES OF CATALYTIC REACTION NETWORKS . . .

PHYSICAL REVIEW E 84, 011927 (2011)

synthesized within the cell, and thus the total set of reactions
constitutes a catalytic network, which was randomly chosen
in this study. From the results of simulations of a variety of
networks, we have obtained the steady state in which a cellgrows with a constant rate.
We have carried out simulations by varying the values of theconcentration of the nutrient outside the cell, i.e., X0 , and therate constant of the phenomenological growth reaction, i.e.,
kg . Results of these simulations show that three distinct phaseswith qualitatively different cellular states exhibiting steady
growth exist, depending on these two parameters. When theexternal nutrient concentration X0 is high, phase I appears fora large value of kg ; this phase is characterized by high cell
growth and influx rates. On the other hand, when the value of
kg is small, phase II appears, characterized by a low growthrate and a high influx rate. At an intermediate value of kg , thesetwo phases coexist. When X0 is low, the influx rate decreases,
and this regime where phases I and II coexist narrows and thencompletely disappears. After the disappearance of phases I andII, phase III comes into existence; this phase is characterized
by low growth and influx rates, regardless of the values ofkg . Further, this phase is characterized by a high efficiency of
conversion of nutrients, which enables an increase in the cellvolume. In addition to these states characterized by a positive
growth rate, we also found phase IV, termed the “death” phase;
this phase is characterized by negative growth rates for low X0 .
When kg is high, it is found that both the death phase and thegrowth phase I or III coexist.
These transitions between the phases are explained using
the mean field equations of the model, which are obtainedby assuming that the path density is sufficiently high. The
mean field equations account well for the transition between
phases I and II, as well as those between phase I or III and
IV as saddle-node bifurcations. The success of the mean field
equations that are independent of each network also explains
why the behaviors observed here are independent of the choiceof networks. The analysis of the mean field equations enable us
to understand the bistability of two phases as mutual inhibition
between the cell volume growth and the concentration of
intracellular chemicals. Because this mutual inhibition does
not depend on the details of the model, we expect that thisbistability will be a general feature of any cells in which
the rate of conversion of nutrient for increasing cell volume
can vary depending on the composition of chemicals in thecell.
In the original simulations of many chemical components,
we have also examined the statistical distribution of theabundance of chemicals. To characterize each cellular phasestatistically, we adopt the distribution of abundances. In phase
I, the abundance of chemicals differs from species to species,
spreading over several orders, and this distribution is far
from that expected from the equilibrium distribution definedby chemical potentials. As the growth rate increases, theabundance distribution of chemicals approaches a certain
power law, known as Zipf’s law. On the other hand, theabundance distribution of chemicals in phase II is close to theequilibrium distribution. Hence the transition from phase II to I
is a bifurcation from a branch near equilibrium to that far fromequilibrium. It should be noted that the power law with respect
to the abundance distribution was also observed in the previous

model, without consideration of energy conversion [12], whilephases II, III, and IV were first observed in this study.
It is interesting to note that the power-law behavior
exhibited in the protocell model in this study is also commonly
observed in the existing cells with regard to gene expression
levels [12]. In particular, it is confirmed in the log phasein bacteria where cells grow efficiently as in the phase
I in our model. Note that the power-law abundances inprotocell models is observed for a state in which a cell grows
efficiently, keeping its chemical composition. The distributionof abundances over chemicals may provide a proper statisticalmeasure for a cellular state.
Because of the simplicity of our model, it is not so easy to
directly compare our results with the behaviors of the existing
cells. However, it is still interesting to note the similarity ofthe phases observed in this study with the log, stationary, anddormant phases. The log phase appears under a nutrient-richcondition, where the cell has a high growth rate with a highinflux rate. The stationary phase appears when a nutrient
becomes exhausted or toxic waste products become excessive
and the cell growth slows. In the dormant phase the cell growthrate is suppressed even under a nutrient-rich condition [5].
These behaviors are reminiscent of phases I, III, and II,
respectively. In this respect, it will be interesting to check
whether there is hysteresis between the transitions from the log
phase to the other phases at the single-cell level and to measure
the efficiency of cell growth. Mutual inhibition between the cell
growth and catalyst abundances to support the bistability of thetwo phases can be examined directly. As mentioned above, statistical distribution of chemical abundances is distinct by eachphase, which may be examined from gene expression analysis.
So far, no reliable microarray data are available for stationary
and dormant phases, which should be important in the future.
The present model does not assume a specific reaction
network and intracellular spatial organizations. Currently,
several attempts are being made to construct an artificial cellthat grows and replicates [24], and the basic structures of thepresent model are often satisfied with such an artificial cell.
Hence it will be interesting to search for phases by measuring
the growth rate and the influx rate for cell growth.
In general, there is no established theory for thermodynamics far from equilibrium. The metabolic efficiency η discussed
above is also not, in a strict sense, the thermodynamic efficiency. However, by restricting the theory to a steadily growing
state, we may hope to develop a thermodynamic theory on thebasis of which the efficiency of cell growth is analyzed [25].
The present model exhibits three distinct phases with regard to
the influx and growth rates, with unique metabolic efficiencies.
Such a thermodynamic theory, if established, will account forthese phases and efficiency of cell growth.
ACKNOWLEDGMENTS

We would like to thank Tetsuya Yomo, Shuji Ishihara,
Chikara Furusawa, and Stuart Kauffman for illuminating
comments and stimulating discussions. The present work
is partially supported by the ERATO Dynamical MicroScale Reaction Environment Project by Japan Science andTechnology Agency (JST). We are grateful to the Japanesepeople for generous support.

011927-7

YOHEI KONDO AND KUNIHIKO KANEKOPHYSICAL REVIEW E 84, 011927 (2011)[1] E. Schrödinger, What is Life? (Cambridge University Press,
Cambridge, England, 1944).
[2] G. Nicolis and I. Prigogine, Self-Organization in Nonequilibrium
Systems (Wiley & Sons, New York, 1977).
[3] J. B. Russell and G. M. Cook, Microbiol. Rev. 59, 48 (1995).
[4] J. Monod, Annu. Rev. Microbiol. 3, 371 (1949).
[5] N. Q. Balaban, J. Merrin, R. Chait, L. Kowalik, and S. Leibler,
Science 305, 1622 (2004).
[6] I. Keren, D. Shah, A. Spoering, N. Kaldalu, and K. Lewis, J.
Bacteriol. 186, 8172 (2004).
[7] M. Eigen and P. Schuster, The Hypercycle (Springer, 1979).
[8] T. Gánti, BioSystems 7, 15 (1975).
[9] R. J. Bagley and J. D. Farmer, in Artif. Life II, edited by C.
G. Langton et al. (Addison Wesley, Redwood City, CA, 1991),
p. 93.
[10] D. Segrè, D. Lancet, O. Kedem, and Y. Pilpel, Origins Life Evol.
Biosphere 28, 501 (1998).
[11] S. Jain and S. Krishna, Proc. Natl. Acad. Sci. USA 98, 543(2001).
[12] C. Furusawa and K. Kaneko, Phys. Rev. Lett. 90 088102 (2003).
[13] C. Furusawa, T. Suzuki, A. Kashiwagi, T. Yomo, and K. Kaneko,
Biophysics 1, 25 (2005).

[14] K. Kaneko, Adv. Chem. Phys. 130, 543 (2005).
[15] C. Furusawa and K. Kaneko, Phys. Rev. E 73, 011912(2006).
[16] T. Rocheleau, S. Rasmussen, P. E. Nielsen, M. N. Jacobi, and H.
Ziock, Philos. Trans. R. Soc. B 362, 1841 (2007).
[17] T. Carletti, R. Serra, I. Poli, M. Villani, and A. Filisetti, J. Theor.
Biol. 254, 741 (2008).
[18] M. Ederer and E. Gilles, Biophys. J. 92, 1846 (2007).
[19] A. Awazu and K. Kaneko, Phys. Rev. E 80, 041931(2009).
[20] S. A. Kauffman, J. Theor. Biol. 119, 1 (1986).
[21] P. F. Stadler, W. Fontana, and J. H. Miller, Physica D 63, 378(1993).
[22] R. Hanel, S. A. Kauffman, and S. Thurner, Phys. Rev. E 72,
036117 (2005).
[23] E. Smith, J. Theor. Biol. 252, 198 (2008).
[24] Protocells: Bridging Nonliving and Living Matter, edited by S.
Rasmussen, M. A. Bedau, L. Chen, D. Deamer, D. C. Krakauer,
N. H. Packard, and P. F. Stadler (MIT Press, Cambridge, MA,2009).
[25] H. V. Westerhoff, K. J. Hellingwerf, and K. V. Dam, Proc. Natl.
Acad. Sci. USA 80, 305 (1983).

011927-8