A Linear Reciprocal Relationship Between Robustness and Plasticity in Homeostatic Biological Networks
aa r X i v : . [ phy s i c s . b i o - ph ] D ec A Linear Reciprocal Relationship Between Robustness and Plasticity in HomeostaticBiological Networks
Tetsuhiro S. Hatakeyama ∗ and Kunihiko Kaneko Department of Basic Science, The University of Tokyo. (Dated: December 17, 2020)In physics of living systems, a search for relationships of a few macroscopic variables that emergefrom many microscopic elements is a central issue. We evolved gene regulatory networks so that theexpression of target genes (partial system) is insensitive to environmental changes. Then, we foundthe expression levels of the remaining genes autonomously increase as a plastic response. Negativeproportionality was observed between the average changes in target and remnant genes, reflectingreciprocity between the macroscopic robustness of homeostatic genes and plasticity of regulatorgenes. This reciprocity follows the lever principle, which was satisfied throughout the evolutionarycourse, imposing an evolutionary constraint.
In recent decades, robustness in biological systems hasbeen studied extensively in systems and quantitative bi-ology [1–3]. Robustness refers to the maintenance of cer-tain features or functions of biological systems againstnoise or changes in the environment [4–6]. Mechanismsfor the robustness of specific gene expression also havebeen intensively studied [7]. In particular, housekeep-ing gene expression levels are believed to be robustlymaintained across different environmental conditions [8].Recent measurements, however, have shown that notall housekeeping genes are robust against environmen-tal changes, but are only partially robust [9, 10]. Thus,investigations are needed to determine the possible limi-tation in the degree of robustness across genes, the con-straints on global gene expression changes, and how thesemechanisms have evolved [11–14].Homeostasis refers to the constancy of “macroscopic”physiological quantities against environmental changes,such as body temperature and blood glucose level. Al-though the mechanism of homeostasis has often been at-tributed to interactions among few organs, many “mi-croscopic” dynamics also play a role, including neuro-transmission and gene expression. Here, we investigatedthe emergence of “macroscopic” homeostasis in a bio-logical network consisting of multiple “microscopic” ele-ments, inspired by recent studies demonstrating the rela-tion between macroscopic thermodynamic quantities andmicroscopic molecular dynamics [15, 16]. Of course, thebiological system is not in equilibrium, and microscopicelements therein are heterogeneous. Nevertheless, theevolved biological states are stable. By virtue of thisstability equivalent to the equilibrium state in thermo-dynamics, some macroscopic laws between robustnessand plasticity can be uncovered, with insensitivity of thehomeostatic component and changeability of the remnantpart in response to external changes. Indeed, previousstudies demonstrated a linear relationship between ro-bustness in the period and plasticity in the phase in thecircadian rhythm [17, 18]. ∗ [email protected] In this study, we explored gene expression dynamicsgoverned by evolving a gene-regulatory-network struc-ture numerically so that the expression level changesof core genes were insensitive to environmental changes.Then, a feedforward structure from the non-core partof the network (i.e., the “regulator genes”) evolved au-tonomously. The robustness in the expression of coregenes and the plasticity in that of regulator genes showeda linear reciprocal relationship; the increase of the formerwas associated with the latter. The proportion coefficientbetween those genes is represented by their number ra-tio, as in the “lever principle”, in which a decrease inthe ratio results in a transition from perfect to partialadaptation, with only a portion of the genes exhibitingrobustness. This result suggests a simple macroscopiclaw for the adaptation characteristic in evolved complexbiological networks.
Core genes Regulator genes
FIG. 1. Schematic representation of the gene regulatorynetwork. Each white circle represents a gene. Genes regu-late the expression of other genes (including self-regulation).Triangular and flat arrowheads represent activating and in-hibitory interactions, respectively.
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 othergenes, including self-regulation (Fig.1). In the presentmodel, these interactions are represented as an interac-tion matrix, J ; when the j th gene activates or inhibits the i th gene, J ij takes on a value of 1 or -1, respectively, andwhen there is no interaction, J ij is 0. Inputs from the en-vironment were further introduced to study adaptationagainst environmental changes, which globally regulatethe expression of all genes. The environmental inputs aremulti-dimensional, and each gene can exhibit a differentdegree of sensitivity to these inputs in two directions,represented as h i = 1 or − α .If signs of α and h i are identical (different), the i th geneis activated (inhibited). We consider on/off-type geneexpression dynamics with a given threshold: if the to-tal input exceeds the threshold value, y th , the genes areturned on. Further, the expression level of each gene canalso fluctuate due to noise. These dynamics are given bythe following stochastic differential equations: dx i dt = 11 + exp( − β ( y i − y th )) − x i + ǫ + η i ( t ) , (1) y i = X j J ij x j √ N + αh i ,< η ( t ) > = 0 , < η i ( t ) η j ( t ′ ) > = σ δ ( t − t ′ ) δ ij , where x i is the expression level of the i th gene, β is thesteepness of gene induction around the threshold (i.e.,when β is sufficiently large, the gene expression dynamicsapproach on/off-type switching [22]), and N is the totalnumber of genes; the interaction term is scaled by √ N considering the scale in random variables. ǫ is a smallspontaneous induction, whose value does not change theresult as long as it is much smaller than 1, and η i ( t ) isthe 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 complexbiological networks. Here, we set y th to 0.3, β to 20.0, ǫ to 0.05, and σ to 0.01.We first investigated the adaptation dynamics involv-ing a large number of components. In general, when theenvironment changes, organisms do not maintain all ofthe components constant but rather need to sustain onlya portion of these essential components. Indeed, in adap-tation experiments, the expression of only a portion ofthe genes in the network could be robustly maintainedagainst environmental change, whereas the expressionlevels of most other genes were altered [23]. This is nat-ural, because maintaining an entire system completelyunchanged is impossible when each element is sensitiveto the environment. To investigate the characteristics ofthese adaptation dynamics, we considered the followingsimple situation: some components behave as a core ofhomeostasis, while others function as regulators to main-tain the core robustly against environmental changes.Accordingly, we designate genes incorporated in the coreas core genes and the others as regulator genes (see Fig.1).We then optimized the network structure J to achieverobustness of the homeostatic core by mimicking the evo-lutionary process [11, 12, 14, 24, 25]. From mutants witha slight change in J , we selected those exhibiting higherrobustness in the expression of core genes to environmen-tal changes for the next generation, as parameterized by α . First, the condition without an environmental stimu-lus was represented by α = 0. The system was then al-lowed to relax to a steady state to obtain the expressionpattern { x st i (0) } , where x st i ( α ′ ) is a steady-state value of x i at α = α ′ . We then changed α to both positive andnegative values ( α and − α ) and let the system in eachcase relax to the steady-state again to obtain { x st i ( ± α ) } .Here, we set α to 1.0. To analyze the robustness andplasticity of the gene expression levels, we calculated theaverage change in the expression level of genes in differentenvironments as follows:∆ X C ( N C ) = X i ∈ core X α ∈{ α , − α } ( x st i ( α ) − x st i (0)) N C , (2)∆ X R ( N R ) = X i ∈ regulator X α ∈{ α , − α } ( x st i ( α ) − x st i (0)) N R . (3)An individual with a smaller ∆ X C has a more robustcore and is assumed to have higher fitness. Then, the k th individual with ∆ X C k can produce an offspring withprobability P k , given as P k = exp( − β evo ∆ X C k ) P l exp( − β evo ∆ X C l ) , (4)where β evo is the strength of the selection pressure. Thenetworks J at the 0th generation are chosen randomly asdescribed below. In each generation, each element in theoffspring’s J is changed among {± , } 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 numberof individuals M to 300. Initially, the elements in J takea value of 1 or -1 with a probability of p link set to 0.1,and take 0 with a probability of 1 − p link . We changedthe fraction of the core genes to the whole genes, N C /N ,from 0.05 to 1.0 and investigated the dependence of thebehavior of evolved gene expression dynamics on N C /N .For all N C values, the network structures evolved todecrease ∆ 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 asteady value in the early generations, which was close tozero; 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 wasonly partial. This value increased with N C . The systemshowed a transition from perfect to partial adaptation at N C = N C ∗ , which lies between 20 and 25. Note thateven when N C was equal to N (i.e., without the regu-latory genes), ∆ X C still decreased slightly throughoutevolution; that is, the networks can show intrinsic ro-bustness without regulators. We define this ∆ X C valuefor 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 x A B Δ X C Δ X R C D
FIG. 2. Evolutionary process of gene regulatory networks.(A, B) Adaptation dynamics of genes of an individual with thehighest fitness before (A: 0th generation) and after (B: 1000thgeneration) evolution. α was changed from 0 to 1 and from 1to -1 at time 100 and 200, respectively. Black and gray linesindicate the time course of the core ( N C = 10) and regulator( N R = 90) genes, respectively. (C) Changes in ∆ X C fromthe 0th to 1000th generations and (D) the corresponding tra-jectory at the ∆ X R –∆ X C plane. All of the trajectories startfrom the same point (∆ X C = ∆ X R = ∆ X ≃ . N C /N : magenta for 0.1, red for 0.2, orange for 0.3, yellow for0.4, lime for 0.5, green for 0.6, cyan for 0.7, blue for 0.8, pur-ple for 0.9, and brown for 1.0. Gray dotted and dashed linesare given by Eq. 5 for N C = 10 and 20, respectively. Fig. 2D). This result implies that under evolutionary se-lection, to increase the robustness of the expression ofcore genes, the plasticity in the expression of regulatorygenes simultaneously increases. Here, the evolutionarytrajectories in the space of ∆ X C and ∆ X R showed nearlylinear behavior (Fig. 2D). Evolution then stopped eitherwhen ∆ X C reached approximately zero or when ∆ X R increased and reached a certain threshold value, ∆ X R ∗ ( ≃ N C was small, inhibitory interactionsfrom the regulators prominently increased, whereas thosefrom the core increased only slightly. By contrast, when N C was large, the number of inhibitory interactions fromthe core also increased significantly. This suggests thatregulation from the regulator is a primary driving forceof homeostasis for small N C , whereas for large N C (i.e.,small N R ), the core itself also functions in maintaininghomeostasis. Indeed, even when all of the interactions Δ X c o r e N C / N E FBAC
ActivateInhibitRN 0.05 1.0-0.03 0.0 0.03 D i ff e r en c e f r o m p li n k N C / NRN 0.05 1.0-0.03 0.0 0.03 D i ff e r en c e f r o m p li n k ActivateInhibit N C / NRN 0.05 1.0 D -0.03 0.0 0.03 D i ff e r en c e f r o m p li n k C → CC → RR → CR → RRNCore to core: C → CRegulator to regulator: R → R ActivateInhibitCore to regulator: C → R 0 1.0N C / N F li pp i ng p r obab ili t y D i ff e r en c e f r o m p li n k Regulator to core: R → C FIG. 3. Interactions between the core and regulator genesin evolved networks with varied N C . (A–D) Difference of thelinking probabilities between two nodes in the evolved net-works from the default value p link . RN indicates the randomnetwork. Each graph shows the linking probabilities (A) fromthe regulator to the core, (B) from the core to the core, (C)from the regulator to the regulator, and (D) from the coreto the regulator. Red and cyan bars represent the linkingprobability for activating and inhibitory interactions, respec-tively. (E) ∆ X of the core without every interaction from theregulator. (F) Flipping probabilities of each node from theoff to on state or from the on to off state after a change inthe sign of h i . Each flipping probability is averaged for everynode. Cyan circles and squares represent the flipping proba-bilities of nodes in the core and the regulator for a change ina node in the regulator, respectively. Red circles and squaresrepresent these flipping probabilities for a change in a nodein the core, respectively. The gray dotted line represents theflipping probability measured for the random network. from the regulators were removed from the evolved net-works, ∆ X C still decreased when N C was large (Fig.3E). In contrast, the number of interactions from thecore or regulator to the regulator changed only slightlycompared with that of the initial random network (Figs.3C and D).In the evolved networks, the propagation of perturba-tion also showed distinct changes from the random net-work. We analyzed how local perturbation to a genepropagates to the entire network; we changed the signof h i for a single gene and then counted the number ofgenes that were flipped between the on and off states.As shown in the cyan and red circles in Fig. 3F, the flip-ping probability of each gene in the core increased withthe increase in N C , which reached the maximal level ataround N C ≃ -0.1 0.0 0.0 0.007( Δ X R - Δ X ) N R / N ( Δ X C - Δ X i n t ) N C / N -0.12 0.0 0.0 0.007( Δ X R - Δ X ) N R / N ( Δ X C - Δ X ) N C / N BA N C N R Robustness( Δ X C - Δ X ) Plasticity( Δ X R - Δ X ) N C N R Robustness( Δ X C - Δ X int ) Plasticity( Δ X R - Δ X ) Intrinsicrobustness
C D
FIG. 4. Relationships between robustness and plasticity.(A, B) Total change in gene expression in the core plottedagainst that in the regulator. Averaged values of ∆ X C and∆ X R through 100 generations from the 900th to 1000th gen-eration are used as the steady-state value. The difference of∆ X C from ∆ X and from ∆ X int is plotted in (A) and (B), re-spectively. (C, D) Lever principle for the robustness-plasticityrelationship. Finally, we analyzed the quantitative relationshipbetween ∆ X C and ∆ X R in the evolutionary steadystate. 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 plot-ted in Fig.4A. For N C < N C ∗ , where perfect adaptationoccurs, the following linear relationship was found: N C (∆ X C − ∆ X ) ≃ − aN R (∆ X R − ∆ X ) (5)with ∆ X C ≃ a is apositive constant ( a ≃ . N C ,only partial adaptation could be achieved, and ∆ X C re-mained finite but was still smaller than ∆ X int , the valuefor N R = 0 (i.e., the case without the regulator). Thedifference ∆ X C − ∆ X int ( <
0) for N C is supported by theplasticity of the regulator; that is, the increment of ∆ X R from the random case. Indeed, for large N C , we foundthe following linear relationship (see Fig. 4C): N C (∆ X C − ∆ X int ) ≃ − aN R (∆ X R − ∆ X ) . (6)Again, to achieve the decrease in ∆ X C , ∆ X R changesmore following the linear rule, where a takes on the samevalue as shown in Eq. 5.Interestingly, the linear relationship in Eq. 5 wasmaintained throughout the evolutionary course. Thetime 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 imposeda constraint at any evolutionary time point [28]. Therefore, we demonstrated linear relationships be-tween robustness in the homeostatic core and plasticityin the regulator in evolved networks. Specifically, whenthe system shows higher robustness, it shows higher plas-ticity. If the fraction of the core is large, the homeostaticcore will achieve only partial (i.e., not perfect) adapta-tion. Nevertheless, the linear relationship holds with thesame proportion coefficient as in the perfect adaptationcase.Although the derivation of the linear relationships(Eqs. 5 and 6) requires further study, an analogy withthe lever principle may provide a more intuitive inter-pretation. For N C < N C ∗ , all genes in the core showperfect adaptation and ∆ X C approaches ∼
0, for whichthe total plasticity in regulator genes N R (∆ X R − ∆ X )compensates for the original change in the core N C ∆ X .Then, if the number of plastic genes required to reach abalance exceeds N R , the plasticity of the regulator genesis not sufficient to cancel out changes in the core, andadaptation is only partial. In the latter scenario, the in-trinsic robustness conferred by regulation from the coreevolved (Fig. 3E), so that “the weight” for compensa-tion is deducted, whereas the action by the regulator ismaintained (Fig. 4D). Therefore, the linear relationshipwith the same coefficient also holds for the case of partialadaptation.The coefficient a in the lever rule is estimated by notingthat at the transition point from perfect to partial adap-tation ( N C = N C ∗ , N R = N R ∗ ), the regulator genes arefully plastic, whereas ∆ X C = 0 is maintained. Then, a = − N C ∗ N R ∗ (0 − ∆ X )(∆ X R ∗ − ∆ X )Noting that N C ∗ ≃ N R ∗ ≃
78 according to Fig. 2,and recalling ∆ X R ∗ ≃ .
471 and ∆ X ≃ . a is esti-mated as a ≃ .
3, which agrees well with the observedvalue.Interestingly, the lever rule (Eq. 5) is also valid for theevolutionary process. Consistency between evolutionarytrajectories and the dependence of N C on the station-ary state in Eq. 5 indicates that evolutionary progresssatisfies the balance between robustness and plasticity.Hence, the same macroscopic law governs both evolution-arily optimized states and their evolutionary trajectories.Thus, the lever principle imposes a fundamental con-straint on homeostasis. Previous analyses of gene ex-pression changes in response to environmental stress re-vealed that the expression levels of some genes changetransiently and then return to the original level, whereasthose of others change continuously [23], correspondingto the core and regulator of our model, respectively.Interestingly, experimental data suggest that the totalchange in gene expression in the steady state is propor-tional or correlated to its transient change, which is sim-ilar to the reciprocity between robustness and plasticityaccording to the lever rule uncovered with our model.Further studies will be required to reveal how the leverrule emerges and if the rule can be generalized to otherhomeostatic behaviors in biology. ACKNOWLEDGMENTS
This research was partially supported by a Grant-in-Aid for Scientific Research (A) 431 (20H00123) and Grant-in-Aid for Scientific Research on Innovative Ar-eas (17H06386) from the Ministry of Education, Culture,Sports, Science and Technology (MEXT) of Japan. [1] N. Barkai and S. Leibler, Nature, 387(6636), 913-917(1997).[2] T.S. Hatakeyama and K. Kaneko, Proc. Natl. Acad. Sci.USA., 109(21), 8109–8114 (2012).[3] J.T. Young, T.S. Hatakeyama, and K. Kaneko, PLoSComput. Biol., 13(3), e1005434 (2017).[4] J.A.G. de Visser et al.
Evolution, 57(9), 1959-1972(2003).[5] A. Wagner, Robustness and evolvability in living systems(Princeton University Press, Princeton NJ, 2005).[6] U. Alon, An Introduction to Systems Biology: DesignPrinciples of Biological Circuits (CRC press, Florida,2006).[7] J. Masel and M. L. Siegal, Trends Genet., 25(9), 395–403(2009).[8] O. Thellin et al.
J. Biotechnol., 75(2-3), 291–295 (1999).[9] E. Eisenberg and E.Y. Levanon, Trends Genet., 29(10),569–574 (2013).[10] N. Nicot, J.F. Hausman, L. Hoffmann, and D. Evers, JExp. Bot., 56(421), 2907–2914 (2005).[11] S. Ciliberti, O.C. Martin, and A. Wagner, PLoS Comput.Biol., 3(2), e15 (2007).[12] K. Kaneko, PLoS ONE 2(5), e434 (2007).[13] S. Nagata and M. Kikuchi, PLoS Comput. Biol., 16(6),e1007969 (2020).[14] L. W. Ancel and W. Fontana, J. Exp. Zool., 288(3), 242-283 (2000).[15] K. Kaneko, Life: An Introduction to Complex SystemsBiology (Springer, Heidelberg and New York, 2006).[16] K. Kaneko and C. Furusawa, Annu. Rev. Biophy., 47, 273–290 (2018).[17] T.S. Hatakeyama and K. Kaneko, Phys. Rev. Lett.,115(21), 218101 (2015).[18] T.S. Hatakeyama and K. Kaneko, Phys. Rev. E, 95(3),030201(R) (2017).[19] L. Glass and S.A. Kauffman, J. Theor. Biol., 39(1), 103–129 (1973).[20] E. Mjolsness, D.H. Sharp, and J. Reinitz, J. Theor. Biol.,152(4), 429–453 (1991).[21] M. Inoue and K. Kaneko, PLoS Comput. Biol., 9(4),e1003001 (2013).[22] β corresponds to the Hill coefficient in a model that isoften adapted in biology.[23] A.P. Gasch et al. , Mol. Biol. Cell, 11(12), 4241–4257(2000).[24] S. Dutta, J.P. Eckmann, A. Libchaber, and T. Tlusty,Proc. Natl. Acad. Sci. USA., 115(20), E4559-E4568(2018).[25] O. Rivoire, Phys. Rev. E, 100(3), 032411 (2019).[26] D.E. Koshland, A. Goldbeter, and J.B. Stock, Science,220–225 (1982).[27] S. Asakura and H. Honda, J. Mol. Biol., 176(3), 349–367(1984).[28] In the case of large N C , the evolutionary trajectory didnot follow the linear relationship because the intrinsicrobustness evolved; in Eq. 6, ∆ X C approaches ∆ X int as N R →
0, whereas the evolutionary trajectory startedfrom ∆ X0