# The Dynamics of Canalizing Boolean Networks

Authors: Elijah Paul; Gleb Pogudin; William Qin; Reinhard Laubenbacher
Journal: Complexity
 (2020)

---

## Abstract

Boolean networks are a popular modeling framework in computational biology to capture the dynamics of molecular networks, such as gene regulatory networks. It has been observed that many published models of such networks are defined by regulatory rules driving the dynamics that have certain so-called canalizing properties. In this paper, we investigate the dynamics of a random Boolean network with such properties using analytical methods and simulations. From our simulations, we observe that Boolean networks with higher canalizing depth have generally fewer attractors, the attractors are smaller, and the basins are larger, with implications for the stability and robustness of the models. These properties are relevant to many biological applications. Moreover, our results show that, from the standpoint of the attractor structure, high canalizing depth, compared to relatively small positive canalizing depth, has a very modest impact on dynamics. Motivated by these observations, we conduct mathematical study of the attractor structure of a random Boolean network of canalizing depth one (i.e., the smallest positive depth). For every positive integerℓ, we give an explicit formula for the limit of the expected number of attractors of length ℓ in an n-state random Boolean network as n goes to infinity.


---

## Body

## 1. Introduction

Dynamic mathematical models are a key enabling technology in systems biology. Depending on the system to be modeled, the data and information available for their construction, the questions to be answered, and different modeling frameworks can be used. For kinetic models, systems of ordinary differential equations have a long tradition. Generally, they will have the very special structure of polynomial equations representing Michaelis–Menten kinetics, even in the case of systems, such as gene regulatory networks, that are not proper biochemical reaction networks. It is this special structure that gives models desirable properties and aids in model analysis. Besides continuous models, a range of discrete models are finding increasingly frequent use, in particular Boolean network models of a broad variety of biological systems, from intracellular molecular networks to population-level compartmental models (see e.g., [1–5]), going back to the work of Kauffman in the 1960s [6–8]. While Boolean network models, a collection of nodes, whose regulation by other nodes is described via a logical rule built from Boolean operators, are intuitive and mathematically simple to describe, their analysis is severely limited by the lack of mathematical tools. It generally consists of simulation results. Any set function on binary strings that takes on binary values can be represented as a Boolean function, so that the class of general Boolean networks is identical to the class of set functions on binary strings of a given length, making any general analysis impossible. The search for special classes of Boolean functions that are broad enough to cover all or most rules that occur in biology, but special enough to allow for mathematical approaches has a long history.It was again Kauffman who proposed a class of functions [7] with properties inspired by the developmental biology concept of canalization, going back to Waddington in the 1940s [9]. There is some evidence that canalizing Boolean functions do indeed appear disproportionately in published models and that the dynamics of Boolean network models consisting of canalizing functions has special properties, in particular a “small” number of attractors. This is important since, in the case of intracellular molecular network models, attractors correspond to the different phenotypes a cell is capable of. Here, again, the majority of available results are obtained by simulating large numbers of such networks. The main question of this paper is as follows: What do the dynamics of a random canalizing Boolean network look like?We approach this question using both computer simulations and analytical methods, with the main result of the paper being Theorem 2, which gives a provable formula for the number of expected attractors of a general Boolean network with a particular canalization property. In addition to providing important information about canalizing Boolean network models, this result can be viewed as a part of a growing body of mathematical results characterizing this class of networks that promises to be as rich as that for chemical reaction network models based on ordinary differential equations.


## 2. Background

The property of canalization for Boolean functions was introduced by Kauffman in [7], inspired by the concept of canalization from developmental biology [9]. A Boolean function is canalizing if there is a variable and a value of the variable such that if the variable takes the value, then the value of the function does not depend on other variables. It was shown that models defined by such functions often exhibit less chaotic and more stable behavior [10, 11]. Nested canalizing functions, obtained by applying the concept of canalization recursively, were introduced in [2]. They form a special subset of canalizing functions and have stable dynamics [11]. We note that there are other important properties shared by Boolean networks arising in modeling (for example, sparsity [7]). In this paper we focus only on canalization and its impact on the dynamics, and one of the natural future directions would be to consider several such properties simultaneously.To cover more models arising in applications, the notion of nested canalizing function was relaxed by Layne et al. [12] by assigning to every Boolean function its canalizing depth. Noncanalizing functions have canalizing depth zero, and nested canalizing functions have the maximal possible canalizing depth equal to the number of variables. Canalizing depth of a Boolean network is defined as the minimum of the canalizing depths of the functions defining the network. In [12], activities and sensitivities of functions with different canalizing depths and stability and criticality of Boolean networks composed from such functions were investigated. It has been observed that Boolean networks of higher canalizing depth tend to be more stable and less sensitive. However, increasing the canalizing depth to the maximum does not improve the stability significantly compared to moderate positive canalizing depth. These observations give a strong indication of the biological utility of canalizing function, even with small canalizing depth.Attractors in Boolean network models can be interpreted as distinct cell types [13, p. 202] and their lengths can be viewed as the variety of different gene expression patterns corresponding to the cell type. Thus, understanding the attractor structure of a random Boolean network defined by functions of a fixed canalizing depth is important for assessing biological relevance of such models. Analytic study of the attractor structure of nested canalizing Boolean networks has been carried out in [11]. For discussion about attractors of length one (i.e., steady state), we refer to [14].


## 3. Our Results

The main question of this paper is as follows:What do the dynamics of a random canalizing Boolean network look like? We approach this question using both computer simulations and analytical methods.In ourcomputational experiments, we generate approximately 30 million random Boolean networks of all possible canalizing depths with the number of variables ranging from 4 to 20. For each of these networks, we determine sizes of all the attractors and basins of attraction and analyze the obtained data. We discover the following:(1)
For a fixed number of variables, the sample mean of the number of attractors and average size of an attractor decrease when the canalizing depth increases(2)
The decrease of the average size of an attractor is much greater than the decrease of the number of attractors as the canalizing depth increases(3)
Both decreases from (8) are substantial when the canalizing depth changes from zero to small canalizing depths, but a further increase of the canalizing depth does not lead to a significant decrease for either the sample means or for the empirical distributions(4)
The relative decrease of the sample mean of the number of attractors and the average attractor size when the canalizing depth changes from zero to one becomes sharper when the number of variables increasesObservations (8) and (A.4) are consistent with the results obtained in [12] for sensitivity and stability. This provides new evidence that Boolean networks of small positive canalizing depth are almost as well-suited for modeling as those with nested canalizing functions, from the point of view of stability. Since there are many more canalizing functions of small positive canalizing depth than nested canalizing functions [15, Section 5], they provide a richer modeling toolbox.Motivated by observation (A.4), we conduct a mathematical study of the attractor structure of a random Boolean network of canalizing depth one (that is, the minimal positive depth). Our main theoretical result, Theorem 2, gives, for every positive integer ℓ, a formula for the limit of the expected number of attractors of length ℓ in a random Boolean network of depth one. The same formulas are valid for a random Boolean network defined by canalizing functions (see Remark 5). In particular, our formulas show that a large random network of depth one, on average, has more attractors of small sizes that an average Boolean network (Remark 6).Formulas similar to the ones in our proofs (e.g., in Lemma A.4) have already appeared in the study of the average number of attractors of a given length in sparse Boolean networks, e.g., [16, equation (2)] and [17, equation (6)]. The results of [16, 17] are based on describing the asymptotic behavior of these formulas in terms of N, the number of nodes in the network, and the asymptotics is of the form ONα. In our case, the average number of attractors of a given length simply approaches a constant as N⟶∞ (that is, O1), but our methods allow us to find the exact value of this constant.The source code we used for generating and analyzing data is available athttps://github.com/MathTauAthogen/Canalizing-Depth-Dynamics. The raw data are available at https://github.com/MathTauAthogen/Canalizing-Depth-Dynamics/tree/master/data.Structure of the Paper. The rest of the paper is organized as follows. Section 4 contains necessary definitions about canalizing functions and Boolean networks. Outlines of the algorithms used in our computational experiments are in Section 5. The main observations are summarized in Section 6. Our main theoretical result about attractors in a random Boolean network of canalizing depth one (Theorem 2) is presented in Section 7. Section 8 contains conclusions. The proofs are located in the Appendix.


## 4. Preliminaries

Definition 1.
ABoolean network is a tuple f=f1,f2,…,fn of Boolean functions in n variables. For a state at=at,1,at,2,…,at,n∈0,1n at time t, we define the state at+1≔fat=at+1,1,…,at+1,n∈0,1n at time t+1 by(1)at+1,1=f1at,1,…,at,n,⋮at+1,n=fnat,1,…,at,n.Definition 2.
(attractors and basins). Letf=f1,…,fn be a Boolean network.(i)
A sequencea1,…,aℓ∈0,1n of distinct states is called an attractor of f if fai=ai+1 for every 1≤i<ℓ and faℓ=a1.(ii)
An attractora1,…,aℓ∈0,1n is called a steady state if ℓ=1.(iii)
LetA=a1,…,aℓ∈0,1nℓ be an attractor of f. The basin of A is the set(2)b∈0,1n∃N:ff…fb…︸Ntimes∈A.Definition 3.
A nonconstant functionfx1,…,xn is canalizing with respect to a variable xi if there exists a canalizing value a∈0,1 such that(3)fx1,…,xi−1,a,xi+1,…,xn≡const.Example 1.
Considerfx1,x2=x1⋅x2 (the product is understood modulo 2, that is, logical AND). It is canalizing with respect to x1 with canalizing value 0 because f0,x2=0 regardless of the value of x2. Analogously, it is canalizing with respect to x2 with canalizing value 0.
Considergx1,x2=x1+x2 (summation is understood modulo 2, that is, logical XOR). It is not canalizing with respect to x1 because(4)g0,x2=x2≠const,g1,x2=x¯2≠const.
The same argument works forx2 as well.Definition 4.
f
x
1
,
…
,
x
n has canalizing depth [15, Definition 2.3] k if it can be expressed as(5)f=b1,xi1=a1,b2,xi1≠a1,xi2=a2,⋮bk,xi1≠a1,xi2≠a2,…,xik−1≠ak−1,xik=ak,g≢bk,xi1≠a1,…,xik≠ak,where(i)
i1,…,ik are distinct integers from 1 to n(ii)
a1,…,ak,b1,…,bk∈0,1(iii)
g is a noncanalizing function in the variables x1,…,xn/xi1,…,xikExample 2.
For example, iffx1,x2,x3=x1+x2x3,(6)fx1,x2,x3=0,x3=0,x1+x2,x3≠0,and x1+x2 is noncanalizing. Therefore, f has canalizing depth 1.Remark 1.
Sinceg in Definition 4 is noncanalizing, every function has a single well-defined canalizing depth. In particular, a function of depth two is not considered to have depth one.Definition 5.
We say that a canalizing Boolean functionfx1,…,xn is nested if f has canalizing depth n, that is, g=0 or g=1 (see Definition 4). For example, fx1,x2,x3=x1x2x3 is nested canalizing because(7)f=0,x3=0,0,x3≠0,x2=0,0,x2,x3≠0,x1=0,1,x1,x2,x3≠0,so the canalizing depth of f is 3, which is equal to n=3.Definition 6.
We say that a Boolean networkf=f1,…,fn has canalizing depth k if f1,…,fn are Boolean functions of canalizing depth k.


## 5. Simulations: Outline of the Algorithms

In our computational experiment, we generated random Boolean networks of various canalizing depths. For each network, we store a list of pairsai,bi, where ai is the size of the ith attractor of the network and bi is the size of its basin. The generated data are available at https://github.com/MathTauAthogen/Canalizing-Depth-Dynamics/tree/master/data. To generate the data, we used two algorithms: one for generating a random Boolean network of a given canalizing depth and one for finding the sizes of attractors and their basins (Algorithm 1).Algorithm 1: Finding the sizes of the attractors and their basins.
In: A Boolean network f=f1,…,fn in n variables
Out: A list of pairs ai,bi, where ai is the size of the ith attractor of f and bi is the size of its basin
(1)
(Network⟶ Graph) Build a directed graph G with 2n vertices corresponding to possible states and a directed edge from a to fa for every a∈0,1n.
(2)
(Attractors) Perform a depth-first search [18, § 22.3] traversal on G viewed as an undirected graph to detect the unique cycle in each connected component, these cycles are the attractors.
(3)
(Basins) For each cycle from Step 2, perform a depth-first search traversal onG with all the edges reversed. The dfs trees will be the basins.
(4)
Return the sizes of the attractors and basins found on Steps 2 and 3.


### 5.1. Generating Random Boolean Functions of a Given Canalizing Depth

[12, Section 5] contains a sketch of an algorithm for generating random Boolean functions that have canalizing depth at least k for a given k. Here, we generate functions of canalizing depth equal to k and take a different approach than [12]. In order to ensure that the probability distribution of possible outputs is uniform, we use the following structure theorem due to He and Macaulay [15].Theorem 1.
(see [15], Theorem 4.5). Every Boolean function fx1,…,xn≢0 can be uniquely written as(8)fx1,…,xn=M1M2⋯Mr−1MrpC+1+1⋯+1+b,where Mi=∏j=1kixij+aij for every 1≤i≤r, pC≢0 is a noncanalizing function, and k=∑i=1rki is the canalizing depth. Each xi appears in exactly one of M1,…,Mr,pC, and the only restrictions on equation (8) are the following “exceptional cases”:(E1) IfpC≡1 and r≠1, then kr≥2(E2) IfpC≡1 and r=1 and k1=1, then b=0Example 3.
Considerfx1,x2,x3,x4=x1x2+1x3x4+x3+x4 can be represented as(9)f=x1+0x2+1x3+1x4+11+1+0,so M1=x1+0x2+1, M2=x3+1x4+1, b=0, k=4, and pC=1. This can be verified by expanding the brackets in the original and new representations of f.
Considergx1,x2,x3,x4,x5=1+x5x1+x2x3+1x4. It can be represented as(10)g=x5+0x3+1x4+0x1+x2+1+1,so M1=x5+0, M2=x4+0x3+1, b=1, k=3, and pC=x1+x2.
Our algorithm is summarized in Algorithms2 and 3 below. Correctness of Algorithm 2 follows from Theorem 1, and correctness of Algorithm 3 can be proved directly by induction on k.Remark 2.
The complexity of Algorithm2 is On2n (see Proposition B.2). Given that the size of the output is O2n, and this is nearly optimal.
We measured the runtimes of our implementation of Algorithm2 on a laptop with a Core i5 processor (1.60 GHz) and 8 Gb RAM. Generating a single function with 20 variables (the largest number we used in our simulations) takes 4.9−5.5 seconds (faster for smaller canalizing depth). On a laptop, our implementation can go up to 24 variables (∼2 minutes to generate a function), and then hits memory limits. One can go further by using a lower level language and more careful packing. However, already a Boolean function in 40 variables would require at least 128 Gb of memory.Algorithm 2: Generating a random Boolean function of a given canalizing depth.
In: Nonnegative integers k and n with k≤n
Out: A Boolean function f in n variables of canalizing depth k such that, for fixed k and n, all possible outputs have the same probability
(1)
In the notation of Theorem1, generate the following:
(a)
random bitsb,a1,…,an∈0,1;
(b)
random subsetX⊂x1,…,xn with X=k;
(c)
random ordered partitionX=X1⊔⋯⊔Xr of X (using Algorithm 2);
(d)
random noncanalizing functionpC≢0 in variables x1,…,xn/X (see Remark 3).
(2)
Form a functionfx1,…,xn using the data generated in Step 1 as in Theorem 1, where Mi involves exactly the variables from Xi for every 1≤i≤r.
(3)
Iff does not satisfy any of the conditions (E1) or (E2), discard it and run the algorithm again. Otherwise, return f.Remark 3.
We generate a random noncanalizing function as follows. We generate a random Boolean function and test for canalization until we generate a noncanalizing one. Then, we return it. Since canalizing functions are rare [15, Section 5], this algorithm is fast enough for our purposes (see Lemma B.1).Algorithm 3: Generating a random ordered partition of a given finite set.
In: A finite set X with X=k
Out: An ordered partition X=X1⊔⋯⊔Xr into nonempty subsets X1,…,Xr such that, for a fixed X, all possible outputs have the same probability
(1)
Computep0,…,pk, where pi is the number of ordered partitions of a set of size i, using the recurrence pj=∑i=0j−1jipj−i, p0=1 (see [19, equation (9)]).
(2)
Generate an integerN uniformly at random from 1,pk.
(3)
Find the minimum integerj between 1 and k such that ∑i=0j−1kipk−i≥N.
(4)
Randomly select a subsetX1⊂X of size j.
(5)
Generate an ordered partitionX2⊔⋯⊔Xr of X/X1 recursively.
(6)
ReturnX1⊔⋯⊔Xr.


## 5.1. Generating Random Boolean Functions of a Given Canalizing Depth

[12, Section 5] contains a sketch of an algorithm for generating random Boolean functions that have canalizing depth at least k for a given k. Here, we generate functions of canalizing depth equal to k and take a different approach than [12]. In order to ensure that the probability distribution of possible outputs is uniform, we use the following structure theorem due to He and Macaulay [15].Theorem 1.
(see [15], Theorem 4.5). Every Boolean function fx1,…,xn≢0 can be uniquely written as(8)fx1,…,xn=M1M2⋯Mr−1MrpC+1+1⋯+1+b,where Mi=∏j=1kixij+aij for every 1≤i≤r, pC≢0 is a noncanalizing function, and k=∑i=1rki is the canalizing depth. Each xi appears in exactly one of M1,…,Mr,pC, and the only restrictions on equation (8) are the following “exceptional cases”:(E1) IfpC≡1 and r≠1, then kr≥2(E2) IfpC≡1 and r=1 and k1=1, then b=0Example 3.
Considerfx1,x2,x3,x4=x1x2+1x3x4+x3+x4 can be represented as(9)f=x1+0x2+1x3+1x4+11+1+0,so M1=x1+0x2+1, M2=x3+1x4+1, b=0, k=4, and pC=1. This can be verified by expanding the brackets in the original and new representations of f.
Considergx1,x2,x3,x4,x5=1+x5x1+x2x3+1x4. It can be represented as(10)g=x5+0x3+1x4+0x1+x2+1+1,so M1=x5+0, M2=x4+0x3+1, b=1, k=3, and pC=x1+x2.
Our algorithm is summarized in Algorithms2 and 3 below. Correctness of Algorithm 2 follows from Theorem 1, and correctness of Algorithm 3 can be proved directly by induction on k.Remark 2.
The complexity of Algorithm2 is On2n (see Proposition B.2). Given that the size of the output is O2n, and this is nearly optimal.
We measured the runtimes of our implementation of Algorithm2 on a laptop with a Core i5 processor (1.60 GHz) and 8 Gb RAM. Generating a single function with 20 variables (the largest number we used in our simulations) takes 4.9−5.5 seconds (faster for smaller canalizing depth). On a laptop, our implementation can go up to 24 variables (∼2 minutes to generate a function), and then hits memory limits. One can go further by using a lower level language and more careful packing. However, already a Boolean function in 40 variables would require at least 128 Gb of memory.Algorithm 2: Generating a random Boolean function of a given canalizing depth.
In: Nonnegative integers k and n with k≤n
Out: A Boolean function f in n variables of canalizing depth k such that, for fixed k and n, all possible outputs have the same probability
(1)
In the notation of Theorem1, generate the following:
(a)
random bitsb,a1,…,an∈0,1;
(b)
random subsetX⊂x1,…,xn with X=k;
(c)
random ordered partitionX=X1⊔⋯⊔Xr of X (using Algorithm 2);
(d)
random noncanalizing functionpC≢0 in variables x1,…,xn/X (see Remark 3).
(2)
Form a functionfx1,…,xn using the data generated in Step 1 as in Theorem 1, where Mi involves exactly the variables from Xi for every 1≤i≤r.
(3)
Iff does not satisfy any of the conditions (E1) or (E2), discard it and run the algorithm again. Otherwise, return f.Remark 3.
We generate a random noncanalizing function as follows. We generate a random Boolean function and test for canalization until we generate a noncanalizing one. Then, we return it. Since canalizing functions are rare [15, Section 5], this algorithm is fast enough for our purposes (see Lemma B.1).Algorithm 3: Generating a random ordered partition of a given finite set.
In: A finite set X with X=k
Out: An ordered partition X=X1⊔⋯⊔Xr into nonempty subsets X1,…,Xr such that, for a fixed X, all possible outputs have the same probability
(1)
Computep0,…,pk, where pi is the number of ordered partitions of a set of size i, using the recurrence pj=∑i=0j−1jipj−i, p0=1 (see [19, equation (9)]).
(2)
Generate an integerN uniformly at random from 1,pk.
(3)
Find the minimum integerj between 1 and k such that ∑i=0j−1kipk−i≥N.
(4)
Randomly select a subsetX1⊂X of size j.
(5)
Generate an ordered partitionX2⊔⋯⊔Xr of X/X1 recursively.
(6)
ReturnX1⊔⋯⊔Xr.


## 6. Simulations: Results

Notation 1.
For a Boolean networkf=f1,…,fn, let Nf and Sf denote the number of the attractors of f and the sum of the sizes of the attractors of f, respectively. We define the average size of an attractor as ASf≔Sf/Nf.


### 6.1. Sample Means ofNf and ASf

For everyn=4,…,20 and every 0≤k≤n, we generate random Boolean networks in n variables of canalizing depth k and compute the mean of Nf and ASf. Figure 1 shows how these means depend on k for n=15 (based on 50,000 samples for each k). The shape of the plots is similar for other values of n we did computation for (that is, n=4,…,20). Note that although both means are decreasing, the decrease of the mean of ASf is more substantial.Figure 1
Dependencies of the sample means ofNf and ASf on the canalizing depth. (a) The number of attractors Nf. (b) Average size of an attractor ASf.
(a)
(b)


### 6.2. Distributions ofNf and ASf

Figure2 shows the empirical distributions of Nf and ASf for n=12 and k=0,1,3,12 based on 300,000 samples for each k. From the plot, we can make the following observations:(i)
The distributions become more concentrated and the peak shifts towards zero whenk increases(ii)
The distributions for nonzero canalizing depths (especially for larger depths) are much closer to each other that to the distribution for zero canalizing depth. This agrees with the plots on Figure1.Figure 2
Empirical distributions ofNf and ASf for n=12 and k=0,1,3,12. (a) Distribution of the number of attractors Nf. (b) Distribution of the average size of an attractor ASf.
(a)
(b)


### 6.3. Relative Decreases

From Figure1, we can observe that, for both Nf and ASf, the sample mean decreases rapidly for small canalizing depths. In order to understand how this decrease behaves for large n, we introduce(11)Nkn≔the sample mean ofNffornvariables and canalizing depthkthe samplemean ofNffornvariables and canalizing depth0.AS
k
n is defined analogously. Figure 3 plots N1n, N2n, N3n, and Nnn and AS1n, AS2n, AS3n, and ASnn as functions of n. From the plots we see that(i)
The relative initial decrease from canalizing depth 0 to canalizing depth 1 becomes even more substantial whenn increases(ii)
The relative decrease from canalizing depth 0 to canalizing depth 3 is already very close to the relative decrease from depth zero to the maximal depth (i.e., nested canalizing functions)Figure 3
Dependence of the relative decreases of the sample means ofNf and ASf on the number of variables n. (a) Relative decrease of the number of attractors. (b) Relative decrease of the average size of an attractor.
(a)
(b)


## 6.1. Sample Means ofNf and ASf

For everyn=4,…,20 and every 0≤k≤n, we generate random Boolean networks in n variables of canalizing depth k and compute the mean of Nf and ASf. Figure 1 shows how these means depend on k for n=15 (based on 50,000 samples for each k). The shape of the plots is similar for other values of n we did computation for (that is, n=4,…,20). Note that although both means are decreasing, the decrease of the mean of ASf is more substantial.Figure 1
Dependencies of the sample means ofNf and ASf on the canalizing depth. (a) The number of attractors Nf. (b) Average size of an attractor ASf.
(a)
(b)


## 6.2. Distributions ofNf and ASf

Figure2 shows the empirical distributions of Nf and ASf for n=12 and k=0,1,3,12 based on 300,000 samples for each k. From the plot, we can make the following observations:(i)
The distributions become more concentrated and the peak shifts towards zero whenk increases(ii)
The distributions for nonzero canalizing depths (especially for larger depths) are much closer to each other that to the distribution for zero canalizing depth. This agrees with the plots on Figure1.Figure 2
Empirical distributions ofNf and ASf for n=12 and k=0,1,3,12. (a) Distribution of the number of attractors Nf. (b) Distribution of the average size of an attractor ASf.
(a)
(b)


## 6.3. Relative Decreases

From Figure1, we can observe that, for both Nf and ASf, the sample mean decreases rapidly for small canalizing depths. In order to understand how this decrease behaves for large n, we introduce(11)Nkn≔the sample mean ofNffornvariables and canalizing depthkthe samplemean ofNffornvariables and canalizing depth0.AS
k
n is defined analogously. Figure 3 plots N1n, N2n, N3n, and Nnn and AS1n, AS2n, AS3n, and ASnn as functions of n. From the plots we see that(i)
The relative initial decrease from canalizing depth 0 to canalizing depth 1 becomes even more substantial whenn increases(ii)
The relative decrease from canalizing depth 0 to canalizing depth 3 is already very close to the relative decrease from depth zero to the maximal depth (i.e., nested canalizing functions)Figure 3
Dependence of the relative decreases of the sample means ofNf and ASf on the number of variables n. (a) Relative decrease of the number of attractors. (b) Relative decrease of the average size of an attractor.
(a)
(b)


## 7. Theory: The Main Result

We will introduce notation needed to state the main theorem. Let us fix a positive integerℓ. For a binary string α∈S≔0,1ℓ, we define(i)
α denotes the number of ones(ii)
α¯ denotes component-wise negation(iii)
sα denotes a cyclic shift to the rightFor binary stringsα,β∈0,1ℓ, we define(12)fα,β≔12β,ifα∨β=β,0,otherwise,gα,β≔14fα,β+fα¯,β+fα,β¯+fα¯,β¯.Then, we define a2ℓ×2ℓ matrix Gℓ by(13)Gℓa,b=ga,sb,where we interpret numbers 1≤aandb≤2ℓ as binary sequences of length ℓ.Theorem 2.
LetAℓ be the limit of the expected number of attractors of length ℓ in a random Boolean network of canalizing depth one (see Definition 6) when the number of variables n goes to infinity. Then,(14)Aℓ=1ℓPGℓ′1,where PGℓ is the characteristic polynomial of matrix Gℓ introduced above. In particular, we have(15)A1=1,A2=23=0.666…,A3=64189=0.3386…,A4=0.2856…,A5=0.2004…,A6=0.1721….Remark 4.
The plots below show that the result of Theorem2 agrees with our simulations (Figure 4).Remark 5.
As explained in Remark A.1, Theorem2 stills holds if we replace a random Boolean network of canalizing depth one with a random Boolean network defined by canalizing functions.Example 4.
Letℓ=2. Then, for example, we have f0,2=f0,1=1/2 and g0,1=g3,1=1/4. In total, we have(16)G2=3/81/41/43/81/81/41/41/81/81/41/41/83/81/41/43/8,PG2t=t4−54t3+12t2.Remark 6.
Theorem2 and Corollary A.1 imply that Aℓ>1/ℓ for every ℓ>1. On the other hand, a direct computation shows that the expected number of attractors of length ℓ in a random Boolean network (without any canalization requirements) is 1/ℓ. This is consistent with our observations from Section 6.1.Remark 7.
A sage script for computing numbersAℓ is available at https://github.com/MathTauAthogen/Canalizing-Depth-Dynamics/blob/master/core/theory.sage.Figure 4
The average number of attractors of fixed length (blue plot) compared to the limiting value from Theorem 2 (orange plot). (a) Length 1. (b) Length 2. (c) Length 3. (d) Length 4.
(a)
(b)
(c)
(d)


## 8. Conclusions

We conducted computational experiments to investigate the attractor structure of Boolean networks defined by functions of varying canalizing depth. We observed that networks with higher canalizing depth tend to have fewer attractors and the sizes of the attractors decrease dramatically when the canalizing depth increases moderately. As a consequence, the basins tend to grow when the canalizing depth increases. These properties are desirable in many biological applications of Boolean networks, so our results give new indications of the biological utility of Boolean networks defined by functions of positive canalizing depth.We proved a theoretical result, Theorem2, which complements the above observation as follows. The theorem implies that a large random Boolean network of canalizing depth one has on average more attractors of small size than a random Boolean network of the same size although it has less attractors in total. This also provides an explanation to the fact that the total size of attractors decreases faster than the number of attractors as the canalizing depth grows.Furthermore, we observed that all the statistics we computed are almost the same in the case of the maximal possible canalizing depth (so-called nested canalizing Boolean networks) and in the case of moderate canalizing depth. This agrees with the results of Layne et al. [12]. This observation elucidates an interesting and powerful feature of canalization: even a very moderate canalizing influence in a Boolean network has a strong constraining influence on network dynamics. It would be of interest to explore the prevalence of these features in published Boolean network models.Finally, we provided evidence that the observed phenomena will occur for Boolean networks with larger numbers of state variables.



---

*Source: 3687961-2020-01-20.xml*