Stability of gene regulatory networks

Homeostasis of protein concentrations in cells is crucial for their proper functioning, and this requires concentrations (at their steady-state levels) to be stable to fluctuations. Since gene expression is regulated by proteins such as transcription factors (TFs), the full set of proteins within the cell constitutes a large system of interacting components. Here, we explore factors affecting the stability of this system by coupling the dynamics of mRNAs and protein concentrations in a growing cell. We find that it is possible for protein concentrations to become unstable if the regulation strengths or system size becomes too large, and that other global structural features of the networks can dramatically enhance the stability of the system. In particular, given the same number of proteins, TFs, number of interactions, and regulation strengths, a network that resembles a bipartite graph with a lower fraction of interactions that target TFs has a higher chance of being stable. By scrambling the $\textit{E. coli.}$ transcription network, we find that the randomized network with the same number of regulatory interactions is much more likely to be unstable than the real network. These findings suggest that constraints imposed by system stability could have played a role in shaping the existing regulatory network during the evolutionary process. We also find that contrary to what one might expect from random matrix theory and what has been argued in the literature, the degradation rate of mRNA does not affect whether the system is stable.


Introduction
Cells require different protein levels to survive in different external environments. The expression of these proteins within the cell are therefore highly regulated. An important regulatory mechanism involves transcription factors (TFs), which are themselves proteins that can either up or down regulate the transcription of mRNAs coding for other proteins by binding to enhancer or promoter regions of the regulated gene [1]. Despite the importance of maintaining desired protein concentrations within cells, factors affecting the stability of these concentrations to perturbations have received little attention.
One approach of studying the stability of such systems with a large number of interacting components was introduced by May in the 1970s in the context of complex ecological communities [2]. The idea is that in a n-species community, the dynamics of the abundances N i of each species may in general be described by a set of ordinary differential equations: for i = 1, 2, ..., n, with corresponding steady-state solution N ss i such that f i ( N ss ) = 0 ∀ i. The dynamics of small perturbations about this steady-state x i (t) = N i (t) − N ss i , when linearized about N ss i , has the form: where A is the Jacobian matrix with elements A ij = ∂fi ∂Nj ss . If all the eigenvalues of A have a negative real part, the system relaxes back to the steady-state upon perturbations and the steady-state is said to be stable; if any of the eigenvalues have a positive real part, the steady-state is unstable as the system will move away from it (exponentially fast) when infinitesimally perturbed. To construct A, one would need to precisely know the functions f i , which is often hard to obtain. May's approach was to model A as a random matrix with independent, identically distributed off-diagonal elements (with mean 0, standard deviation σ,

Results
The model Gene expression involves two major steps: transcription and translation (Fig.1a). Transcription is the process in which mRNA is synthesized by RNA polymerase using DNA as a template. The transcription rate of a gene i therefore depends on the number of RNA polymerases n and its effective gene copy number g i which takes into account both its copy number and how strongly RNA polymerase can bind to the promoter of that gene [11]. Due to the presence of TFs, g i ( c) can in general depend on the set of protein concentrations c (Fig.1a). We assume that multiple TFs acting on the same gene act independently, with their effects stacking multiplicatively. This allows for both "OR"-and "AND"-gate-like combinatorial effects [12], and can emerge from a thermodynamic model of TF binding (SI Section ??). Therefore, we adopt the following form for transcriptional regulation throughout the paper: where g i0 is the gene copy number of i if it were unregulated, and γ ij controls the type and strength of regulation i.e. how much gene expression of i changes in the presence of the TF j. In particular, γ ij > 0 if j up-regulates i and −1 ≤ γ ij < 0 if j down-regulates i. For each regulatory interaction, we assume that the fold-change Ω ij is drawn from a uniform distribution between 1 and Ω max , such that γ ij = Ω ij − 1 if γ ij > 0 (up-regulating) 1 Ωij − 1 if γ ij < 0 (down-regulating) (4) since this would allow g i (c j ) to increase (if j up-regulates i) or decrease (if j down-regulates i) by a factor of Ω ij in the limit of high c j . In SI section ??, we show that the main results do not depend on the particular distribution P (Ω) used. Motivated by experimental measurements of the relationship between transcription factor input and gene expression output showing a sigmoidal functional form of f ij (c j ) [13,14], we take it to be a Hill function with n ij > 0. Following Ref. [11], we assume a threshold number n c of RNA polymerases above which the gene copy number is limiting the transcription rate (Fig.1b). When this is the case, transcription rate is proportional to g i and is independent of n. If instead n < n c , it is the RNA polymerases that are limiting, in which case the different genes have to compete for the limited pool of RNA polymerases. The transcription rate of a gene i is then proportional to both n and the fraction of RNA polymerases working on that gene, the gene allocation fraction: Denoting the number of different genes by N , the dynamics of mRNA m i for i = 1, ...N can therefore be described by the following equation: where k m characterizes the transcription rate of a single RNA polymerase, τ m is the mRNA lifetime, and n s is the maximum number of RNA polymerases per gene. Similarly for the process of translation where ribosomes make proteins using mRNA as a template, the translation rate depends on the number of ribosomes r and the mRNA copy number m i . As for RNA polymerases, there is also a threshold number of ribosomes r c above which mRNA number is limiting and below which ribosomes are limiting (Fig.1c). The dynamics of protein numbers p i for i = 1, ..., N (with n = p N −1 corresponding to RNA polymerases and r = p N corresponding to ribosomes) are therefore given by: where k p characterizes the translation rate of a single ribosome, τ p is the protein lifetime, and r s is the number of ribosomes per mRNA when ribosomes are in excess. Depending on whether the RNA polymerases and ribosomes are limiting, there are 4 different cellular phases (Fig.1b, c). The regime where n ≥ n c and r ≥ r c (phase 3 of the model, where the production rate of mRNAs and proteins are proportional to gene and mRNA copy numbers respectively) has been widely studied [15][16][17], but has been shown to be inconsistent with experimental observations in wild type cells showing exponential growth of protein levels [18,19]. Instead, the regime where n < n c and r < r c (phase 1 of the model) is the one where wild type fission yeast [18] and mammalian cells appear to be in [19]. We therefore focus on this phase for the rest of the paper. Note, however, that the phase 3 regime has been experimentally observed in defective budding yeast and mammalian cells that are excessively large [20], whereas the regime where RNA polymerases are in excess (n ≥ n c ) while ribosomes are limiting (r < r c ) (phase 2 of the model) has been observed in mutant fission yeast [18]. We will address these two phases in the SI. The regime where n < n c and r ≥ r c (phase 4 of the model) is biologically unrealistic as ribosomes are typically more expensive to make compared to other proteins and hence having excess ribosomes while RNA polymerases are limited would be inefficient [21,22]. This regime is therefore not considered.
It will be convenient to consider the dynamics for the concentrations of mRNAs c mi = mi V and proteins c i = pi V . In bacteria [23,24] and mammalian cells [25], the volume of the cell V is approximately proportional to the total protein mass. Hence, we assume for simplicity that each protein has the same mass and Effective gene copy number:  set the cell density to be 1, such that V = i p i . The dynamics for concentrations in phase 1 are then given by: where c mT = i c mi is the total concentration of all mRNAs and 1 τ = 1 τm − 1 τp is the difference between mRNA and protein degradation rates (which can be positive or negative).

Effects of network features and topology on stability of the system
To study how properties of the transcriptional regulatory network affect the stability of the system, we first consider the regime where the lifetime of mRNAs is much shorter than that of proteins, which is typically true for wild-type cells [26]. In this limit of fast mRNA degradation, the relaxation dynamics of mRNA is much faster than that of proteins such that dcmi dt ≈ 0 at all times. Eliminating the fast process (by substituting the steady-state mRNA concentrations c mi = kmcn kpcr+ 1 τ φ i ( c) obtained from Eqn. 9 into Eqn.10), the dynamics of protein concentrations can be written as a set of N ODEs: The stability of the system therefore depends only on the eigenvalues of the N × N Jacobian matrix A = k p c ss r (M − I), where we define the interaction matrix with the steady-state protein concentrations given by c ss i = φ i ( c ss ) (from Eqn. 11). Denoting λ M as the eigenvalues of M, the system is stable as long as the maximal real part of these eigenvalues λ M,rmax is smaller than 1 (such that all eigenvalues of A have a negative real part). It is therefore useful to understand the structure of M by breaking it into two parts using Eqn. 6: where captures the direct interactions between proteins, while is a rank-1 matrix that captures the indirect interactions arising from competition for ribosomes. It can be shown that both the structure of M (Eqn.13) and the fact that stability only depends on M still hold in the other phases, despite the exact equations for protein dynamics being different (see SI Section ??). Therefore, even though the simulations in the rest of this section are carried out in phase 1, our findings and conclusions also apply to the other phases.

Stability of the system scales with
√ N for random regulatory networks.
We start by exploring the stability of 'fully random' regulatory networks, which we take to be our null model. Since the maximum eigenvalue of a random matrix depends on the standard deviation of its elements, we first carry out a naive estimate of how the elements of M scale with N . With g i (c) given by Eqn. 3, Biologically, TF concentrations are often comparable to the values of dissociation constants K d for DNA binding [26]. Therefore, since c j ∼ 1/N , we also choose K ij ∼ 1/N (Eqn.5), which would allow cells to maintain the full range of gene expression response. From Eqn. 5, this implies that f ij ∼ O(1) and ∂fij ∂cj ∼ N , and hence M 1 and M 2 also scale with N (Eqns. 14, 15). We therefore expect M ij ∼ O(1) (Eqn. 13), and hence (from RMT), for λ M,rmax to scale approximately as √ N for random interaction networks. λ M,rmax also increases with the strength of the interactions γ, implying that the system will become unstable either when N exceeds a critical number or the regulation strength becomes too high. However, this argument neglects correlations between the elements of M, which could potentially be relevant. In fact, we will see in the later sections that the structure of M (Eqn. 13) plays an important role in influencing the stability of the system. Therefore, to test if this scaling relation holds, we constructed networks of a specified interaction density ρ by randomly selecting ρN 2 interactions from the N (N − 1) possibilities (where we have assumed that ribosomes cannot act as TFs), and choose half of the interactions to be up-regulating with the remaining half being down-regulating.
By taking the ensemble average over the randomly drawn networks, we indeed recover the √ N scaling (Fig. 2a), which is also robust to the fraction of up-and down-regulatory interactions (see SI Section ??, Fig. ??a) and the distribution of fold-changes P (Ω) (see SI Section ??, Fig. ??). For sufficiently large N or Ω max , we can no longer find the fixed point of the system. Nevertheless, by simulating the dynamics, we find that for interaction networks of a given N and ρ, we get oscillatory, followed by chaotic behaviour as Ω max is increased (Fig. 2b). Similar phenomena have also been described and analyzed in models of neural networks [27] and ecological systems [28].
However, transcriptional regulatory networks are typically not random. Instead, they are enriched for distinct structural features such as the following motifs: feedforward loops (FFL), single input module (SIM) and dense overlapping regulons (DOR) which do not contain any loops besides autoregulatory ones [1,29]. In the next few subsections we therefore explore the effects of network topology on the system stability. Figure 2: Stability of random interaction networks. (a) For random interaction networks (red markers,'rand'), the maximal real part of the eigenvalues of the interaction matrix λ M,rmax scales with √ N . Surprisingly, for random directed acyclic networks (blue markers,'DAG'), λ M,rmax also scales approximately with √ N . In both of these cases, increasing the interaction strength from Ω max = 1.5 (circles) to Ω max = 2 (triangles) increases λ M,rmax . These results suggest that the system will become unstable (i.e. log 10 (λ M,rmax ) exceeds 0, indicated by the black dashed line) when N or Ω max becomes too large. Each data point is obtained from an average of 10 randomly drawn networks. Each random network is constructed by randomly selecting ρN 2 interactions from N (N − 1) possibilities, with half of the interactions chosen to be up-regulating and the remaining half to be down-regulating. Construction of DAGs is described in (c). For each regulatory interaction, fold change is chosen uniformly between 1 and Ω max . [Other parameters: ρ = 0.01, n = 1.] (b) When systems go out of stability, they exhibit oscillatory (left, Ω max = 20) followed by chaotic behavior (right, Ω max = 200) as interaction strengths are increased. [Other parameters: N = 200, ρ = 0.2, n = 1, fully random network] (c) Random directed acyclic networks are constructed by randomly drawing connections between proteins (red circles represent TFs, blue circles represent non-TFs). If a drawn connection creates a loop (e.g. the grey arrow with a cross on it), it is rejected.
Random directed acyclic networks can also be unstable.
Since transcription networks as a whole resemble directed acyclic graphs (DAGs) [1,29], we explore the stability of such networks.
In systems where the Jacobian matrix reflects the presence of direct interactions between components, the elements of the Jacobian matrix A ij is 0 if j does not influence or regulate i. In such cases, if there are no interaction loops involving 2 or more components (e.g. E regulates F which also regulates E), A can be written as a triangular matrix for such a DAG and the eigenvalues are the diagonal elements of the matrix i.e. the self-regulation loops. The system is therefore stable if there are no auto-activation among the components i.e. there are no positive elements along the diagonal of A.
In our case, the presence of indirect interactions captured by the additional M 2 matrix (Eqn. 13) implies that even if the regulation network is a DAG, the stability of the system is not determined solely by the selfregulation loops. Instead, we find that if we draw DAGs randomly (constructed by adding a connection only if the resultant network is still acyclic, Fig.2c), even if there are no self interactions, the largest eigenvalue still scales approximately with √ N , suggesting that it is still possible for such a network to go unstable. Nevertheless, there is a negative offset in λ M,rmax compared to the fully random case (Fig. 2a), implying that the lack of loops does help to stabilize the system.
Bipartite structure can maintain stability of large networks.
A commonly found motif in the Escherichia coli sensory transcription networks is the dense-overlapping regulons (DORs) which consist of a set of regulators that combinatorially control a set of output genes [1,29,30]. There are several of these DORs in E. coli, each with hundreds of output genes, and they appear to occur in a single layer i.e. there is no DOR at the output of another DOR. Such a structure can be thought of as a bipartite graph in which there are 2 types of nodes representing transcription factors (TFs) and nontranscription factors (non-TFs), and every directed edge go from a TF to a non-TF. Since such graphs do not contain any regulatory loops (and are therefore also DAGs), we expect them to be more stable than random networks. However, they are a specific subset of DAGs in which none of the TFs are themselves regulated. This is also a key difference between these networks and bipartite, mutualistic networks commonly studied in ecological models [9,10]. In this subsection, we investigate the stability of such networks.
To study this problem, we first group proteins into 2 categories: q TFs and N − q non-TFs, such that for any general network the components of the Jacobian matrix have the following structure: where T 1 (T 2 ) is a q × q matrix representing the direct (indirect) effect of TFs on TFs while R 1 (R 2 ) is a (N − q) × q matrix representing the direct (indirect) effect of TFs on non-TFs, with their elements defined previously (Eqn. [13][14][15]. The non-zero eigenvalues of M are therefore the eigenvalues of the sub-matrix Q with elements: When the network is sparse, each TF only regulates a small fraction of the total number of genes. Since c ss ∼ 1/N , the strength of indirect interactions are therefore typically much weaker than that of direct interactions (i.e. the non-zero elements of M 2 are much smaller in magnitude than that of M 1 , Eqns.14,15).
Bipartite Random DAG Ω = 1.5 Figure 3: Stability of bipartite networks. (a) When constructing a bipartite interaction network, we group the proteins into TFs (red circles) and non-TFs (blue circles), and only allow directed regulatory interactions to go from a TF to a non-TF. (b) For bipartite networks, there is a critical value of P neg (that is slightly greater than 0.5) below which λ M,rmax = 0 and above which λ M,rmax > 0. In the regime where λ M,rmax = 0 (which can be considered to be 'deeply' stable since it is furthest from the point λ M,rmax = 1 where the system becomes unstable), this value of λ M,rmax stays the same even when N (star markers to circles) or interaction strengths Ω max (star markers to squares) are increased. (c) When there is an equal fraction of up/down regulatory interactions P neg = 0.5, λ M,rmax is independent of both N and Ω max for bipartite networks. This is in contrast to fully random networks (red markers) and random DAGs (blue markers) where the system approaches the instability limit (λ M,rmax = 1) as N or Ω max (circles to triangles) is increased. This implies that a bipartite network structure can maintain and enhance the stability of the system as N or Ω max is increased. [Parameters: ρ = 0.01, n = 1, number of TFs q = 0.1N .] When constructing random bipartite networks, we only allow TFs to regulate non-TFs (Fig.3a), implying that T 1 = 0. The matrix Q therefore only consists of weak indirect interactions, and we expect the maximal eigenvalue to be smaller than that of random networks and DAGs. Moreover, since in this case Q is of rank-1, it has a unique real eigenvalue λ Q,b which can be shown to be (see SI Section ??): where ∂loggT ∂ci = N j=1 c j ∂loggj ∂ci as defined in Eqn.15 are the elements of the M 2 matrix (and therefore small when the interaction density is low). The maximum eigenvalue of the interaction matrix M is then given by λ M,b = max(λ Q,b , 0), since 0 is also an eigenvalue of M (see Eqns. 17,18).
This expression (Eqn.20) implies that unlike for fully random networks and random DAGs, the stability of bipartite networks can depend strongly on the ratio of up-and down-regulating interactions (see SI Section ??). In particular, there is a limit on the total strength of down-regulation (relative to that of upregulation) for the system to be stable. For example, if the majority of the interactions are up-regulating, λ Q,b should be negative and hence λ M,b must be 0. On the other hand, λ M,b must be positive when the fraction of down-regulations is sufficiently high. This tendency for inhibitory (activating) interactions to destabilize (stabilize) the system comes from the indirect effect that a regulator has on itself: a slight increase in the concentration of an inhibitor from its steady-state value will reduce the gene copy number and hence mRNA levels of the regulated gene. This in turn causes the mRNA fraction and hence protein concentrations of the inhibitor to increase further. This positive feedback also exists in the other phases, although its physical origin may be different (see SI Section ??, Fig. ??b).
Indeed, by numerically constructing multiple instances of a bipartite network and varying the fraction of inhibitory interactions P neg , we find that λ M,b = 0 when P neg is below a critical value that is approximately (but slightly greater than) 0.5 (Fig. 3b). Importantly, within this regime, the value of λ M,b = 0 is independent of both N and the strength of interactions Ω max (Fig. 3b, c). This suggests that such a bipartite network structure can help to maintain and enhance the stability of the system, especially for large N and Ω max .
Scrambling the interactions of E. coli transcriptional regulatory network can destabilize the system.
Real transcription networks, however, are not strictly bipartite graphs -there are autoregulatory elements as well as transcription factors that regulate other transcription factors. To investigate how relevant network stability is to biological networks, we obtained the E. coli transcriptional regulatory network from ref. [31]. The network consists of u = 5654 regulatory interactions (of which u p = 3187 are up-regulating), with q = 211 TFs regulating N = 2274 genes. We compared its stability with that of randomly constructed networks with the same N , density of interactions ρ = u N 2 ≈ 0.0011, and ratio of positive (activating) to negative (inhibitory) regulation.
We first explored two different ways of scrambling the original network: (1) randomly choosing u directed connections out of the N (N − 1) possible connections, and (2) fixing the number of TFs q and randomly choosing u directed connections out of qN possibilities. The second method of scrambling is motivated by the fact that q N and the stability of the system is governed solely by the q × q matrix Q representing how TFs affect TFs (Eqn. 19). For each drawn interaction network, we randomly choose u p of the interactions to be up-regulating (γ ij > 0) and the rest to be down-regulating (γ ij < 0). We draw the fold-change Ω ij of each regulatory interaction from a uniform distribution between 1 and Ω max = 1000. This choice of Ω max is motivated by the fact that TFs have been shown experimentally to change target protein levels by 100-1000 fold [13].
We find that with the real network, the system always converges to a stable fixed-point regardless of the regulation strengths (Fig. 4a). In contrast, for the randomly constructed networks (both with and without keeping q fixed), the probability of the system becoming unstable drastically increases when the interactions become too strong (Fig. 4a). This loss of a stable fixed point can give rise to either an oscillatory (Fig. 4b) or chaotic behaviour (Fig. 4c). This suggests that for typical regulation strengths and density, the interaction network cannot be random, and that certain structural features of real networks are important for stability.

Network stability depends on the density of TF-TF interactions
Since it is the maximal eigenvalue of the q × q sub-matrix Q (Eqn. 19) that determines the stability of the system, and direct regulatory interactions are typically stronger than the indirect background effects, we expect a higher density of direct interactions among TFs to destabilize the system. This suggests that what matters for stability is not only the number of TFs and the total number of regulatory interactions, but also the fraction of those interactions that target TFs.
We therefore analyzed the composition of regulatory interactions in the E. coli transcription network, and found that there are (i) u s = 134 self-regulations (of which 42 are activating), (ii) u t = 373 TF-other TF interactions (of which 201 are activating), and (iii) u n = u − u s − u t = 5148 TF-nonTF interactions (of which 2944 are activating) (Fig. 5a). In comparison, the scrambling method that maintained both the number of TFs and the total number of interactions gives a smaller number of self-interactions ( u s = 2.5) and a larger number of direct TF-other TF interactions ( u t = 522) than in the real network. To investigate if this could be the origin of the enhanced stability of the E. coli regulatory network, we tried another scrambling method with the composition of the interactions kept fixed. In particular, after setting the first q = 211 (out of N = 2274) proteins to be TFs, we randomly drew the numbers of interaction pairs within the 3 categories (self, TF-otherTF, and TF-nonTF) by choosing each TF and its target separately. The sign of the interactions are then randomly assigned while maintaining the fraction of positive/negative interactions within each of these categories. We find that this scrambling procedure, which fixes the composition of regulatory interactions (in addition to N , q and ρ), significantly increases the probability of the network having a stable fixed point (Fig. 4a).
Direct interactions among TFs can either be auto-regulatory loops or TFs regulating other TFs. We explored the effects of both of these factors, and found that assuming up-and down-regulations to be equally likely, a random network is almost always stable when the density of TF-other TF interactions ρ q = ut q(q−1) is sufficiently low (Fig.5b). Above this threshold value of ρ q , the probability of the system not exhibiting a stable steady-state increases with ρ q (Fig.5b). This effect is observed regardless of the number of self-interactions or whether u n is kept fixed (Fig.5b).
While this implies that systems with a small number of TF-TF interactions are almost always stable, it does not mean that having a high density of TF-TF interactions will necessarily lead to an unstable system. This can be seen from the fact the the probability of the system being stable does not drop sharply with ρ q (Fig. 5b) − there are still systems with a relative high density of TF-TF interactions that are still stable. This suggests that in the high ρ q regime, the details of the interactions become important. For such a network with a large number of TF-TF interactions to be stable, the type and strength of those interactions will need to be more fine-tuned.
This phenomenon that a small ρ q promotes stability is consistent with the stability of bipartite networks (ρ q = 0) and the fact that direct regulatory interactions are typically much stronger than the indirect background interactions. Nevertheless, since Q (which has contributions from both T 1 and T 2 , Eqn. 19) is not a sparse matrix even when ρ q is small, we do not expect the maximal eigenvalue λ M,rmax to scale with ρ q the way it does for a q × q random matrix with density ρ q . Indeed, we find numerically that the presence of T 2 can affect λ M,rmax (Fig. ??), suggesting that the indirect coupling between proteins can also play a role in influencing the stability of the system.

Effect of degradation rates on protein level stability
So far, we have been working in the limit of fast mRNA degradation, where the stability of the system is governed only by the interaction matrix M (Eqn. 12). In this regime, since M is independent of degradation rates 1/τ m and 1/τ p (see Eqns. 12, 6, 3), these do not affect whether the system is stable. The relaxation rates are also independent of τ m and τ p , with the relaxation rate in the absence of interactions given by (from Eqn. 11): However, it is not clear if this insensitivity (of both stability and relaxation rates) to τ m and τ p still holds outside of the τ m τ p regime. Within the framework of RMT, a more negative self-regulation term typically increases the relaxation rate and hence has a stabilizing effect [2]. Here, we ask if this is the case by investigating how mRNA and protein degradation rates affect the stability of the system and its relaxation timescale. In particular, can faster mRNA degradation rates help to stabilize a system that would otherwise be unstable if mRNAs degrade too slowly?
Values of mRNA and protein degradation rates do not affect whether the system is stable.
To investigate how the degradation rates of proteins and mRNAs affect the stability of the system when τ m is not too small, here we consider the full set of 2N equations (Eqns. 9, 10) and study how the eigenvalues of the (2N × 2N ) Jacobian matrix J varies with τ m and τ p .
To compare the relaxation rates of the full system with the protein relaxation rates when there are no interactions, we work with the transformed Jacobian matrix: For an arbitrary regulatory network with a corresponding interaction matrix M (Eqn.12), we find that the eigenvaluesλ ofJ are given by (see SI Section ??): where λ M are the eigenvalues of M as before, and ω is a dimensionless quantity given by: which reflects the difference between mRNA and protein degradation rates 1 τ = 1 τm − 1 τp . Since on average cell volume increases exponentially with rate (see Eqn. 8): On the other hand, when τ m τ p , the eigenvalue spectrum approximately resembles two circular regions, one corresponding to the dynamics of mRNAs and one for that of proteins. In this limit, increasing mRNA degradation rate only shifts the eigenvalues for the mRNA 'sector' to more negative values, leaving the maximal real part of the eigenvalues approximately unchanged, ω = 0.5 (green) vs ω = 0.4 (blue).
a growing cell has to satisfy the condition 1 τpkpφr < 1. Therefore, since τ m ≥ 0, we have ω ≥ −1. The expression forλ (Eqn. 23) therefore implies that the system is stable if and only if λ M,rmax ≤ 1, regardless of the value of τ m and τ p (Fig 6a). We find that despite differences in the details of the model, this conclusion still holds in the other phases (see SI Section ??).
Therefore, unlike what has been argued in the literature and what one might expect from RMT, changing mRNA nor protein degradation rates has no effect on whether the overall system is stable. If steady-state protein concentrations are unstable because λ M,rmax is too large (e.g. when interactions are too strong), increasing mRNA or protein degradation rates can never help to stabilize the system. Importantly, this finding also implies that our results for how structural features of the transcription network affects stability holds outside the regime of fast mRNA degradation, since stability only depends on M.
Increasing mRNA degradation rate can improve response times, but only up to some limit.
Besides system stability, another quantity of biological interest is the response time of the system to perturbations, which is especially relevant for cells experiencing changes in nutrient conditions [32,33]. Since this relaxation timescale is determined by the slowest eigenvalue of the Jacobian matrix, here we discuss how the maximal real part of the eigenvaluesλ rmax changes with τ .
The expression forλ (Eqn.23) implies that when the system is stable (λ M,rmax < 1), the rate at which the system relaxes to steady-state initially increases as ω increases from −1, but eventually plateau off − in the ω → ∞ limit (where τ m τ p ),λ → λ M − 1 (Eqn. 23, Fig. 6a). This implies that there is some benefit to having fast mRNA degradation in terms of response times, but once mRNA degrades much faster than proteins, further increasing mRNA degradation rate no longer affects the response time of the system. The eigenvalue spectrum in this τ m τ p limit appears to consist of two circular regions, one for the dynamics of mRNAs and the other for that of proteins (Fig. 6b), reminiscent of the RMT's circular law. Increasing τ m only shifts the eigenvalues corresponding to the mRNA sector and hence does not affectλ rmax . This is consistent with the fact that when τ m τ p , the dynamics of the overall system is governed only by the protein 'sector' (Eqn.11). Therefore, the slowest relaxation rate back to steady-state levels depends only on M and increasing mRNA degradation rate no longer improves the response time.

Discussion
In systems with a large number of interacting components, the question of stability is often an important one, as results from random matrix theory (RMT) predict instability when the system size N becomes too large or interactions become too strong. In the context of gene expression, transcriptional regulation is crucial for cells to adapt to different environmental conditions by changing their gene expression levels. It is therefore important for transcriptional regulatory networks (TRNs) to be able to accommodate a large number of regulatory interactions without the system going unstable. However, we find here that similar to the intuition provided by RMT, λ ∼ √ N for a fully random regulation network, suggesting that the system will go unstable as the number of genes exceeds a threshold. In fact, based on typical values for the density of actual regulatory networks and interaction strengths, we find that the system has a high probability of being unstable if the TRN is randomly constructed.
Besides the number of genes, and the density and strengths of interactions, there are other factors that can affect the stability of the system, one of which is the network topology. This aspect is particularly relevant in this system since TRNs are far from being random but instead consist of recurring motifs. While the properties of these specific motifs have been widely studied and shown to be important for specific functions such as adaptation, robustness, and fast response to environmental changes [1,29,30], how they contribute to the overall stability of the network is less clear. We find here that global structural features of the network, which are fundamentally shaped by many of these motifs, can play a huge role in determining the stability of the system. In particular, given the same number of proteins, TFs, interaction density and regulation strengths, a network that resembles a bipartite graph with a lower density of TF-otherTF interactions ρ q has a higher chance of being stable. The significance of ρ q fundamentally arises because of two main factors: (i) the eigenvalues of the Jacobian matrix and hence the stability of the system about its steady-state are governed only by the TF sector (i.e. how perturbations in TF concentrations affect TFs), and (ii) for a sparse regulatory network, the indirect background interactions arising from competition for ribosomes between different genes are typically much weaker than the direct regulatory interactions.
TRNs are also known to be scale-free, having a power-law out-degree distribution. This is consistent with the fact that most TFs only regulate a small number of genes, but there are TFs that regulate a very large number of genes ('master regulators'). Within a more abstract model of gene regulatory dynamics, the presence of these outgoing hubs has been shown to significantly increase the probability of the system reaching a stable target phenotype when the interaction strengths are allowed to vary while the network topology is kept fixed [34]. Here, we find that having a low ρ q can already significantly stabilize the system without the need to control the degree distributions. Nevertheless, having just a few master regulators may contribute to the network having a low ρ q if for instance most of the regulations on TFs are carried out by the master regulators (and non-master regulators predominantly regulate non-TFs).
Besides structural features of the network, another factor that could affect stability is the degradation rates of mRNA and proteins. Based on RMT, one may expect faster degradation to stabilize the system. This has in fact been argued to be the case [5,6]. However, by taking into account the dynamics of protein concentrations and how it couples to the dynamics of mRNA levels, we find that this is not the case. Instead, the stability of the system depends solely on the regulatory network and the strengths of those regulations − if the system is unstable, it will be unstable regardless of how fast mRNA or protein degrades. This highlights the importance of taking into account key aspects of the interactions (through the form of the dynamical equations) when analyzing the stability of large coupled systems, similar in spirit to studies of ecological models where explicitly considering interactions mediated through competition for nutrients can give drastically different results from assuming random pairwise interactions between species [7].
From an evolutionary perspective, there are many possible factors (such as the range of gene expression levels, environmental conditions, response time [32,33], level of unwanted crosstalk [35], etc.) that drive the addition or removal of regulatory connections. Our findings suggest that in addition to these considerations, another fundamental factor is the stability of the overall network. For example, there could be many ways of achieving a certain task such as allowing the cell to switch between two desired gene expression levels in two different nutrient conditions, but the only ones that can survive are those that also maintain the stability of the system. In other words, stability of the system may have played a role in shaping current existing regulatory networks through the evolutionary process. Our approach can therefore provide insights into the design and evolutionary constraints for a functional regulatory network, which may potentially be useful for guiding the construction of synthetic genetic circuits [36][37][38].

Stability of gene regulatory networks (Supporting Information)
1 Example of a thermodynamic model of RNA polymerases and transcription factors binding to DNA We consider the scenario where a gene has 1 promoter site and L regulatory sites, each corresponding to a binding site for a different transcription factor. Let q i = ci Ki = e −β( i−µi) be the binding affinities of each site i = 0, 1, ..., L, where c i is the concentration of the protein i, K i and i are respectively the dissociation constant and binding energy between protein i and site i, and µ i is the chemical potential of i. We choose the index i = 0 to represent binding of RNA polymerase to the promoter and the other indices represent TF binding to the corresponding regulatory site.
The state of the system is then given by σ with σ i = {0, 1} representing whether the i th binding site is occupied. We allow pairwise interactions between RNAP and each of the TFs, but neglect any pairwise interactions among the TFs, such that the free energy E of any state is given by: where w 0j ≥ 0 captures the strength and nature of the pairwise interaction between a bound RNAP and a bound TF j. Specifically, w > 1 indicates a positive interaction (with the TF up-regulating gene expression), w = 1 indicates no interaction, while w < 1 indicates a repulsive interaction. The limit where w 0j = 0 corresponds to the case where the TF is a steric inhibitor i.e. binding of j blocks RNAP from binding to the promoter. Denoting Z ON (Z OF F ) as the sum over the weights of all possible RNAP-bound 'ON' (RNAP-unbound 'OFF') configurations, the equilibrium probability P b of RNAP binding to the promoter is given by where P b0 = q0 1+q0 is the probability of RNA polymerase being bound to the promoter in the absence of any transcriptional regulation (L = 0), and F reg ( c), which is the regulatory function which captures the effect of TFs on the the binding of RNAP to the promoter, is given by: with the approximation taken in the limit of low RNAP concentrations q 0 L i=1 1+qi 1+w0iqi and q 0 1. This model is therefore an example of how a multiplicative form for F reg ( c) can arise, and serves as a motivation for our choice of regulatory function for the effective gene copy number (which we assume to be proportional to the probability of RNAP binding to promoter). Even though in this model the Hill coefficient is 1 for the effect of individual TFs, one could imagine higher Hill coefficients if there are cooperative effects in the binding of each TF to its binding site.
where k m and k p are constants characterizing the transcription and translation rates of a single RNA polymerase and ribosome respectively, φ i ( c) = gi( c) j gj ( c) is the gene allocation fraction (with g i ( c) being the effective copy number of gene i), c mT = i c mi and 1 τ = 1 τm − 1 τp , with τ m and τ p being the lifetimes of mRNA and proteins respectively.
The corresponding steady-state concentrations are given by: Since by definition i φ i = 1, the total steady-state mRNA concentration is We denote the total number of genes by N , and choose the index N − 1 to represent polymerases (the number of which we also denote as n) and the N th index to represent ribosomes (the number of which we also denote by r). The Jacobian of the full coupled mRNA-protein system is a 2N × 2N matrix where we have made the assumption that RNAPs and ribosomes cannot act as transcription factors.
Letλ be the eigenvalues ofJ with corresponding eigenvectors − This therefore provides a relation between each value ofλ and its corresponding eigenvalue of ∂φ ∂c , which we denote by λ M . Since λ M is independent of τ , we can find how τ affectsλ for any given λ M : where ω = 1 τ kpc * r .

Phase 2:
The regime where RNAPs are in excess and ribosomes are limiting (n ≥ n c , r < r c ) Whenever it is the gene copy numbers (instead of RNAPs) that are limiting (n ≥ n c ), the transcription rate is no longer proportional to the number of RNAPs, and hence it is the mRNA numbers m i rather than their concentrations that are kept at steady-state levels within the cell. We therefore analyze the the dynamics for m i and c i which in this case are given by: where m T = j m j is the total number of mRNAs, n s is the maximal number of RNA polymerases a single gene can accommodate, and the other variables are as defined previously. The corresponding steady-state mRNA and protein levels are: where we note that as before the steady-state protein concentrations are independent of the degradation lifetimes.
Following the same approach as in the previous section, we define the scaled Jacobian matrix (Eqn. 9), where now the elements ofJ are given byÃ = − 1 such that The eigenvaluesλ ofJ are hence given bỹ where ω 2 = 1 τmkpc ss r − 1, and as before λ M denote the eigenvalues of ∂φ ∂c , which are independent of τ . This equation is the same as that in phase 1 (Eqn. 13) with ω replaced by ω 2 .
Therefore, in both of these phases we get similar dependence of the stability of the system on degradation rates -the system is always stable as long as λ M < 1 and unstable if λ M > 1, regardless of the values of ω or ω 2 .
Phase 3: The regime where both RNAPs and ribosomes are in excess (n ≥ n c , r ≥ r c ) In this regime, the dynamics of m i and c i are given by: wherek m = k m n s , andk p = kprs V is dependent on cell volume which is linearly increasing over time. It is useful to define the growth rate per unit volume µ 3 , which is given by: At steady-state, and while these are constant over the whole cell cycle, the rate at which the system goes back to steady-state levels after a perturbation depends on its current volume at that point in time.
The Jacobian matrix of this system can again be written as J = A B C D , where A = − 1 τm I, D = −k 1 m T I, B ij =k m ∂gi ∂cj , and C ij =k p (δ ij − c i ). Unlike phases 1 and 2, here D depends on m T =k m g T τ m which is a function of τ m .
The eigenvalues λ of J can be found from with where K =k 0k1 g T . Unlike the other phases, here we choose not to scale J by the diagonal elements of D since we are investigating the effect of τ m on the eigenvalues and D itself depends on τ m . We therefore have where as before, λ M are the eigenvalues of the interaction matrix M = ∂φ ∂c . This implies that the system becomes marginally stable (λ → 0), for both τ m → 0 and τ m → ∞, i.e. in both these limits, even if the system is stable, it takes a long time for it to relax back to its steady-state when perturbed. This suggests that there is an intermediate regime of τ m for which the system responds fast to perturbations away from steady-state. This 'Goldilocks effect' arises because when τ m is large, the restoring force for mRNA numbers is small, while for small τ m , the restoring force for protein concentrations is small.
If we were to consider the scaled relaxation ratesλ = λ b0,3 , where b 0,3 = Kτ m is the relaxation rate for proteins when there are no transcriptional regulation (which is also the growth rate per unit volume µ 3 in the limit τ p → ∞), thenλ where This expression is the same as that in phases 1 and 2 (Eqns. 13,21), with ω and ω 2 now replaced by ω 3 . Therefore, as before, the system is always stable as long as λ M < 1 and unstable if λ M > 1, regardless of the value of ω 3 .

Effect of sign of regulatory interactions on stability
In this section, we investigate how the relative fraction of up-and down-regulatory interactions affect the maximal eigenvalue λ M,rmax of the interaction matrix.
We find that for random and DAG networks, the fraction of up-regulating interactions p up does not significantly affect λ M,rmax (Fig. S1a). However, for bipartite networks (which do not have any direct interactions between TFs), having only down-regulating interactions (p up = 0) increases λ M,rmax dramatically compared to the scenario of having p up = 0.5 (Fig. S1b). This is consistent with the tendency for inhibitory (activating) regulations to destabilize (stabilize) the system, which comes from the indirect effect that a regulator has on itself: a slight increase in the concentration of an inhibitor from its steady-state value will reduce the gene copy number and hence mRNA levels of the regulated gene. When ribosomes are limiting (phases 1 and 2), this causes the mRNA fraction and hence protein concentrations of the inhibitor to increase further. In phase 3, the reduction in mRNA levels of the regulated gene reduces the rate at which proteins are made. This slowing down of the increase in cell volume causes the inhibitor protein concentration to increase. This effect is much smaller in the case of random and DAG networks because their stability is dominated by the stronger, direct interactions among TFs. 4 Effect of distribution of fold-change Ω ij on stability In this section, we investigate the effect that the distribution of fold-changes Ω ij of the regulatory interactions has on the maximal eigenvalue λ M,rmax of the interaction matrix.
In the main text, all the simulations were carried out with Ω drawn from a uniform distribution. For any fixed value of Ω max , having P (Ω) ∼ 1 Ω (such that the logarithm of Ω is uniformly distributed [1,2]) would result in a lower Ω and a higher fraction of weaker interactions. Nevertheless, we find that the qualitative behavior of how λ M,rmax scales with N remains unchanged (Fig. S2).
where c i is the concentration of TF i, and a i = ∂loggT ∂ci . Since this is a rank-1 matrix, it only has one eigenvalue λ b with corresponding eigenvector v such that for all i = 1, 2, ..., q. This implies that Therefore and v i = c i .
6 Effect of density of TF-otherTF interactions ρ q on maximum eigenvalue λ M,r max Here, we investigate how the density of TF-otherTF interactions ρ q affects the maximal eigenvalue λ M,rmax of the interaction matrix. We find that without any auto-regulation loops, increasing ρ q increases λ M,rmax , which is consistent with our observation that the probability of the system going unstable increases when ρ q is too large (Fig. ??b).
These values of λ M,rmax can be higher than the maximum eigenvalue of the corresponding matrix Q 1 consisting only of the direct interactions i.e. Q 1,ij = c i ∂loggi ∂cj (Fig. S3), especially for small values of ρ q , suggesting that the indirect interactions can potentially play a role in affecting the stability of the system. In fact, in the limit where ρ q = 0 (i.e. bipartite network), stability is only determined by these indirect interactions. Nevertheless, these indirect interactions are much weaker than the direct interactions, which accounts for the stability of the system at low ρ q .