Canalization reduces the nonlinearity of regulation in biological networks

Biological networks, such as gene regulatory networks, possess desirable properties. They are more robust and controllable than random networks. This motivates the search for structural and dynamical features that evolution has incorporated into biological networks. A recent meta-analysis of published, expert-curated Boolean biological network models has revealed several such features, often referred to as design principles. Among others, the biological networks are enriched for certain recurring network motifs, the dynamic update rules are more redundant, more biased, and more canalizing than expected, and the dynamics of biological networks are better approximable by linear and lower-order approximations than those of comparable random networks. Since most of these features are interrelated, it is paramount to disentangle cause and effect, that is, to understand which features evolution actively selects for, and thus truly constitute evolutionary design principles. Here, we compare published Boolean biological network models with different ensembles of null models and show that the abundance of canalization in biological networks can almost completely explain their recently postulated high approximability. Moreover, an analysis of random N–K Kauffman models reveals a strong dependence of approximability on the dynamical robustness of a network.


Introduction
Biological systems are frequently represented as networks, which describe the interactions between different biological entities such as genes, proteins or metabolites.For instance, gene regulatory networks (GRNs) describe how a collection of genes governs key processes within a cell.A static biological network is completely described by a wiring diagram, which contains nodes (e.g., genes) and edges between nodes, which can be undirected (e.g., in protein-protein interaction networks), directed and even signed (e.g., in gene regulatory networks).Static networks are, however, insufficient to obtain accurate insights into the often complex, non-linear dynamics of biological networks [1].Dynamic biological networks possess the additional information how each node is regulated by the set of its regulators.Popular dynamic modeling frameworks include differential equation models and discrete models.While the former harbors the potential for quantitative predictions, it requires a substantial amount of data for an accurate inference of its many kinetic parameters.Therefore, many modelers prefer discrete models and their qualitative predictions.Boolean networks constitute the simplest type of discrete model.Here, each node takes on only two values and time is discretized as well.The two values can be interpreted as low and high concentration, unexpressed and expressed genes or proteins, etc. Particularly for GRNs, Boolean networks have become increasingly popular.Over 160 Boolean GRN models have been curated by experts in their respective fields -most over the course of the last twelve years [2].These models range in size from 3 to 302 nodes, and describe various processes in many species and kingdoms of life.
Over the last few decades, a number of interesting features of biological networks have been identified.At the structural, "wiring diagram" level, biological networks are sparsely connected with an average degree of about 2.5, and are enriched for certain network motifs such as coherent feed-forward loops and complex feedback loops, particularly those that contain many negative interactions [2,3].Dynamically, most biological networks operate at the critical edge between order and chaos [2,4,5].For random N − K Kauffman networks, it is well-established that the network dynamics are generally ordered whenever 2Kp(1 − p) < 1 and chaotic whenever 2Kp(1 − p) > 1; at 2Kp(1 − p) = 1, a phase transition happens [6,7].Here, K is the average degree of the network, while p describes the bias of picking a one in the Boolean function's truth table; the unbiased case corresponds to p = 0.5, and the absolute bias can be quantified by 2|0.5 − p| ∈ [0, 1], or alternatively by 1 − 4p(1 − p) ∈ [0, 1], with 0 corresponding in both cases to the unbiased case.Networks with ordered dynamics typically possess few and short attractors, while chaotic dynamics are characterized by the presence of many long attractors [8].
The dynamic update rules of Boolean biological network models are also remarkable.They are highly canalizing, redundant and biased [2,9,10].Canalization is a widely used term in biology.First coined by developmental geneticist Waddington in the 1940s [11], it refers to the tendency of developmental processes to follow particular trajectories, despite internal and external perturbations [12].In other words, it refers to low variation in phenotypes despite potentially high variation in genotypes and the environment [13].Correspondingly, Kauffman introduced Boolean canalizing functions as suitable update rules to describe the gene regulatory logic [14].A canalizing function possesses a canalizing variable, which, when it receives its canalizing input, determines the output of the function, irrespective of all other inputs.If the subfunction which is evaluated when the canalizing variable does not receive its canalizing input is also canalizing, the function is 2-canalizing, etc [15].If all n variables of a function become eventually canalizing, the function is n-canalizing, also known as nested canalizing [16].The number of variables which become eventually canalizing is known as the canalizing depth [15].Every non-zero Boolean function possesses a unique standard monomial form, from which the canalizing depth and the number of variables in each "layer" of canalization can be directly derived [17,18].As the number of variables increases, canalizing and especially nested canalizing functions become increasingly rare [19][20][21].It is therefore very surprising that almost all rules in published Boolean biological network models are canalizing and even nested canalizing [2,9].
Another recently discovered feature of biological Boolean network models is the high approximability of their dynamics by linear and low-order continuous Taylor approximations of the Boolean update rules [22].Here, the mean approximation error (MAE) is defined as follows: Each update rule of a given Boolean network is replaced by a continuous Taylor approximation of a defined order.The MAE describes the mean squared error between the long-term state of the Boolean network and the long-term state of the continuous approximation when starting from a random initial state (see Methods for details).Manicka et al found that biological networks were consistently more approximable (i.e., had lower MAE values) than random networks with the same wiring diagram (i.e., matching degree distribution) and matching update rule bias [22].
Many of the described remarkable features of biological networks are interrelated and correlated.For instance, canalizing Boolean functions are on average more redundant and biased than random functions [2].In this paper, we show that the described increased approximability of biological networks can be fully explained by the abundance of canalization, which was not considered in [22].We further show that the approximability of a Boolean network depends mostly on its dynamic regime, which in turn depends on its update rules (that is, average degree, bias and amount of canalization) [2,23].A network with ordered dynamics (i.e., few and short attractors) tends to possess much more approximbale dynamics than a network with chaotic dynamics.

Results and discussion
To test the hypothesis that the increased canalization in biological networks explains their increased approximability, we compared the approximability of published expertcurated biological networks with several ensembles of random null models, similar to [22].All random networks possessed the same wiring diagram as the respective biological network.The authors in [22] considered an "unconstrained" null model, where each biological update rule was replaced by a non-constant random Boolean function (of the same degree), and a "constrained" model (null model type 1 in this study), which additionally matched the bias of each biological update rule.Neither model accounted for the high degree of canalization in biological networks.We therefore Fig. 1 Canalization explains the high approximability of biological networks..The distribution of mean approximation errors is shown for the biological networks (orange) and three different types of random null networks (shades of blue), which match different characteristics (bias and/or canalizing depth) of the biological network.Each box depicts the interquartile range (IQR), each whisker extends to the most extreme value within 1.5 * IQR from the box, and each horizontal line within a box depicts the median.For a fixed approximation order (1-3, x-axis), differences between the MAE distribution of the biological and the random networks are assessed using the two-sided Wilcoxon signed-rank test.Fig. 2 contains scatterplots showing the MAE values of all biological networks and their random null models.
considered two additional null models, one which matches the degree and canalizing depth of each biological update rule (null model type 2), and one which matches degree, canalizing depth and bias (null model type 3; see Methods for details).After excluding highly similar biological models and those with a maximal degree of eleven or more, we compared the approximability of 110 published expert-curated biological Boolean network models [2] and the three different ensembles of null models.As in [22], we found that random networks of type 1 were less approximable (Fig. 1).However, random networks that accounted for the increased canalization (null models of type 2 and type 3) exhibited similar levels of approximability as the biological networks.Interestingly, the higher the order of employed approximation the more significant were the differences in the MAE distributions between biological and random networks (Fig. 2).Third-order Taylor approximations recovered the dynamics of biological networks slightly better than those of random networks with matched degree, bias and canalizing depth.Overall, these results show that the approximability of biological networks can be almost entirely explained by their high degree of canalization, measured by the canalizing depth.
However, a related question, which has implications for the control of Boolean networks [24], remains: Why can the dynamics of biological networks be approximated so well by low-order and even linear continuous Taylor approximations?We hypothesized that the approximability of a Boolean network is strongly correlated with its dynamical robustness, which is typically measured by the average sensitivity [7] and Derrida values [25,26].These metrics describe how a small perturbation affects the network over time.If the perturbation gets on average smaller after each node has been synchronously updated once, the system operates in the ordered regime; if, on average, it increases in size, the system is in the chaotic regime, and if it remains, on average, of similar size, the system exhibits criticality.All biological systems that have thus far been modeled as Boolean networks operate close to the critical edge between order and chaos [2,4,5].This is likely because most update rules in biological networks are nested canalizing -in fact, biological networks are even particularly enriched for insensitive nested canalizing functions (NCFs) [2] -and the expected average sensitivity of an NCF in any number of variables is 1.On the contrary, the average sensitivity of random Boolean functions with degree k and bias p is 2kp(1 − p).That is, it increases as the number of inputs increases and decreases as the function becomes more biased (where p = 0.5 corresponds to the unbiased case).Boolean networks governed by such random functions thus exhibit a phase transition at 2kp(1 − p) = 1 [6,7].
To test which features of a biological network make it highly approximable, we computed Spearman correlations (ρ) between the mean approximation errors of the 110 biological networks and several dynamics-related properties (Fig. 3).Highly connected networks proved less approximable (ρ > 0.6).This is likely due to the fact that a continuous Taylor approximation of order n matches a Boolean function with k ≤ n variables perfectly everywhere.Thus, the higher the average degree, ⟨K⟩, of a Boolean network, the lower is the chance for perfect matches.Across the three approximation orders, the average degree was even slightly more negatively correlated with network approximability than the average effective degree, ⟨K e ⟩, defined in [10].This is somewhat surprising because the latter, which takes into account the importance of Boolean inputs, is a much stronger predictor of the dynamical robustness of a Boolean network, measured by its mean average sensitivity [2,23].In line with this, the strongest predictor of the mean average sensitivity of a Boolean network, ⟨K e ⟩⟨p(1 − p)⟩, as well as the mean average sensitivity itself were both not strongly correlated with the approximability of a Boolean network, with the correlation becoming insignificant for higher-order approximations.On the contrary, the proportion of Boolean rules in a biological network, which are nested canalizing, was fairly strongly correlated with the approximability, for all orders of approximation (|ρ| > 0.4).The higher this proportion the more approximable was the network.Canalizing rules, especially those with a low sensitivity, are typically fairly biased.In line with the result on the proportion of NCFs, more biased networks proved more approximable (|ρ| > 0.5).Biological Boolean rules with a higher number of inputs tend to be more biased [5].Interestingly, the covariance between p(1 − p) and the in-degree was the only property that became more correlated with approximability at higher approximation orders.
Metrics that explicitly describe dynamic aspects of a Boolean network also exhibited interesting correlations with the approximability.Assuming, as in the computation of approximability [22], a synchronous update of all nodes, we obtained, through simulation, for each biological network a lower bound of the number of attractors, as well as the approximate mean length of the attractors, the proportion of steady state attractors and the entropy of the basin sizes (see Methods).While the third-order approximability was not correlated with any of these metrics, networks with more attractors, a lower proportion of steady state attractors and higher entropy possessed dynamics that were less approximable at first and second order.This is surprising, since the presence of many long attractors, and concomitant high entropy, is associated with Boolean networks that operate in the chaotic regime [27].
To rule out potential confounders such as differences in network size, average degree as well as degree distribution, we considered modified N − K Kauffman networks, first defined in [28].In these random networks of size  degree K.The Boolean update rule of each node is generated by drawing 2 K times randomly with replacement from {0, 1} with probability 1 − p and bias p, respectively.We further required the wiring diagram of each network to be strongly connected since the dynamics decouple otherwise [29].Networks with a higher absolute bias exhibited more approximable dynamics (Fig. 4).Moreover, sparse networks (i.e., with low indegree) were on average more approximable.Interestingly, the MAE did not always decrease as the approximation order increased.For unbiased networks with high indegree (e.g., K = 5, p = 0.5), the MAE was very close to the maximally observed value of 0.25, even when using forth-order Taylor approximations.Since a Boolean function with K inputs is perfectly matched everywhere by a continuous Taylor approximation of order K, the MAE values were zero in these cases.If only J < K of the inputs of a Boolean function are essential (a Boolean input is non-essential if a change in this input never changes the output of the function.For example, f (x, y) = x has a non-essential input y), then the Jth order Taylor approximation already provides a perfect match.To rule out a potentially confounding effect created by perfect matches, we required, in a sensitivity analysis, all update rules to be non-degenerated, i.e., to contain only essential variables (Fig. S2).Most MAE For strongly connected 15-node Boolean networks with a constant in-degree (y-axis) governed by random update functions generated with a certain bias (x-axis), the mean error is shown when approximating their dynamics using different order Taylor polynomials (subplots).Each cell depicts the MAE across 50 networks, and the same networks were used to estimate the MAE using first-order to fourth-order Taylor polynomials.Results from an equivalent analysis where the functions are required to be essential in all its variables are shown in Fig. S2.
values were slightly higher, likely due to the higher effective degree.Qualitatively, the results were, however, very similar.Combining all 2000 random networks (100 each for combinations of constant indegree K ∈ {2, 3, 4, 5} and bias p ∈ {0.1, 0.2, 0.3, 0.4, 0.5}), we computed, as before, the Spearman correlation between MAE values and metrics that explicitly describe network dynamics.The dynamical robustness of a network, measured by the mean average sensitivity, was strongly positively correlated with first-, second-and third-order MAE values (ρ > 0.75; Fig. 5).Given that the average sensitivity of random Kauffman networks is 2Kp(1 − p) [7], this agrees qualitatively with the results from Fig. 4. Also in line is the finding that random networks are more approximable if they have few and short attractors, a high proportion of steady states, and low entropy in the distribution of the basin sizes.These four properties characterize networks that operate mostly in the ordered and critical dynamical regime.As observed for the biological networks, the correlations were consistently weaker when considering higher-order approximations.Note however that, by design of the computational experiment, 25% (50%) of the networks perfectly match their second-order (third-order) approximation, which certainly contributed to weaker correlations.To study the effect of canalization on the nonlinearity of regulation in more detail, we modified the random networks such that the update rules were restricted to specific classes of functions.First, we compared the approximability of random networks governed by 4-variable functions with different minimal canalizing depth (see Methods).While networks without required canalization were hardly approximable (MAE ≈ 0.25), the restriction to canalizing update rules gave rise to more approximable dynamics (Fig. 6a).Functions with a higher canalizing depth are however on average also less sensitive [30] and exhibit a higher absolute bias.
While the canalizing depth provides a crude measure of the amount of canalization in a Boolean function, more detailed information is contained in the canalizing layer structure [17,18,30].To investigate this, we compared the approximability of random networks, each governed entirely by 4-variable NCFs but with different layer structure.Networks governed by NCFs with layer structure k 1 = 4, e.g., an AND-NOT function x 1 ∧ x2 ∧ x3 ∧ x 4 , are highly approximable (Fig. 6b).On the other hand, networks governed by NCFs with layer structure k 1 = 1, k 2 = 3, e.g., functions such as x 1 ∨ (x 2 ∧ x 3 ∧ x 4 ), are much less approximable.Again, as the approximability of these networks decreases, the sensitivity of the underlying NCFs increases and the absolute bias decreases [30].

Conclusion
The idea of a probabilistic generalization of Boolean logic dates back all the way to George Boole [31].In this manuscript, we study in depth a recent implementation a b Fig. 6 The approximability of Boolean network dynamics depends on canalization.Each boxplot shows the distribution of the mean approximation error for 50 strongly connected N = 15-node Boolean networks with a fixed in-degree of K = 4 and a variable degree of canalization, characterized by (a) the minimal canalizing depth of each update rule (x-axis).In (b), all functions are nested canalizing (i.e., have canalizing depth 4) but the canalizing layer structure differs.The order of the Taylor polynomial used for the approximation is depicted by color.Fourth-order Taylor polynomials match the functions perfectly, the mean approximation error is 0. Each box extends across the interquartile range (IQR), whiskers extend to the lowest data point still within 1.5 IQR of the lower quartile, and the highest data point still within 1.5 IQR of the upper quartile, and black circles show outliers.
of this idea: using continuous Taylor approximations of Boolean functions to approximate the dynamics of a Boolean network.We show that the high approximability of biological networks, first postulated in [22], can be almost entirely explained by the abundance of canalization in biological networks.We conjecture that the remaining higher approximability of biological networks is due to the reported increased occurrence of insensitive canalizing rules in biological networks [2].Through a computational analysis of random networks, we show that the dynamical robustness of a network strongly influences its approximability: Networks with low mean average sensitivity, operating in the ordered and critical dynamical regime and characterized by few and short attractors, possess generally more approximable dynamics.In line with this, networks governed by canalizing or even nested canalizing functions that are highly biased and insensitive to perturbations proved more approximable.
Fully disentangling the relative contribution of the related properties canalization, bias, and sensitivity on approximability constitutes one of several open questions.Moreover, it remains to be investigated how well non-perfect continuous approximations of Boolean networks perform in the context of predicting control targets or specific dynamical features.A more technical question is whether Boolean functions that can be well approximated by low-order continuous extensions give rise to more approximable Boolean networks.base versions in the analysis.We approximated the entropy of the basin sizes as Finally, we used s as the lower bound of the number of attractors.In a network with many attractors, we almost certainly fail to discover all attractors when starting from only 1000 random states.However, all attractors with a large basin size are discovered with high probability.
For the random Boolean networks of fixed size n = 15, analyzed in Figs. 4, 5 6, S2, we computed the entire state space.All dynamics-related metrics, including the number of network attractors, are therefore exact in these analyses.

Taylor polynomials of Boolean functions
Since f is a continuous-variable function, we can consider different orders of approximation for f using its Taylor expansion.As described in [22], f is a square-free polynomial and its Taylor expansion is finite.More specifically, the n th order approximation will match f perfectly, and if only m < n inputs of f are essential, then the m th order approximation already matches f perfectly.
For a given α = (α 1 , . . ., α n ) ∈ {0, 1} n and x ∈ [0, 1] n , we define with the convention that ∂ 0 i f ≡ f .For p ∈ [0, 1] n , we have If p = ( 1 2 , . . ., 1  2 ), which represents the unbiased selection of each variable, then f (p) equals the output bias of f , as shown in [22].The Taylor decomposition yields different approximations of a Boolean function by restricting the sum in Equation 1to α with |α| ≤ m ≤ n.The Taylor polynomial of order m is given by Approximability of a Boolean network by continuous extensions 1} n be a Boolean network.We define the m th order approximation of F to be )), . . ., max(0, min (1, where the update functions of F (m) are the m th order Taylor approximations of the update functions of F , f i as defined in Equation 2, rescaled to the interval [0, 1].With this, we can define the mean approximation error (MAE) as the mean squared error between the long-term state of the Boolean network and the long-term state of its continuous approximation.That is, where F ∞ (x 0 ) and F (m),∞ (x 0 ) describe the long-term state of the Boolean network F and its m th order approximation, respectively.In practice, we approximated the MAE, using the Python library boolion [22], by updating both F and F (m) synchronously 25 times and using 1000 random initial values.

Canalization
This study employs several mathematical concepts related to canalization.By [14], a Boolean function f (x 1 , . . ., x n ) : {0, 1} n → {0, 1} is canalizing if there exists a canalizing variable x i , a canalizing input a ∈ {0, 1} and a canalized output b ∈ {0, 1} such that If the subfunction g is also canalizing, then f is 2-canalizing, etc.More generally, f is k-canalizing, where 1 ≤ k ≤ n, with respect to the permutation σ ∈ S n , inputs a 1 , . . ., a k , and outputs b 1 , . . ., b k if Here, f C = f C (x σ(k+1) , . . ., x σ(n) ) is the core function, a Boolean function on n − k variables.When f C is not canalizing, then the integer k is the canalizing depth of f [15].If k = n (i.e., if all variables are become eventually canalizing), then f is a nested canalizing function (NCF) [16].By [17], every nonzero Boolean function f (x 1 , . . ., x n ) can be uniquely written as where each M i = ki j=1 (x ij + a ij ) is a non-constant extended monomial, p C is the core polynomial of f , and k = r i=1 k i is the canalizing depth.Each x i appears in exactly one of {M 1 , . . ., M r , p C }.The layer structure of f is the vector (k 1 , k 2 , . . ., k r ) and describes the number of variables in each layer M i [18,30].

Random null models of biological networks
We compared biological Boolean network models to three ensembles of null models that matched different characteristics of the biological networks, as shown in Fig. 1.All null models matched the in-degree of the biological networks.Null models 1 and 3 matched, in addition, the bias of each biological update rule, while null models 2 and 3 matched the canalizing depth.
Let F = (f 1 , . . ., f n ) be a biological Boolean network model.For each f i , we first simplified the function to only include essential variables, yielding fi : {0, 1} k → {0, 1}, where k is the number of essential variables, i.e., the in-degree.While this step was omitted in [22], it appears important for an unbiased comparison given that close to 2% of regulators in biological networks are non-essential [2].We then computed the number of ones in the truth table of fi , denoted q, and fi 's canalizing depth d, following [18].
To obtain a random Boolean function g (for null model 1) with the same bias as fi and arbitrary canalizing depth, we simply selected a random subset Ω ⊆ {0, 1} k of size |Ω| = q, and set To obtain a random Boolean function g (for null model 2) with exact canalizing depth d and arbitrary bias, we randomly selected d out of fi 's k essential variables, arranged them in a random order, and randomly selected for each of the d variables a canalizing input value a ∈ {0, 1} and a canalized output value b ∈ {0, 1} (see Equations 4,5).Finally, we randomly selected a core function g C : {0, 1} k−d → {0, 1}, ensuring that g C depends on all k − d variables and that g C is not canalizing, by repeating this random selection process until both conditions were met.We then filled the truth table of g, as outlined in Equation 5.This entire procedure has already been implemented in the Python library canalizing function toolbox, published along with [2].
To obtain a random Boolean function g (for null model 3) with the same bias as fi and the same canalizing depth d, we followed the same procedure as for null model 2, with two exceptions.First, we did not randomly select the canalized output values b 1 , . . ., b d but instead used the canalized output values of fi .Otherwise, it is impossible to obtain the same bias.Second, we randomly selected a core function g C of g that has the same number of ones as the core function of fi (following the same approach as for null model 1).

Random Boolean networks
To generate a random Boolean network F = (f 1 , . . ., f N ) (modified N − K Kauffman network), we first generated a random directed graph of N nodes (the wiring diagram), where each node has K incoming edges.We ensured the graph is simple (i.e., does not contain self-edges/auto-regulations).We further ensured the graph is strongly connected since the dynamics decouple otherwise [29].
To obtain the random Boolean update rules f 1 , . . ., f N , we randomly selected, for the networks analyzed in Figs. 4, 5, any Boolean function g : {0, 1} K → {0, 1}.In a sensitivity analysis, reported in Fig. S2, we ensured that g is non-degenerated, i.e., that all variables of g are essential, by repeating the random selection until this condition was met.For the random Boolean networks with constant degree 4 and minimal canalizing depth d ∈ {0, 1, 2, 4}, analyzed in Fig. 6a, we followed a very similar procedure as for null model 2 (see above), with one exception: We allowed the core function to be canalizing so that the realized canalizing depth may be larger than d.For the random nested canalizing Boolean networks with constant degree 4 and different layer structure, analyzed in Fig. 6b, we followed again a very similar procedure as for null model 2, with the exception that the layer structure determines the canalized output values, b 1 , . . ., b 4 [30].For strongly connected 15-node Boolean networks with a constant in-degree (y-axis) governed by random non-degenerated update functions generated with a certain bias (x-axis), the mean error is shown when approximating their dynamics using different order Taylor polynomials (subplots).Each cell depicts the mean approximation error across 50 networks, and the same networks were used to estimate the mean approximation error using first-order to fourth-order Taylor polynomials.Results from an equivalent analysis where the functions are allowed to contain non-essential variables are shown in Fig. 4.

Fig. 2
Fig.2Mean approximation errors of biological networks and their random null models.For a fixed approximation order (1-3, columns), differences between the MAE values of the 110 biological and the different random networks (rows) are shown, in addition to the Spearman correlation coefficient, ρ.A summary of this data is shown in Fig.1.

Fig. 3
Fig. 3 Predictors of approximability of biological networks.Pairwise Spearman correlation between the first-, second-and third-order mean approximation errors and various network properties across the 110 published biological networks, ordered by the mean correlation.< • > denotes the mean, p = output bias, K = number of variables, Ke = effective connectivity, Cov = covariance of p(1 − p) and K.The pairwise Spearman correlations between all shown properties are in Fig. S1.

Fig. 4
Fig.4Effect of bias and in-degree on the approximability of the dynamics of Boolean networks.For strongly connected 15-node Boolean networks with a constant in-degree (y-axis) governed by random update functions generated with a certain bias (x-axis), the mean error is shown when approximating their dynamics using different order Taylor polynomials (subplots).Each cell depicts the MAE across 50 networks, and the same networks were used to estimate the MAE using first-order to fourth-order Taylor polynomials.Results from an equivalent analysis where the functions are required to be essential in all its variables are shown in Fig.S2.

Fig. S2
Fig.S2Effect of bias and in-degree on the approximability of the dynamics of nondegenerated Boolean networks.For strongly connected 15-node Boolean networks with a constant in-degree (y-axis) governed by random non-degenerated update functions generated with a certain bias (x-axis), the mean error is shown when approximating their dynamics using different order Taylor polynomials (subplots).Each cell depicts the mean approximation error across 50 networks, and the same networks were used to estimate the mean approximation error using first-order to fourth-order Taylor polynomials.Results from an equivalent analysis where the functions are allowed to contain non-essential variables are shown in Fig.4.
N , each node has constant