PRL 108, 208103 (2012)

week ending
18 MAY 2012PHYSICAL REVIEW LETTERS

Adaptation to Optimal Cell Growth through Self-Organized Criticality
Chikara Furusawa
Quantitative Biology Center, RIKEN, 6-2-3 Furuedai, Suita, Osaka 565-0874, Japan

Kunihiko Kaneko
Research Center for Complex Systems Biology, University of Tokyo, Komaba, Meguro-ku, Tokyo 153-8902, Japan(Received 22 October 2011; published 15 May 2012)
A simple cell model consisting of a catalytic reaction network is studied to show that cellular states are
self-organized in a critical state for achieving optimal growth; we consider the catalytic network dynamics
over a wide range of environmental conditions, through the spontaneous regulation of nutrient transport
into the cell. Furthermore, we find that the adaptability of cellular growth to reach a critical state depends
only on the extent of environmental changes, while all chemical species in the cell exhibit correlatedpartial adaptation. These results are in remarkable agreement with the recent experimental observations of
the present cells.
DOI: 10.1103/PhysRevLett.108.208103

PACS numbers: 87.17.Aa, 82.39.k, 87.10.e

Cells flexibly change their intracellular states withrespect to environmental changes, so that they can sustain
cellular growth under a variety of environmental conditions. Experimental analysis of microorganisms suggested
that metabolic reactions can be organized to achieve optimal or nearly optimal growth under various environmental
conditions [1]. Although environmental adaptation is generally explained by using signal transduction machineries
as lac operon [2], it is difficult to imagine that such
regulatory networks are prepared by evolution against unexpected environmental changes. Hence, it is natural toexpect that cells will exhibit a rather generic, ubiquitous
process of environmental adaptation, even in their primitive form. Given that all cellular processes are a result ofcomplex chemical dynamics of catalytic reactions, to discuss the generic adaptation process of (primitive) cells, it isimportant to study how such a system with complex catalytic reaction networks can adapt to environmental changes
and maintain optimal or nearly optimal growth.
In this Letter, we show that a simple active transport ofnutrients enables generic adaptation toward a state with
optimal growth, by simulating a simple cell model with acatalytic reaction network [3]. In our previous study [4],we showed that recursive cellular reproduction with optimal growth is achieved by controlling the flow rate of thenutrients, and we found that such optimal growth isachieved at a ‘‘critical’’ cellular state, in which the rank
distribution of chemical abundance obeys a power law withexponent 1 (i.e., Zipf’s law). Here, by introducing asimple active transport of nutrients, we show that thecellular dynamics are able to sustain optimal cellular
growth by self-organization into a critical state [5].
In this study, we adopt a simple model of a reproducing
cell with intracellular catalytic reaction networks thattransform nutrient chemicals into other chemical species
that exhibit catalytic activity [4,6]. Let us consider a cell
0031-9007=12=108(20)=208103(5)consisting of N chemical species. The state of this cell canbe represented by a set of chemical concentrations(x0 ; x1 ; . . . ; xN1 ), where xi is the intracellular concentration of chemical species i. Depending on whether an
enzymatic reaction occurs from i to j, catalyzed by someother chemical ‘, the reaction path is connected as ði þ
‘ ! j þ ‘Þ. The rate of increase in xj and decrease in xi
through this reaction is given by xi x‘ . For the sake ofsimplicity, all the reaction coefficients have been set to 1.
The connection paths of this catalytic network have beenchosen randomly, where for all combinations of substrate iand product j the reaction path is randomly connected witha fixed probability .
Next, some chemicals are supplied from the environment as nutrients. The transport of the chemical nutrients
into the cell is activated with the aid of some otherchemicals, i.e., ‘‘transporters.’’ Here, the uptake flux ofthe ith chemical from the environment is assumed to be
represented by Fi ¼ Sxwi  , where the wi th chemical acts
as the transporter for chemical i and xwi is the concentration of the wi th chemical. Parameter S denotes theconcentration of the nutrient chemicals outside the celland is a transport constant. Parameter  represents thehigher-order effect in the nutrient transportation. If thisparameter is greater than 1, the transporter catalysis is ofa higher order, which corresponds to the formation ofmultimeric complex of the transporter. It should be notedthat the cell volume increases with the abundance ofchemicals within the cell. For the sake of simplicity, we
assume that the cell volume is directly proportional to thetotal abundance of chemicals in the cell. The concentration of each chemical species is diluted with an increasein the cell volume by the nutrient transport.
P This assumption corresponds to imposing restriction N1i¼0 xi ¼ const,
which can be set to 1 by appropriately normalizing theconcentration variables.

208103-1

Ó 2012 American Physical Society

N1
N1
X
X
dxi
¼
Wji‘ xj x‘ 
Wij0 ‘0 xi x‘0dtj;‘¼0j0 ;‘0 ¼0N1X

Fj ðþFi Þ;

(1)

j¼0

where Wij‘ is 1 if reaction i þ ‘ ! j þ ‘ occurs; otherwise, it is 0. The third term with
P the sum of nutrient uptake
flux gives the constraint of N1i¼0 xi ¼ 1 because of thevolume growth. The last term is added only for the nutrients, and it represents their transportation into a cell fromthe environment. For simplification purposes, we consider
a cell with only one nutrient chemical and a correspondingtransporter, whose concentrations are x0 and x1 , respectively. Here, the cell growth rate is equal to the (total)
nutrient uptake F0 ¼ Sx1 . We also assume that the nutrientchemical does not exhibit catalytic activity and cannot be
the transporter.
Before describing the results of the present model, we
note that, in the earlier model that did not include a transporter as in the case when  ¼ 0, we found that thetransition occurred at a certain critical value of S, beyond
which the cell could no longer grow continuously [4].
When the nutrient uptake rate exceeds the critical value,chemical reactions that transform the nutrient into otherchemical substances, including catalysts, cannot follow thenutrient uptake. In such a case, the reaction network sustaining ‘‘metabolism’’ becomes unstable, which in turn
eventually ceases the cell growth. This finding indicatesthat the critical nutrient uptake rate corresponds to the
maximal capacity of the reaction network to assimilate
the nutrient for recursive reproduction. Results of our study
indicate that, at the critical nutrient uptake rate, the chemical abundance distribution obeys a power law, which ismaintained by a hierarchical organization of enzymatic
reactions. Furthermore, it should be noted that a maximal
growth rate is achieved at this critical nutrient uptake rate.
Here, we numerically study the present model with  1 by considering a variety of randomly chosen networks. Figure 1 shows a plot of the nutrient uptake F0 (thatcorresponds to the growth rate) at the steady state against
parameter S for the cases with  ¼ 2. For most cases, thesystem reaches a unique fixed-point attractor (see Fig. S1[7] for exceptional cases). As shown in this figure, thenutrient uptake F0 with a low concentration (S < 103 )increases proportionally with S, indicating that the concentration of transporter x1 is almost unchanged. In contrast,
when the concentration of the nutrients in the environment
is high (S > 1), the nutrient uptake is sustained at a certain
level, irrespective of S. This result implies that theconcentration of the transporter decreases with an increase
in S, in proportion to S1= . On the other hand, x0 increases

with S, and it was numerically found that x1 / ð1  x0 Þ ,where  depends on each network wherein  ranges from
2 to 5 (see the inset in Fig. 1). It is important to note that, in
this region, the nutrient uptake is almost identical to the
optimal nutrient uptake achieved only at the criticalnutrient uptake for  ¼ 0 (denoted by the dotted line).
In other words, over a wide range of parameter S, thecellular dynamics including the transporter concentration
are self-regulated to achieve the critical state with an
optimal-growth rate. This self-regulation of the transporter
concentration is always achieved for any network and forany value of   1, as long as the environmental nutrient
concentration is higher than some value (Fig. S1 [7]).
The phenomenon of self-regulation to achieve criticality
is explained as follows. When S is increased and when the
nutrient uptake exceeds the critical level for reproduction,
the nutrient concentration in a cell is too large; therefore,
there is no room for producing all other non-nutrientchemicals [4]. Hence, the concentration of the transporter
also decreases, which in turn causes the nutrient uptake
rate to decrease down to the critical value. Here, thenutrient uptake regulates both the flow of chemical reaction
from the nutrient and the dilution of chemical concentration due to the cellular growth. By achieving a balance
between the production and dilution of the chemicals by
negative feedback through the transporter concentration,
the adaptive dynamics toward a state with a nearly optimal
growth are self-organized.

(a)
10-2

10-3

(b)
-4

10

10-4
x1

To summarize these processes, the dynamics of theintracellular concentration of the ith chemical, where i isin the range i ¼ 0; 1; . . . ; N  1, are represented as xi

week ending
18 MAY 2012PHYSICAL REVIEW LETTERS

nutrient uptake F0

PRL 108, 208103 (2012)

-5

10

-5

10

10-6

0.2

-6

10

-4

10

-2

0.5

1

2

10

1-x0

10-6
0

10
10
10
nutrient concentration S

4

FIG. 1 (color online). Adaptive dynamics of nutrient uptake.
(a) Nutrient uptake F0 plotted as a function of the environmental
nutrient concentration S. The results of ten randomly generated
reaction networks are overlaid. The dotted line denotes the
maximal nutrient uptake for the case  ¼ 0, beyond which thecell can no longer grow continuously. The parameters were set as
N ¼ 10 000,  ¼ 2:5  102 , and  ¼ 2. (b) The concentrationof transporter x1 is plotted as a function of the total concentrationof non-nutrient chemicals 1  x0 . The relationships obtained byten random networks shown in (a) are presented. The solid line
represents a slope of 2 for reference.

208103-2

PHYSICAL REVIEW LETTERS

PRL 108, 208103 (2012)

Figure 2 shows a plot of the rank-ordered distributions of
the chemical concentrations obtained by a specific networkfor  ¼ 2, where the ordinate shows the chemical concentration xi and the abscissa shows the rank determined by
the magnitude of xi . When the nutrient concentration S islow, the internal reaction dynamics progress faster than theflow of the nutrient from the environment, all the chemicals
have nearly identical concentrations, and the influence of S
is limited to a finite range of chemical species. For S > 10,
where the nutrient uptake is close to the critical level, thedistributions start to approximately follow a power law.
Although the less abundant chemicals (e.g., rank >5000)
exhibit some deviation, the exponent of the power law isclose to 1 [8]. Indeed, such a power law was observed at
the critical point with the maximal nutrient uptake in ourprevious study for  ¼ 0. This power law distribution withexponent 1 is attributed to the hierarchical organizationof catalytic reactions, in which chemicals with low concentrations are produced from those with high concentrations [4]. In the present model, this power law holds aslong as S is sufficiently high so that the nutrient uptake
is self-regulated to a critical level. This confirms the selforganization of a cellular state to a critical state with
optimal nutrient uptake.
This self-regulation to the critical state is also estimatedby a simple ‘‘mean-field’’ calculation, by assuming thehierarchical flow of chemical reactions. Here, we consider
a hierarchical layer of chemical species from the nutrient,as L0 ! L1 ! L2 ! . . . Lk , where L0 denotes the nutrientchemical within the cell and Li denotes a set of chemical
species synthesized from Li1 chemicals. We introduce mithat denotes the mean concentration of Li chemicals andassume that the concentrations of all enzymatic chemicals
can be approximated by an identical mean-field value010

S=1x10 -3
S=1x10 -2
S=1x101-1
S=1x10 2
S=1x10 3
S=1x10-1chemical concentration

10

-2

10

-3

10

10-4
10-5
10-6
0

10

1

10

2

10
rank

3

10

4

10

FIG. 2 (color online). Rank-ordered concentration distributionof chemical species. Distributions for different values of S aresuperimposed by using different colors. The solid line has aslope of 1 for reference. The other parameter values are
identical to those shown in Fig. 1.

week ending
18 MAY 2012

ð1  m0 Þ=N, the mean concentration of all non-nutrientchemicals [9]. Based on this mean-field approximation andby assuming that the transporter chemicals belong to thekth layer, we obtain
dm0 =dt ¼ Smk  ðM=NÞð1  m0 Þm0  Smk m0 ;(2)

dmj =dt¼ðM=NÞ½ð1m0 Þmj1 ð1m0 Þmj Smk mj ;(3)where j ¼ 1; . . . ; k and M ¼ N denotes the mean numberof reaction paths from each chemical. The last term in thisequation represents the dilution of each chemical, attributed
to the flow of nutrients. Then, we assume the steady stateof chemical concentrations, i.e., dmj =dt ¼ 0 for j ¼
0; 1; . . . ; k. From dm0 =dt ¼ 0, we obtain F0 ¼ Smk ¼
ðM=NÞM0 ¼ m0 . By substituting it into Eq. (3) andassuming dmj =dt ¼ 0, we get mj ¼mj1 ð1m0 Þ. By iteratively applying this relationship, we get mk ¼m0 ð1m0 Þk .
From these equations, the steady values of F0 , mk , and m0
are obtained as functions of S. In particular, when S is large,
mk / ð1  m0 Þk follows; this result is consistent with thenumerical observation shown in Fig. 1(b) [10]. Furthermore,
we obtain F0  , indicating that the nutrient uptake F0
is independent of the nutrient concentration when S islarge. Indeed, this relationship F0   does not differ
much from the numerical observation, as shown in Fig. 1.
In the Supplemental Material [7], we provide an explanation
for the power law distribution of chemical abundance withexponent 1, as shown in Fig. 2, by considering theconservation of chemicals and by using the above meanfield approximation.
The self-regulation of the nutrient uptake rate to thecritical value indicates that adaptive dynamics restore thenutrient uptake to the original level after environmental
changes [11]. Figure 3 shows the adaptive dynamics of thenutrient uptake F0 against the change in S from S0 to S00 .
The response profile of F0 to the environmental changes
depends only on ratio S00 =S0 . This type of adaptive response, termed fold-change detection, was recently confirmed in cellular responses [12]. Simple models of gene
regulation or catalytic reaction involving few chemicalcomponents have been proposed to show such a type of
adaptive dynamics [13,14]. In the present case, this foldchange detection is also shown from the mean-field equation in the large S limit. For example, let us consider therelaxation dynamics of Eqs. (2) and (3) with k ¼ 2 toward
the fixed point [m0  1  ð=SÞ1=2 , m2  ð=SÞ1=
ð1  =ð=SÞ1=2 ]. (Because m1 ¼ 1  m0  m2 , there
are only two independent variables.) For large S, theeigenvalues of the Jacobi matrix around the fixed point
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
are given by ð1  1  8=Þ=2 and are independentof S. On the other hand, when S is changed from S0 to S00 ,the initial change in the nutrient uptake F0 ¼ Smk issolved from Eq. (1) to give 1=F0 ðtÞ  ðS0 =S00 Þ=F0 ð0Þ ¼

208103-3

PHYSICAL REVIEW LETTERS

PRL 108, 208103 (2012)
0.05

100
200
400

nutrient uptake F0

0.04

200
400
800

0.03

0.02

0.01
-2000

0

2000

4000
time

6000

8000

10000

FIG. 3 (color online). Adaptation dynamics of the nutrient
uptake F0 after the environmental nutrient concentration S ischanged. At t ¼ 0, S is increased twice, and the responses
corresponding to S ¼ 100 ! 200, S ¼ 200 ! 400, and S ¼
400 ! 800 are overlaid. The other parameter values are identicalto those shown in Fig. 1.

ð1  S0 =S00 Þt (for small t and large S0 ). Thus, the timecourse F0 ðtÞ depends only on the fold change S00 =S0 ; therefore, fold-change detection is demonstrated by using themean-field approximation.
This perfect adaptation of the nutrient uptake (i.e., thecellular growth rate), as shown in Fig. 3, is supported by alarge number of catalytic reactions through negative feedback control by a transporter. Besides the growth rate, theconcentrations of all the chemical components (that areneither a nutrient nor a transporter) are changed to support
the collective adaptive dynamics of the growth rate, against
the change in S. It should be noted that, even though all thechemical components do not show perfect adaptation to
return to the original value as in the case of the growth rate,
in general, all the chemical components show partial adaptation [15], as shown in Fig. S2 [7]. It is interesting to notethat such a partially adaptive response is observed for many
gene expressions against environmental changes [16].
In conclusion, we have demonstrated that a cellular
system consisting of a catalytic network shows optimal
growth against environmental changes when nutrients areactively transported into the system with the help of acatalyst. It should be emphasized that the results reported
herein are independent of the details of the modeling. They
are also true for different kinetic forms in the catalytic
reaction, different topologies of the reaction network suchas the scale-free [17], and distributed reaction coefficients.
For example, we confirmed that the use of MichaelisMenten reaction kinetics with distributed coefficients
does not change the present results, as long as a certain
fraction of the reactions is in the linear regime. In fact,
existing organisms would adopt sophisticated machineries
in metabolic reactions, transcription, and translation.
However, the self-organizing behavior for the adaptive

week ending
18 MAY 2012

growth regulation proposed herein emerges as long as thenutrient uptake dynamics are related to the intracellular
dynamics and thus are expected to generally exist behindsuch sophisticated regulatory reaction processes. In fact, it
is well known that the nutrient uptake in microorganisms is
generally regulated by the intracellular reaction dynamics,
such as a phosphotransferase system [18], which maintains
a balance between the nutrient uptake and the metabolic
activity. Furthermore, the synthetic activity of transporter
proteins is also linked to the metabolic activity, and suchprocesses generally include the higher-order effect as
dimer and tetramer formation, proving that parameter  2 in our model. Furthermore, it is interesting to notethat this adaptation to optimal growth, power law in abundance distribution (Zipf’s law), adaptive response with
fold-change detection, and partial adaptation of many
gene expressions are indeed generally observed in thepresent cells, supporting our hypothesis. We expect thatthe present self-organized criticality to the optimal-growthstate will ubiquitously underlie the flexible and robust
adaptation in cellular systems.
The authors thank T. Yomo, M. Inoue, Y. Kondo, and S.
Ishihara for useful discussions. The present study is partially supported by the ERATO project, JST, and Grants-inAid for Scientific Research (No. 23680030, No. 23650152,
and No. 21120004) from MEXT Japan.

[1] R. U. Ibarra, J. S. Edwards, and B. O. Palsson, Nature(London) 420, 186 (2002); S. S. Fong and B. O. Palsson,Nat. Genet. 36, 1056 (2004); Y. Shinfuku, N. Sorpitiporn,
M. Sono, C. Furusawa, T. Hirasawa, and H. Shimizu,
Microb. Cell Fact. 8, 43 (2009).
[2] F. Jacob and J. Monod, J. Mol. Biol. 3, 318 (1961).
[3] A. Filisetti, R. Serra, T. Carletti, M. Villani, and I. Poli,
Eur. Phys. J. B 77, 249 (2010); S. Jain and S. Krishna,
Proc. Natl. Acad. Sci. U.S.A. 98, 543 (2001); D. Segre, D.
Ben-Eli, and D. Lancet, Proc. Natl. Acad. Sci. U.S.A. 97,
4112 (2000).
[4] C. Furusawa and K. Kaneko, Phys. Rev. Lett. 90, 088102(2003).
[5] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. A 38, 364(1988); D. Dhar and R. Ramaswamy, Phys. Rev. Lett. 63,
1659 (1989); H. J. Jensen, Self-Organized Criticality:
Emergent Complex Behavior in Physical and Biological
Systems (Cambridge University Press, Cambridge,
England, 1998); A. Awazu and K. Kaneko, Phys. Rev. E
80, 010902(R) (2009);[6] C. Furusawa and K. Kaneko, Phys. Rev. E 73, 011912(2006); Y. Kondo and K. Kaneko, Phys. Rev. E 84, 011927(2011).
[7] See Supplemental Material at http://link.aps.org/
supplemental/10.1103/PhysRevLett.108.208103
forFigs. S1 and S2 and supplemental text.
[8] The slopes in Fig. 2 obtained by the least squares fit are1:04, 1:12, and 1:21 for S ¼ 101 , 102 , and 103 ,respectively. For larger S, they might be slightly smaller

208103-4

PRL 108, 208103 (2012)PHYSICAL REVIEW LETTERS

than 1, which is predicted by the mean-field approximation discussed later, presumably due to the finite size
effect of the reaction network.
[9] Hence, the mean-field version differs from a onedimensional chain of catalytic reactions, where adaptation
to the critical state is not observed.
[10] In the simulations shown in Fig. 1, the transporter chemical
is located at L2 or L3 . The reason why the slopes in Fig. 1(b) are sometimes greater than 3 is presumably because of
the hierarchical structures of the catalysts, which are omitted in the mean-field calculation in Eqs. (2) and (3).
[11] D. E. Koshland, Jr., A. Goldbeter, and J. B. Stock, Science
217, 220 (1982); F. Oosawa and Y. Nakaoka, J. Theor.
Biol. 66, 747 (1977); S. Asakura and H. Honda, J. Mol.
Biol. 176, 349 (1984); U. Alon, M. G. Surette, N. Barkai,and S. Leibler, Nature (London) 397, 168 (1999).
[12] L. Goentoro and M. W. Kirschner, Mol. Cell 36, 872 (2009);
C. Cohen-Saidon, A. A. Cohen, A. Sigal, Y. Liron, and U.
Alon, Mol. Cell 36, 885 (2009); M. D. Lazova, T. Ahmed,
D. Bellomo, R. Stocker, and T. S. Shimizu, Proc. Natl.
Acad. Sci. U.S.A. 108, 13 870 (2011).
[13] O. Shoval, L. Goentoro, Y. Hart, A. Mayo, E. Sontag, andU. Alon, Proc. Natl. Acad. Sci. U.S.A. 107, 15995 (2010).

week ending
18 MAY 2012[14] M. Inoue and K. Kaneko, Phys. Rev. Lett. 107, 048301(2011)[15] Adaptive dynamics of the cellular state is generally classified into ‘‘perfect adaptation,’’ in which a variable inconcern returns to the original value after a response, and‘‘partial adaptation,’’ in which the variable shows a change
to decrease the response, but it does not return to theoriginal value completely. See Koshland, Goldbeter, andStock’s paper in Ref. [11] for further details.
[16] A. P. Gasch, P. T. Spellman, C. M. Kao, O. Carmel-Harel,
M. B. Eisen, G. Storz, D. Botstein, and P. O. Brown, Mol.
Biol. Cell 11, 4241 (2000); S. Stern, T. Dror, E. Stolovicki,
N. Brenner, and E. Braun, Mol. Syst. Biol. 3, 106 (2007).
[17] When the transporter is directly produced from thenutrient, feedback to the nutrient uptake often leads tooscillatory behavior, even in the case when   2.
However, as long as the transporter is not directly produced from the nutrient, the adaptive dynamics described
here, as well as the power law in abundances, generally
appears regardless of the position of the transporter in thereaction network.
[18] J. Deutscher, C. Francke, and P. W. Postma, Microbiol.
Mol. Biol. Rev. 70, 939 (2006).

208103-5

Supplementary Text: A sketch of an explanation on Zipf ’s law from the meanfield analysis

The origin of power-law in abundances is due to hierarchical organization of abundancesof chemicals. Although complete explanation on the organization of power-law is out of thescope of the simple “mean-ﬁeld” calculation, here we provide a sketch for the explanation.
From the mean-ﬁeld analysis in the text, we get that the mean concentration of a chemical
at k-th layer obeys mk = m0 (1 − m0 )k . On the other hand, at each k-th layer, there are∼ (ρN )k chemical species. Hence, the ranking of the chemical at k-th layer, denoted by
rk , in the order of abundances increases as rk ∼ (ρN )k when ρN is enough large, and thus
k = log(rk )/ log(ρN ). From these equations, we obtain m(rk ) = m0 (1 − m0 )log(rk )/ log(ρN ) ,where m(rk ) represents the chemical concentration of rk -th ranked chemical. Thus,
log mk = log m0 − α log(rk )(S1)with α = − log(1 − m0 )/ log(ρN ), so that the abundances against the ranking follows apower law with the exponent α.
Note that the total abundances of chemicals at each layer is roughly given by m0 C K with
C = (1 − m0 )(ρN ). The total sum of the abundances up to the maximal layer K should be
1 − m0 , so that
m0 C

(C K − 1)∼ (1 − m0 ).
C −1

(S2)

Now, if C > 1, the left hand side of eq.(S2) diverges as K increases, and eq.(S2) cannot besatisﬁed. On the other hand, if C < 1, eq.(S2) in the large K limit would imply m0 C/(1 −C) = 1 − m0 , which could be satisﬁed only if ρN = 1, i.e., no branching at all between
layers. Thus, for the self-consistency condition to satisfy eq.(S2) for the multiple branchingcase (i.e., ρN > 1), C = (1 − m0 )ρN ∼ 1 would be demanded, so that α ∼ 1. Hencethe Zipf’s law follows. Basically it is a result of the hierarchical organization of chemical
reaction and the conservation law. Here, we note that this is a “self-consistency” argument,
and whether the hierarchical organization in the mean-ﬁeld approximation is unique in theoriginal equation is not answered.
This approximation would be oversimpliﬁed, since in reality, the concentration of eachchemical in a given layer is not identical, but diﬀers by each element, depending on the abundances of the chemical that catalyzes its synthesis. Then, we need to take the hierarchical1organization not only as the substrate but also as the catalyst into account. Although thisis beyond simple mean-ﬁeld approximation, it might be expected that the above power-law
abundances is also preserved within each layer, as the abundance of the catalyst follows
the power law according to the above argument. Thus, Zipf’s law across all chemicals is
expected, even though complete derivation for it would possibly require advanced techniquesuch as renormalization group.

2Supplementary Figures

nutrient uptake F0

10-2

10-3

10-4

10

-5

α=1
α=2
α=3
α=4

10-6
-6

10

10

-4

-2

0

10
10
nutrient concentration S

2

10

10

4

Fig. S 1: Adaptive dynamics of the nutrient uptake under diﬀerent values of parameter α. Thisparameter α represents the higher order eﬀect in the nutrient transportation, which is introduced
to investigate the eﬀect of the strength of negative feedback on the adaptive behavior. In the ﬁgure,scaled nutrient uptake F0 = S(N x1 )α is plotted for easy comparison of results with diﬀerent α. As
shown in the ﬁgure, when α ≥ 2, the nutrient uptake F0 is maintained close to the critical value,denoted by the dotted line, for wide ranges of the environmental nutrient concentration S. When
α = 1, the nutrient uptake decreases for most networks as S increases to a considerably high level,
where the attractor ceases to be a ﬁxed-point attractor. The concentrations of chemicals exhibit
periodic or chaotic oscillation, which reduces the (average) nutrient uptake from the optimal value.
This result suggested that the introducing a higher order eﬀect increases the parameter range
within which the self-regulation to the optimal growth is achieved.

3final concentration change

0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-0.8

-0.6
-0.4
-0.2
0
short-time concentration change

0.2

Fig. S 2: Relationship between the short-term response and the ﬁnal concentration change after
the environmental nutrient concentration S is changed. The other parameter values are identicalto those used in Fig.1. The short-term response is obtained as the concentration ratio of chemicals
between before (t < 0) and soon after changing S as S = 100 → 200 (t = 1500 in Fig. 3). The ﬁnal
concentration change is calculated as the ratio between before changing S and after the dynamics
stabilize to a ﬁxed-point in the new environment. The log10 -transformed concentration ratios areplotted. The solid line has a slope of 1 for reference. Recall that if perfect adaptation occurs,
the points lie on the horizontal axis, whereas if adaptation does not occur, the points lie on thereference line with a slope of 1. Here, the points lie between these two lines, implying partial
adaptation; these points are ﬁtted around a slope of 0.66, suggesting that all chemicals show acommon trend in partial adaptation.

4