A linear reciprocal relationship between robustness and plasticity in homeostatic biological networks

Abstract

In physics of living systems, a search for relationships of a few macroscopic variables that emerge from many microscopic elements is a central issue. We evolved gene regulatory networks so that the expression of core genes (partial system) is insensitive to environmental changes. Then, we found the expression levels of the remaining genes autonomously increase to provide a plastic (sensitive) response. A feedforward structure from the non-core to core genes evolved autonomously. Negative proportionality was observed between the average changes in core and non-core genes, reflecting reciprocity between the macroscopic robustness of homeostatic genes and plasticity of regulator genes. The proportion coefficient between those genes is represented by their number ratio, as in the “lever principle”, whereas the decrease in the ratio results in a transition from perfect to partial adaptation, in which only a portion of the core genes exhibits robustness against environmental changes. This reciprocity between robustness and plasticity was satisfied throughout the evolutionary course, imposing an evolutionary constraint. This result suggests a simple macroscopic law for the adaptation characteristic in evolved complex biological networks.

Introduction

In recent decades, robustness in biological systems has been studied extensively in systems and quantitative biology [1–3]. Robustness refers to the maintenance of certain features or functions of biological systems against noise or changes in the environment [4–6]. Mechanisms for the robustness of specific gene expression also have been intensively studied [7]. In particular, housekeeping gene expression levels are believed to be robustly maintained across different environmental conditions [8]. Recent measurements, however, have shown that not all housekeeping genes are robust against environmental changes, but are only partially robust [9, 10]. Thus, investigations are needed to determine the possible limitation in the degree of robustness across genes, the constraints on global gene expression changes, and how these mechanisms have evolved [11–14].

Homeostasis refers to the constancy of “macroscopic” physiological quantities against environmental changes, such as body temperature and blood glucose level. Although the mechanism of homeostasis has often been attributed to interactions among few organs, many “microscopic” dynamics also play a role, including neurotransmission and gene expression. Here, we investigated the emergence of “macroscopic” homeostasis in a biological network consisting of multiple “microscopic” elements, inspired by recent studies demonstrating the relation between macroscopic thermodynamic quantities and microscopic molecular dynamics [15, 16]. Of course, the biological system is not in equilibrium, and microscopic elements therein are heterogeneous. Nevertheless, the evolved biological states are stable. By virtue of this stability similar to the equilibrium state in thermodynamics, some macroscopic laws between robustness and plasticity can be uncovered, with insensitivity of the homeostatic component and changeability of the remnant part in response to external changes. Indeed, previous studies demonstrated a linear relationship between robustness in the period and plasticity in the phase in the circadian rhythm [17, 18].

In this study, we explored gene expression dynamics governed by evolving a gene-regulatory-network structure numerically so that the expression level changes of core genes were insensitive to environmental changes. Then, a feedforward structure from the non-core part of the network (i.e., the regulator genes) evolved autonomously. The robustness in the expression of core genes and the plasticity in that of regulator genes showed a linear reciprocal relationship; the increase of the former was associated with the latter. The proportion coefficient between those genes is represented by their number ratio, as in the “lever principle”, in which a decrease in the ratio results in a transition from perfect to partial adaptation, with only a portion of the genes exhibiting robustness. This result suggests a simple macroscopic law for the adaptation characteristic in evolved complex biological networks.

Model and methods

Gene regulatory networks are among the most well-known examples of complex biological networks [11–13, 19–21], in which each gene activates or inhibits other genes, including self-regulation (Fig 1). In the present model, these interactions are represented as an interaction matrix, J; when the jth gene activates or inhibits the ith gene, J_ij takes on a value of 1 or -1, respectively, and when there is no interaction, J_ij is 0. Inputs from the environment were further introduced to study adaptation against environmental changes, which globally regulate the expression of all genes. The environmental inputs are multi-dimensional, and each gene can exhibit a different degree of sensitivity to these inputs in two directions, represented as h_i = 1 or −1 for simplicity. Sign of h_i is randomly chosen and fixed during each evolutionary process. The environmental changes are represented by a single parameter, α. If signs of α and h_i are identical (different), the ith gene is activated (inhibited). We consider on/off-type gene expression dynamics with a given threshold: if the total input exceeds the threshold value, y_th, the genes are turned on. Further, the expression level of each gene can also fluctuate due to noise. These dynamics are given by the following stochastic differential equations:(1)where x_i is the expression level of the ith gene, β is the steepness of gene induction around the threshold (i.e., when β is sufficiently large, the gene expression dynamics approach on/off-type switching, and β corresponds to the Hill coefficient in a model that is often adapted in biology), and N is the total number of genes; the interaction term is scaled by considering the scale in random variables. ϵ is a small spontaneous induction, whose value does not change the result as long as it is much smaller than 1, and η_i(t) is the Gaussian white noise in the gene expression level. Of note, our model is quite similar to a neural network [20], and we can easily extend our results to other complex biological networks. Here, we set y_th to 0.3, β to 20.0, ϵ to 0.05, and σ² to 0.01.

[Figure omitted. See PDF.]

Each white circle represents a gene. Genes regulate the expression of other genes (including self-regulation). Triangular and flat arrowheads represent activating and inhibitory interactions, respectively.

We first investigated the adaptation dynamics involving a large number of components. In general, when the environment changes, organisms do not maintain all of the components constant but rather need to sustain only a portion of these essential components. Indeed, in adaptation experiments, the expression of only a portion of the genes in the network could be robustly maintained against environmental change, whereas the expression levels of most other genes were altered [22]. This is natural, because maintaining an entire system completely unchanged is impossible when each element is sensitive to the environment. To investigate the characteristics of these adaptation dynamics, we considered the following simple situation: some components behave as a core of homeostasis, while others function as regulators to maintain the core robustly against environmental changes. Accordingly, we designate genes incorporated in the core as core genes and the others as regulator genes (see Fig 1).

We then optimized the network structure J to achieve robustness of the homeostatic core by mimicking the evolutionary process [11, 12, 14, 23, 24]. From mutants with a slight change in J, we selected those exhibiting higher robustness in the expression of core genes to environmental changes for the next generation, as parameterized by α. First, the condition without an environmental stimulus was represented by α = 0. The system was then allowed to relax to a steady state to obtain the expression pattern , where is a steady-state value of x_i at α = α′. Steady-state value was calculated by averaging x_i over 1,000 time steps after the relaxation for each genotype at each generation. We then changed α to both positive and negative values (α₁ and −α₁) and let the system in each case relax to the steady-state again to obtain . Here, we set α₁ to 1.0. To analyze the robustness and plasticity of the gene expression levels, we calculated the average change in the expression level of genes in different environments as follows:(2) (3)

An individual with a smaller ΔX^C shows higher robustness in the core and is assumed to have higher fitness. Reproduction was asexual. Then, the kth individual with can produce an offspring with probability P_k, given as(4)where β_evo is the strength of the selection pressure. The networks J at the 0th generation are chosen randomly as described below. In each generation, each element in the offspring’s J is changed among {±1, 0} with probability p_mut. We set β_evo and p_mut to 40.0 and 0.01, respectively.

Note that here the steady state without an environmental stimulus did not change much by evolution, so that the evolution of insensitivity against environmental changes to be observed was not due to the change in the steady state. This is because we set the standard deviation of η and ϵ much smaller than the threshold. When α is zero, all genes are activated only by the noise given by η and spontaneous induction given by ϵ. Then, all genes are normally turned off and fluctuate around x_i = ϵ almost independently of J; since is 0.005, inputs from 60 genes (60% of the total genes) would have to concentrate on a given gene to exceed the threshold and to turn on a gene constantly without α. Such concentration hardly occurs due to the sparseness of networks. Hence, the average gene expression without input did not change through an evolutionary process, and all genes were kept to be off initially. In addition, since the threshold value was fixed, the sensitivity of genes against environmental changes was kept constant. Indeed, in Fig 2A (before evolution) and Fig 2B (after evolution), gene expression levels fluctuated around ϵ without input (time = 0–100), and those of the core transiently increased in response to changes in α (at time = 100 and 200), indicating that those genes were kept to be sensitive against environmental changes. This setting might seem biologically unrealistic, but we here focus on the evolution of responses against environmental changes and not on the nature of the steady-state gene expression pattern unlike previous studies (for example [25]). If we adopt constant selection for an arbitrary expression pattern, it will be difficult to distinguish the effect of the evolution of responses from that of the steady state. Thus, we adopt the present simplest settings, in which gene expression pattern is quite stable throughout evolution, and only responses can evolve.

[Figure omitted. See PDF.]

(A, B) Adaptation dynamics of genes (x_i(t)) of an individual with the highest fitness before (A: 0th generation) and after (B: 1000th generation) evolution. α was changed from 0 to 1 and from 1 to -1 at time 100 and 200, respectively. Black and gray lines indicate the time course of the core (N^C = 10) and regulator (N^R = 90) genes, respectively. (C) Changes in ΔX^C from the 0th to 1000th generations and (D) the corresponding trajectory at the ΔX^R–ΔX^C plane. All of the trajectories start from the same point (ΔX^C = ΔX^R = ΔX₀ ≃ 0.462). Different color lines indicate evolutionary trajectories with different N^C/N: magenta for 0.1, red for 0.2, orange for 0.3, yellow for 0.4, lime for 0.5, green for 0.6, cyan for 0.7, blue for 0.8, purple for 0.9, and brown for 1.0. Gray dotted and dashed lines are given by Eq 5 for N^C = 10 10 and 20, respectively.

In this study, we set N to 100 and the total number of individuals M to 300, unless otherwise noted. Initially, the elements in J take a value of 1 or -1 with a probability of p_link set to 0.1, and take 0 with a probability of p_link. We changed the fraction of the core genes to the whole genes, N^C/N, from 0.05 to 1.0 and investigated the dependence of the behavior of evolved gene expression dynamics on N^C/N. We repeated the simulation 100 times for each parameter set and plotted the average of those.

Results and discussion

For all N^C values, the network structures evolved to decrease ΔX^C (see Fig 2A and 2B as an example), whereas the final evolved state depended on N^C (Fig 2C). When N^C was sufficiently small, ΔX^C reached a steady value in the early generations, which was close to zero; that is, the core genes showed perfect adaptation [26, 27]. In contrast, when N^C was large, its steady-state value was larger than zero; that is, adaptation was only partial. This value increased with N^C. The system showed a transition from perfect to partial adaptation at N^C = N^C*, which lies between 20 and 25. It is noteworthy that the relaxation was slowed down for N^C/N = 0.2 to 0.6, in particular between N^C/N = 0.2 and 0.3. In statistical physics, it is commonly observed that the relaxation of systems to the equilibrium is slowed down in the vicinity of the transition point. Then, the observed slowing down will be the sign of critical slowing down and provide evidence of a transition between perfect and partial adaptation regimes. Note that even when N^C was equal to N (i.e., without the regulatory genes), ΔX^C still decreased slightly throughout evolution; that is, the networks can show intrinsic robustness without regulators. We define this ΔX^C value for the case of N^C = N as ΔX_int.

Interestingly, as ΔX^C decreased during evolution, ΔX^R increased almost monotonically in all cases (see Fig 2D). This result implies that under evolutionary selection, to increase the robustness of the expression of core genes, the plasticity in the expression of regulatory genes simultaneously increases. Here, the evolutionary trajectories in the space of ΔX^C and ΔX^R showed nearly linear behavior (Fig 2D). Evolution then stopped either when ΔX^C reached approximately zero or when ΔX^R increased and reached a certain threshold value, ΔX^R* (≃ 0.471). Of note this threshold value depends on the stochasticity in gene expression dynamics. If there is no fluctuation, the steady state value of each x_i is ϵ without environmental input, and 1 + ϵ under the input α with the same sign to h_i. Because half of the genes have positive h_i, and another half have a negative one, the threshold value would be estimated as 0.5 (= (1 + ϵ − ϵ)²/2.0). However, because of fluctuation and a constraint to keep gene expression positive, the steady state value of x_i without input will be higher than ϵ; x_i can reach about 0.3 by fluctuation but is truncated at 0 at the bottom, as shown in Fig 2A and 2B, where we set ϵ as 0.05. Then, the noise in gene expression is asymmetric in the vicinity of x_i = 0. Thus, the threshold value cannot reach 0.5 and will be slightly lower than 0.5. This suggests that there is an upper limit at which the regulator can buffer the changes in the core. If the buffering capacity is reached through evolution, the compensation by the regulators is not sufficient to allow for perfect adaptation of the core.

Here, values of ΔX^C (evolved from 0.462 to 0 at most) and ΔX^R (evolved from 0.462 to 0.473) were asymmetrical, which was due to the setting to make gene expression be zero without input and to be one only under sufficiently strong input. Since α₁ was set to be larger than the threshold y_th, half of the genes whose h₁ had the same sign as α were strongly activated. Then, at the beginning of evolution, when each gene regulatory network was random, only some, but not many, genes were suppressed if inhibitory regulations were concentrated on them. We set α₁ and y_th to 1 and 0.3, respectively, and then if seven inhibitory interactions were concentrated, core genes would be suppressed. Those genes’ expression levels would be intermediate between 0 and 1. Indeed, in Fig 2A, before evolution, some genes showed intermediate expression level, whereas those of almost all activated genes were around 1. Thus, when the environmental input was sufficiently strong, there was an asymmetry between the numbers of fully activated genes (whose expression level is almost 1) and intermediately activated genes (whose expression level is less than 1). When systems evolved to make robust core genes, inhibition of a large number of fully activated genes mainly contributed to ΔX^C. In contrast, activation of a small number of intermediately activated genes mainly contributed to ΔX^R. Therefore, the asymmetry between ΔX^C and ΔX^R emerged.

Evolved networks have distinct structures. The number of inhibitory interactions to the core genes from both the core and regulator increased (Fig 3A and 3B). In particular, when N^C was small, inhibitory interactions from the regulators prominently increased, whereas those from the core increased only slightly. By contrast, when N^C was large, the number of inhibitory interactions from the core also increased significantly. This suggests that regulation from the regulator is a primary driving force of homeostasis for small N^C, whereas for large N^C (i.e., small N^R), the core itself also functions in maintaining homeostasis. Indeed, even when all of the interactions from the regulators were removed from the evolved networks, ΔX^C still decreased when N^C was large (Fig 3E). ΔX^core started to decrease at around N^C ≃ 25–30. In the perfect adaptation region, the core genes were perfectly suppressed by regulatory genes, whereas if the regulation from the regulatory genes was lacking, the core genes never showed adaptation (below N^C ≃ 30). In contrast, in the partial adaptation region, the core genes regulated each other because the number of regulatory genes was insufficient. Then, these showed partial robustness without the regulation from the regulatory genes (above N^C ≃ 30). In contrast, the number of interactions from the core or regulator to the regulator changed only slightly compared with that of the initial random network (Fig 3C and 3D).

[Figure omitted. See PDF.]

(A–D) Difference of the linking probabilities between two nodes in the evolved networks from the default value p_link. RN indicates the random network. Each graph shows the linking probabilities (A) from the regulator to the core, (B) from the core to the core, (C) from the regulator to the regulator, and (D) from the core to the regulator. Red and cyan bars represent the linking probability for activating and inhibitory interactions, respectively. (E) ΔX of the core without every interaction from the regulator. (F) Flipping probabilities of each node from the off to on state or from the on to off state after a change in the sign of h_i. Each flipping probability is averaged for every node. Cyan circles and squares represent the flipping probabilities of nodes in the core and the regulator for a change in a node in the regulator, respectively. Red circles and squares represent these flipping probabilities for a change in a node in the core, respectively. The gray dotted line represents the flipping probability measured for the random network.

In the evolved networks, the propagation of perturbation also showed distinct changes from the random network. We analyzed how local perturbation to a gene propagates to the entire network; we changed the sign of h_i for a single gene and then counted the number of genes that were flipped between the on and off states. As shown in the cyan and red circles in Fig 3F, the flipping probability of each gene in the core increased with the increase in N^C, which reached the maximal level at around N^C ≃ 30. This will be related to the transition. If N^C increased, the contribution of each regulatory gene to repress genes in the core needed to be increased, and then perturbation of one regulatory gene had a great influence. However, if N^C was further increased above the transition point, the mutual regulation within the core genes turned to be relatively stronger than the regulation from the regulatory genes, then the influence of a flip of one regulatory gene decreased. Thus, we expect that the influence of perturbation will show the peak around the transition point from perfect to partial adaptation.

Note that the number of links to the regulator did not change (Fig 3C and 3D). Nevertheless, the flipping probability of each gene in the regulator decreased throughout evolution (see the cyan and red squares in Fig 3F). This result indicates that the expression of each regulator gene behaves more independently than those in the random network, which could increase the plasticity in the regulator genes, as shown in Fig 2D.

Then, we analyzed the quantitative relationship between ΔX^C and ΔX^R in the evolutionary steady state. The total change in gene expression in the core, N^C(ΔX^C − ΔX₀), and in the regulator, N^R(ΔX^R − ΔX₀), where ΔX₀ is ΔX of the random network, is plotted in Fig 4A. For N^C < N^C*, where perfect adaptation occurs, the following linear relationship was found:(5)with ΔX^C ≃ 0 due to perfect adaptation, where a is a positive constant (a ≃ 14.3). In contrast, for large N^C, only partial adaptation could be achieved, and ΔX^C remained finite but was still smaller than ΔX_int, the value for N^R = 0 (i.e., the case without the regulator). The difference ΔX^C − ΔX_int (< 0) for N^C is supported by the plasticity of the regulator; that is, the increment of ΔX^R from the random case. Indeed, for large N^C, we found the following linear relationship (see Fig 4B):(6)

[Figure omitted. See PDF.]

(A, B) Total change in gene expression in the core plotted against that in the regulator. Averaged values of ΔX^C and ΔX^R through 100 generations from the 900th to 1000th generation are used as the steady-state value. The difference of ΔX^C from ΔX₀ and from ΔX_int is plotted in (A) and (B), respectively. Green dotted lines are drawn to fit the points for N^C > 30 in (B). Two lines pass through the origin and have the same slope for both (A) and (B). (C, D) Lever principle for the robustness-plasticity relationship.

Again, to achieve the decrease in ΔX^C, ΔX^R changes more following the linear rule, where a takes on the same value as shown in Eq 5.

Interestingly, the linear relationship in Eq 5 was maintained throughout the evolutionary course. The time course of (ΔX^C, ΔX^R) satisfied Eq 5 as long as N^C < N^C* (see gray dotted and dashed lines in Fig 2D). This indicates that the linear relationship imposed a constraint at any evolutionary time point. Note that, in the case of large N^C, the evolutionary trajectory did not follow the linear relationship because the intrinsic robustness evolved; in Eq 6, ΔX^C approaches ΔX_int as N^R → 0, whereas the evolutionary trajectory started from ΔX₀. If the evolution of the regulations from the regulatory genes were much faster than the intrinsic robustness, the linear relationship would hold in the early stage. However, the timescale of evolution depends on N^C; if N^C is much larger than N^R, the number of intra-connections among the core genes is much larger than that of interconnections between core and regulatory genes and they can change with higher probability by the mutation in each generation. Then, intra-connections will typically change faster than the interconnections and evolve faster. This will be why we cannot observe the clear linear relationship for large N^C.

Therefore, we demonstrated linear relationships between robustness in the homeostatic core and plasticity in the regulator in evolved networks. Specifically, when the system shows higher robustness, it shows higher plasticity. If the fraction of the core is large, the homeostatic core will achieve only partial (i.e., not perfect) adaptation. Nevertheless, the linear relationship holds with the same proportion coefficient as in the perfect adaptation case.

Although the derivation of the linear relationships (Eqs 5 and 6) requires further study, an analogy with the lever principle may provide a more intuitive interpretation. For N^C < N^C*, all genes in the core show perfect adaptation and ΔX^C approaches ∼0, for which the total plasticity in regulator genes N^R(ΔX^R − ΔX₀) compensates for the original change in the core N^C Δ X₀ (Fig 4C). Then, if the number of plastic genes required to reach a balance exceeds N^R, the plasticity of the regulator genes is not sufficient to cancel out changes in the core, and adaptation is only partial. In the latter scenario, the intrinsic robustness conferred by regulation from the core evolved (Fig 3E), so that “the weight” for compensation is deducted, whereas the action by the regulator is maintained (Fig 4D). Therefore, the linear relationship with the same coefficient also holds for the case of partial adaptation.

The coefficient a in the lever rule is estimated by noting that at the transition point from perfect to partial adaptation (N^C = N^C*, N^R = N^R*), the regulator genes are fully plastic, whereas ΔX^C = 0 is maintained. Then,

Noting that N^C* ≃ 22, N^R* ≃ 78 according to Fig 2, and recalling ΔX^R* ≃ 0.471 and ΔX₀ ≃ 0.462, a is estimated as a ≃ 14.3, which agrees well with the observed value.

Interestingly, the lever rule (Eq 5) is also valid for the evolutionary process. Consistency between evolutionary trajectories and the dependence of N^C on the stationary state in Eq 5 indicates that evolutionary progress satisfies the balance between robustness and plasticity. Hence, the same macroscopic law governs both evolutionarily optimized states and their evolutionary trajectories.

Finally, we investigated the system-size dependence of the robustness by changing the total number of genes and numerically evolved the system (Fig 5). We found that the robustness of the core, ΔX^C − ΔX⁰, was scaled by (N^C/N)N^3/4; if only the ratio of N^C to N were the important control parameter, robustness should be scaled by (N^C/N)N⁰, or if only the absolute number of the core genes were important, it should be scaled by (N^C/N)N¹. The actual scaling found numerically is intermediate between these two cases; the maximum core size that shows perfect adaptation increases with the total number of genes, but not in proportion to it, but with the power 3/4. This suggests that the balance between core and regulatory genes depends not simply on the number of genes.

[Figure omitted. See PDF.]

Dependence of ΔX^C − ΔX⁰ on (A) N^C/N and (B) (N^C/N)N^3/4 for systems with various numbers of genes.

The lever principle may impose a fundamental constraint on homeostasis. Previous analyses of gene expression changes in response to environmental stress revealed that the expression levels of some genes change transiently and then return to the original level, whereas those of others change continuously [22], corresponding to the core and regulator of our model, respectively. Interestingly, experimental data suggest that the total change in gene expression in the steady state is proportional or correlated to its transient change, which is similar to the reciprocity between robustness and plasticity according to the lever rule uncovered with our model. Further studies will be required to reveal how the lever rule emerges and if the rule can be generalized to other homeostatic behaviors in biology.

Supporting information

S1 Data.

https://doi.org/10.1371/journal.pone.0277181.s001

(ZIP)

S2 Data.

https://doi.org/10.1371/journal.pone.0277181.s002

(ZIP)

Citation: Hatakeyama TS, Kaneko K (2023) A linear reciprocal relationship between robustness and plasticity in homeostatic biological networks. PLoS ONE 18(1): e0277181. https://doi.org/10.1371/journal.pone.0277181

About the Authors:

Tetsuhiro S. Hatakeyama

Roles: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

E-mail: [email protected]

Affiliation: Department of Basic Science, University of Tokyo, Meguro-ku, Tokyo, Japan

ORICD: https://orcid.org/0000-0003-2088-9153

Kunihiko Kaneko

Roles: Conceptualization, Funding acquisition, Supervision, Validation, Writing – review & editing

Affiliations: Department of Basic Science, University of Tokyo, Meguro-ku, Tokyo, Japan, Center for Complex Systems Biology, Universal Biology Institute, University of Tokyo, Komaba, Tokyo, Japan, The Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark

ORICD: https://orcid.org/0000-0001-6400-8587

References

1. Barkai N, Leibler S. Nature. 1997; 387(6636): 913–917. pmid:9202124

2. Hatakeyama TS, Kaneko K. Proc Natl Acad Sci USA. 2012; 109(21): 8109–8114. pmid:22566655

3. Young JT, Hatakeyama TS, Kaneko K. PLoS Comput Biol. 2017; 13(3): e1005434. pmid:28288155

4. de Visser JAG et al. Evolution. 2003; 57(9): 1959–1972. pmid:14575319

5. Wagner A. Robustness and evolvability in living systems. Princeton(NJ): Princeton University Press; 2005.

6. Alon U. An Introduction to Systems Biology: Design Principles of Biological Circuits. Florida: CRC press; 2006.

7. Masel J, Siegal ML. Trends Genet. 2009; 25(9): 395–403. pmid:19717203

8. Thellin O et al. J Biotechnol. 1999; 75(2-3): 291–295. pmid:10617337

9. Eisenberg E, Levanon EY. Trends Genet. 2013; 29(10): 569–574. pmid:23810203

10. Nicot N, Hausman JF, Hoffmann L, Evers D, J Exp Bot. 2005; 56(421): 2907–2914. pmid:16188960

11. Ciliberti S, Martin OC, Wagner A. PLoS Comput Biol. 2007; 3(2): e15. pmid:17274682

12. Kaneko K. PLoS ONE. 2007; 2(5): e434. pmid:17502916

13. Nagata S, Kikuchi M. PLoS Comput Biol. 2020; 16(6): e1007969. pmid:32598360

14. Ancel LW, Fontana W. J Exp Zool. 2000; 288(3): 242–283. pmid:11069142

15. Kaneko K. Life: An Introduction to Complex Systems Biology. Heidelberg and New York: Springer; 2006.

16. Kaneko K, Furusawa C. Annu Rev Biophy. 2018; 47: 273–290.

17. Hatakeyama TS, Kaneko K. Phys Rev Lett. 2015; 115(21): 218101. pmid:26636874

18. Hatakeyama TS, Kaneko K. Phys Rev E. 2017; 95(3): 030201(R). pmid:28415261

19. Glass L, Kauffman SA. J Theor Biol. 1973; 39(1): 103–129. pmid:4741704

20. Mjolsness E, Sharp DH, Reinitz J. J Theor Biol. 1991; 152(4): 429–453. pmid:1758194

21. Inoue M, Kaneko K. PLoS Comput Biol. 2013; 9(4): e1003001. pmid:23592959

22. Gasch AP et al. Mol Biol Cell. 2000; 11(12): 4241–4257. pmid:11102521

23. Dutta S, Eckmann JP, Libchaber A, Tlusty T. Proc Natl Acad Sci USA. 2018; 115(20): E4559–E4568. pmid:29712824

24. Rivoire O, Phys Rev E. 2019; 100(3): 032411. pmid:31640027

25. Siegal ML, Bergman A. Proc Natl Acad Sci USA. 2002; 99(16): 10528–10532. pmid:12082173

26. Koshland DE, Goldbeter A, Stock JB. Science. 1982; 220–225. pmid:7089556

27. Asakura S, Honda H. J Mol Biol. 1984; 176(3): 349–367. pmid:6748079