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

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.

In recent decades, robustness in biological systems has been studied extensively in systems and quantitative biology [1][2][3].Robustness refers to the maintenance of certain features or functions of biological systems against noise or changes in the environment [4][5][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][12][13][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 equivalent 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].* hatakeyama@complex.c.u-tokyo.ac.jp 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.

Core genes
Regulator genes Schematic representation of the gene regulatory network.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.
Gene regulatory networks are among the most wellknown examples of complex biological networks [11][12][13][19][20][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.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: 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 [22]), and N is the total number of genes; the interaction term is scaled by √ N 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, 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 σ 2 to 0.01.
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 [23].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,24,25].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 {x st i (0)}, where x st i (α ′ ) is a steady-state value of x i at α = α ′ .We then changed α to both positive and negative values (α 1 and −α 1 ) and let the system in each case relax to the steady-state again to obtain {x st i (±α 1 )}.Here, we set α 1 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: An individual with a smaller ∆X C has a more robust core and is assumed to have higher fitness.Then, the kth individual with ∆X C k can produce an offspring with probability P k , given as 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.In this study, we set N to 100 and the total number of individuals M to 300.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 1 − 2p 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 .For all N C values, the network structures evolved to decrease ∆X C (see Fig. 2A and B 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 steadystate 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.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 /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 and 20, respectively. 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).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.Evolved networks have distinct structures.The number of inhibitory interactions to the core genes from both the core and regulator increased (Figs. 3A and B).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).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 (Figs.3C and D).
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.
Note that the number of links to the regulator did not change (Fig. 3C and D).Nevertheless, the flipping prob-ability 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.
Relationships between robustness and plasticity.(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 ∆X0 and from ∆Xint is plotted in (A) and (B), respectively.(C, D) Lever principle for the robustness-plasticity relationship.
Finally, 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 0 ), and in the regulator, N R (∆X R − ∆X 0 ), where ∆X 0 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: 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. 4C): 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 [28].
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 0 ) compensates for the original change in the core N C ∆X 0 .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 evolutionoptimized states and their evolutionary trajectories.Thus, the lever principle imposes 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 [23], 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.

FIG. 2 .
FIG. 2. Evolutionary process of gene regulatory networks.(A,B) Adaptation dynamics of genes 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 = ∆X0 ≃ 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 and 20, respectively.

FIG. 3 .
FIG.3.Interactions between the core and regulator genes in evolved networks with varied N C .(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 hi.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.