The Dynamics of Canalizing Boolean Networks

TL;DR

This paper explores how canalizing properties influence the dynamics of Boolean networks, revealing that higher canalizing depth generally leads to fewer and smaller attractors, with implications for biological modeling.

Contribution

It provides analytical and simulation-based insights into the impact of canalizing depth on Boolean network dynamics, including explicit formulas for attractor counts.

Findings

Higher canalizing depth results in fewer attractors.

Attractors are smaller and basins are larger with increased canalizing depth.

High canalizing depth has a modest impact on attractor structure.

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 $\ell$, we give an explicit formula for the limit of the expected number of attractors of length $\ell$ in an $n$-state random Boolean network as $n$ goes to infinity.

Equations

$$a_{t+1,1}=f_1(a_{t,1},\dots ,a_{t,n}),$$

$$⋮$$

$$a_{t+1,n}=f_n(a_{t,1},\dots ,a_{t,n}).$$

$$\{\mathbf{b}\in \{0,1{\}}^n|\exists N:\underset{N\text{ times}}{\underbrace{\mathbf{f}(\mathbf{f}(\dots (\mathbf{f}(}}\mathbf{b})\dots ))\in A\}.$$

$$f(x_1,\dots ,x_{i-1},a,x_{i+1},\dots ,x_n)\equiv \text{const}.$$

$$g(0,x_2)=x_2\ne \mathrm{const}\text{ and }g(1,x_2)={\bar{x}}_2\ne \mathrm{const}.$$

$$f=\{\begin{array}{ll}{b}_1& x_{i_1}=a_1\\ b_2& x_{i_1}\ne a_1,x_{i_2}=a_2\\ ⋮& \\ b_k& x_{i_1}\ne a_1,x_{i_2}\ne a_2\dots x_{i_{k-1}}\ne a_{k-1},x_{i_k}=a_k\\ g\not\equiv b_k& x_{i_1}\ne a_1,\dots ,x_{i_k}\ne a_k,\end{array}$$

$$f(x_1,x_2,x_3)=\{\begin{array}{ll}0,& x_3=0\\ x_1+x_2,& x_3\ne 0\end{array}$$

$$f=\{\begin{array}{ll}0& x_3=0\\ 0& x_3\ne 0,x_2=0\\ 0& x_2,x_3\ne 0,x_1=0\\ 1& x_1,x_2,x_3\ne 0\end{array}$$

$$f(x_1,\dots ,x_n)=M_1(M_2(\cdots (M_{r-1}(M_r{p}_C+1)+1)\cdots )+1)+b,$$

$$f=((x_1+0)(x_2+1))(((x_3+1)(x_4+1))(1)+1)+0,$$

$$g=(x_5+0)(((x_3+1)(x_4+0))(x_1+x_2)+1)+1,$$

$$N_k(n):=\frac{\text{the sample mean of N(f) for n variables and canalizing depth k}}{\text{the sample mean of N(f) for n variables and canalizing depth 0}}.$$

$$f(\alpha ,\beta ):=\{\begin{array}{l}\frac{1}{2^{|\beta |}},\text{ if }\alpha \vee \beta =\beta ,\\ 0,\text{otherwise}\end{array}\text{ and }g(\alpha ,\beta ):=\frac{1}{4}(f(\alpha ,\beta )+f(\bar{\alpha },\beta )+f(\alpha ,\bar{\beta })+f(\bar{\alpha },\bar{\beta })).$$

$$(G_{\ell }{)}_{a,b}=g(a,s(b)),$$

$$A_{\ell }=\frac{1}{\ell P_{G_{\ell }}^{\prime }(1)},$$

$$A_1$$

$$A_2$$

$$A_3$$

$$G_2=\left(\begin{array}{cccc}3/8& 1/4& 1/4& 3/8\\ 1/8& 1/4& 1/4& 1/8\\ 1/8& 1/4& 1/4& 1/8\\ 3/8& 1/4& 1/4& 3/8\end{array}\right)\text{ and }{P}_{G_2}(t)=t^4-\frac{5}{4}{t}^3+\frac{1}{2}{t}^2.$$

$$S_{i,j}:=\{({\alpha }_1,\dots ,{\alpha }_{\ell })\in S\mid {\alpha }_i={\alpha }_j\}.$$

$$\sum _{\alpha \in S}g(\alpha ,s(\beta ))=1.$$

$$\sum _{\alpha \in S}h(\alpha ,\beta )=\sum _{k=0}^{|\beta |}\left(\genfrac{}{}{0ex}{}{|\beta |}{k}\right)\frac{1}{2^{|\beta |}}=1⟹\sum _{\beta \in S}g(\alpha ,\beta )=1.$$

$$\sum _{\alpha \in S_{i,j}}g(\alpha ,s(\beta ))⩽\frac{2^{\ell +2}-1}{2^{\ell +2}}.$$

$$\sum _{\alpha \in S_{i,j}}g(\alpha ,s(\beta ))⩽\left(\sum _{\alpha \in S}g(\alpha ,s(\beta ))\right)-g(\gamma ,s(\beta ))⩽1-\frac{1}{2^{\ell +2}}.$$

$$\mathbf{a}!:=a_1!\cdot \dots \cdot a_n!,{\mathbf{a}}^{\mathbf{b}}:=a_1^{b_1}\cdot \dots \cdot a_n^{b_n},|\mathbf{a}|:=a_1+\dots +a_n.$$

$$C(A{)}_n:=\sum _{\begin{array}{c}\mathbf{m}\in {\mathbb{Z}}_{⩾0}^s\\ |\mathbf{m}|=n\end{array}}\frac{(A\mathbf{m}{)}^{\mathbf{m}}}{n^n}.$$

$$y(z)=1-\epsilon +\frac{1}{3}{\epsilon }^2-\dots ,\text{ where }\epsilon =\sqrt{2-2ez}.$$

$$\frac{1}{y(z)}=1+\epsilon +\frac{2}{3}{\epsilon }^2-\dots ,\text{ where }\epsilon =\sqrt{2-2ez}.$$

$$\sum\limits_{\mathbf{m}\in\mathbb{Z}_{\geqslant 0}^{s}}\frac{(A\mathbf{m})^{\mathbf{m}}}{\mathbf{m}!}x^{|\mathbf{m}|}\exp\bigl{(}-x\sum\limits_{i,j}m_{j}a_{i,j}\bigr{)}=\frac{1}{\det|E-xA|}.$$

$$\sum _{\mathbf{m}\in {\mathbb{Z}}_{⩾0}^s}\frac{(A\mathbf{m}{)}^{\mathbf{m}}}{\mathbf{m}!}{x}^{|\mathbf{m}|}{e}^{-x|\mathbf{m}|}=\frac{1}{\det |E-xA|}.$$

$$\sum _{\mathbf{m}\in {\mathbb{Z}}_{⩾0}^s}\frac{(A\mathbf{m}{)}^{\mathbf{m}}}{\mathbf{m}!}{z}^{|\mathbf{m}|}=\frac{1}{\det |E-y(z)A|}.$$

Full text

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 $\ell$, we give an explicit formula for the limit of the expected number of attractors of length $\ell$ in an $n$-state random Boolean network as $n$ goes to infinity.

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, and the questions to be answered, 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 properly 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., [4, 15, 18, 22, 19]), going back to the work of S. Kauffman in the 1960s [11, 12, 13]. 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 S. Kauffman who proposed a class of functions [12] with properties inspired by the developmental biology concept of canalization, going back to C. Waddington in the 1940s [24]. 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: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 7.1, 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 S. Kauffman in [12], inspired by the concept of canalization from developmental biology [24]. 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 [16, 10]. Nested canalizing functions, obtained by applying the concept of canalization recursively, were introduced in [15]. They form a special subset of canalizing functions and have stable dynamics [16]. We note that there are other important properties shared by Boolean networks arising in modeling (for example, sparsity [12]). 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, Dimitrova, and Macaulay [17] by assigning to every Boolean function its canalizing depth. Non-canalizing 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 [17], 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 [14, 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 done in [16]. For discussion about attractors of length one (i.e., steady state), we refer to [23].

3 Our results.

The main question of this paper is:What do the dynamics of a random canalizing Boolean network look like? We approach this question using both computer simulations and analytical methods.

In our computational 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. (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. (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. (3)

Both decreases from (1) 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. (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 increases.

The observations (1) and (3) are consistent with the results obtained in [17] 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 [9, Section 5], they provide a richer modeling toolbox.

Motivated by the observation (3), 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 7.1, gives, for every positive integer $\ell$, a formula for the limit of the expected number of attractors of length $\ell$ in a random Boolean network of depth one. The same formulas are valid for a random Boolean network defined by canalizing functions (see Remark 7.2). 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 7.3).

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. [20, Eq. (2)] and [5, Eq. (6)]. The results of [20] and [5] are based on describing the asymptotic behavior of these formulas in terms of $N$, the number of nodes in the network, and the asymptotics are of the form $\mathrm{O}(N^{\alpha})$. In our case, the average number of attractors of a given length simply approaches a constant as $N\to\infty$ (that is, is $\mathrm{O}(1)$), 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 at https://github.com/MathTauAthogen/Canalizing-Depth-Dynamics. The raw data is 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 7.1) is presented in Section 7. Section 8 contains conclusions. The proofs are located in the Appendix.

4 Preliminaries

Definition 4.1**.**

A Boolean network is a tuple $\mathbf{f}=(f_{1},f_{2},\ldots,f_{n})$ of Boolean functions in $n$ variables. For a state $\bm{a}_{t}=(a_{t,1},a_{t,2},\ldots,a_{t,n})\in\{0,1\}^{n}$ at time $t$, we define the state $\bm{a}_{t+1}:=\mathbf{f}(\bm{a}_{t})=(a_{t+1,1},\ldots,a_{t+1,n})\in\{0,1\}^{n}$ at time $t+1$ by

[TABLE]

Definition 4.2** (Attractors and basins).**

Let $\mathbf{f}=(f_{1},\ldots,f_{n})$ be a Boolean network.

• •

A sequence $\bm{a}_{1},\ldots,\bm{a}_{\ell}\in\{0,1\}^{n}$ of distinct states is called an attractor of $\mathbf{f}$ if $\mathbf{f}(\bm{a}_{i})=\bm{a}_{i+1}$ for every $1\leqslant i<\ell$ and $\mathbf{f}(\bm{a}_{\ell})=\bm{a}_{1}$.

• •

An attractor $\bm{a}_{1},\ldots,\bm{a}_{\ell}\in\{0,1\}^{n}$ is called a steady state if $\ell=1$.

• •

Let $A=(\bm{a}_{1},\ldots,\bm{a}_{\ell})\in\left(\{0,1\}^{n}\right)^{\ell}$ be an attractor of $\mathbf{f}$. The basin of $A$ is the set

[TABLE]

Definition 4.3**.**

A nonconstant function $f(x_{1},\ldots x_{n})$ is canalizing with respect to a variable $x_{i}$ if there exists a canalizing value $a\in\{0,1\}$ such that

[TABLE]

Example 4.1**.**

Consider $f(x_{1},x_{2})=x_{1}\cdot x_{2}$ (the product is understood modulo 2, that is, logical AND). It is canalizing with respect to $x_{1}$ with canalizing value [math], because $f(0,x_{2})=0$ regardless of the value of $x_{2}$. Analogously, it is canalizing with respect to $x_{2}$ with canalizing value [math].

Consider $g(x_{1},x_{2})=x_{1}+x_{2}$ (summation is understood modulo 2, that is, logical XOR). It is not canalizing with respect to $x_{1}$, because

[TABLE]

The same argument works for $x_{2}$ as well.

Definition 4.4**.**

$f(x_{1},\ldots,x_{n})$ has canalizing depth [9, Definition 2.3] $k$ if it can be expressed as

[TABLE]

where

• •

$i_{1},\ldots,i_{k}$ are distinct integers from $1$ to $n$;

• •

$a_{1},\ldots,a_{k},b_{1},\ldots,b_{k}\in\{0,1\}$;

• •

$g$ is a noncanalizing function in the variables $\{x_{1},\ldots,x_{n}\}\setminus\{x_{i_{1}},\ldots,x_{i_{k}}\}$.

Example 4.2**.**

For example, if $f(x_{1},x_{2},x_{3})=(x_{1}+x_{2})x_{3}$,

[TABLE]

and $x_{1}+x_{2}$ is noncanalizing. Therefore $f$ has canalizing depth 1.

Remark 4.1**.**

Since $g$ in Definition 4.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 4.5**.**

We say that a canalizing Boolean function $f(x_{1},\ldots,x_{n})$ is nested if $f$ has canalizing depth $n$; that is, $g=0$ or $g=1$ (see Definition 4.4). For example, $f(x_{1},x_{2},x_{3})=x_{1}x_{2}x_{3}$ is nested canalizing because

[TABLE]

so the canalizing depth of $f$ is 3, which is equal to $n=3$.

Definition 4.6**.**

We say that a Boolean network $\mathbf{f}=(f_{1},\ldots,f_{n})$ has canalizing depth $k$ if $f_{1},\ldots,f_{n}$ 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 pairs $(a_{i},b_{i})$, where $a_{i}$ is the size of the $i$-th attractor of the network, and $b_{i}$ is the size of its basin. The generated data is 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.

5.1 Finding the sizes of the attractors and their basins

5.2 Generating random Boolean functions of a given canalizing depth

[17, 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 [17]. In order to ensure that the probability distribution of possible outputs is uniform, we use the following structure theorem due to He and Macaulay [9].

Theorem 5.1** ([9, Theorem 4.5]).**

Every Boolean function $f(x_{1},\ldots,x_{n})\not\equiv 0$ can be uniquely written as

[TABLE]

where $\displaystyle M_{i}=\prod_{j=1}^{k_{i}}(x_{i_{j}}+a_{i_{j}})$ for every $1\leqslant i\leqslant r$, $p_{C}\not\equiv 0$ is a noncanalizing function, and $k=\sum\limits_{i=1}^{r}k_{i}$ is the canalizing depth. Each $x_{i}$ appears in exactly one of $\{M_{1},\dots,M_{r},p_{C}\}$, and the only restrictions on Eq. (1) are the following “exceptional cases”:

1. (E1)

If $p_{C}\equiv 1$ and $r\neq 1$, then $k_{r}\geq 2$; 2. (E2)

If $p_{C}\equiv 1$ and $r=1$ and $k_{1}=1$, then $b=0$.

Example 5.1**.**

Consider $f(x_{1},x_{2},x_{3},x_{4})=x_{1}(x_{2}+1)(x_{3}x_{4}+x_{3}+x_{4})$ can be represented as

[TABLE]

so $M_{1}=(x_{1}+0)(x_{2}+1)$, $M_{2}=(x_{3}+1)(x_{4}+1)$, $b=0$, $k=4$, and $p_{C}=1$. This can be verified by expanding the brackets in the original and new representations of $f$.

Consider $g(x_{1},x_{2},x_{3},x_{4},x_{5})=1+x_{5}(x_{1}+x_{2})(x_{3}+1)x_{4}$. It can be represented as

[TABLE]

so $M_{1}=(x_{5}+0)$, $M_{2}=(x_{4}+0)(x_{3}+1)$, $b=1$, $k=3$, and $p_{C}=x_{1}+x_{2}$.

Our algorithm is summarized in Algorithms 2 and 3 below. Correctness of Algorithm 2 follows from Theorem 5.1, and correctness of Algorithm 3 can be proved directly by induction on $k$.

Remark 5.1**.**

The complexity of Algorithm 2 is $\mathrm{O}(n2^{n})$ (see Proposition B.2). Given that the size of the output is $\mathrm{O}(2^{n})$, this is nearly optimal.

We measured the runtimes of our implementation of Algorithm 2 on a laptop with a Core i5 processor (1.60GHz) and 8Gb 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 ($\sim 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 128Gb of memory.

Remark 5.2**.**

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 [9, Section 5], this algorithm is fast enough for our purposes (see Lemma B.1).

6 Simulations: results

Notation 6.1**.**

For a Boolean network $\mathbf{f}=(f_{1},\ldots,f_{n})$, let $N(\mathbf{f})$ and $S(\mathbf{f})$ denote the number of the attractors of $\mathbf{f}$ and the sum of the sizes of the attractors of $\mathbf{f}$, respectively. We define the average size of an attractor as $\operatorname{AS}(\mathbf{f}):=\frac{S(\mathbf{f})}{N(\mathbf{f})}$.

6.1 Sample means of $N(\mathbf{f})$ and $\operatorname{AS}(\mathbf{f})$

For every $n=4,\ldots,20$ and every $0\leqslant k\leqslant n$, we generate random Boolean networks in $n$ variables of canalizing depth $k$ and compute the mean of $N(\mathbf{f})$ and $\operatorname{AS}(\mathbf{f})$. Figure 1 shows how these means depend on $k$ for $n=15$ (based on $50000$ samples for each $k$). The shape of the plots is similar for other values of $n$ we did computation for (that is, $n=4,\ldots,20$). Note that although both means are decreasing, the decrease of the mean of $\operatorname{AS}(\mathbf{f})$ is more substantial.

6.2 Distributions of $N(\mathbf{f})$ and $\operatorname{AS}(\mathbf{f})$

Figure 2 shows the empirical distributions of $N(\mathbf{f})$ and $\operatorname{AS}(\mathbf{f})$ for $n=12$ and $k=0,1,3,12$ based on $300000$ samples for each $k$. From the plot, we can make the following observations.

• •

The distributions become more concentrated and the peak shifts towards zero when $k$ increases.

• •

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 Figure 1.

6.3 Relative decreases

From Figure 1, we can observe that, for both $N(\mathbf{f})$ and $\operatorname{AS}(\mathbf{f})$, the sample mean decreases rapidly for small canalizing depths. In order to understand how this decrease behaves for large $n$, we introduce

[TABLE]

$\operatorname{AS}_{k}(n)$ is defined analogously. Figure 3 plots $N_{1}(n)$, $N_{2}(n)$, $N_{3}(n)$, $N_{n}(n)$ and $\operatorname{AS}_{1}(n)$, $\operatorname{AS}_{2}(n)$, $\operatorname{AS}_{3}(n)$, $\operatorname{AS}_{n}(n)$ as functions of $n$. From the plots we see that

• •

the relative initial decease from canalizing depth [math] to canalizing depth $1$ becomes even more substantial when $n$ increases;

• •

the relative decrease from canalizing depth [math] to canalizing depth $3$ is already very close to the relative decrease from depth zero to the maximal depth (i.e., nested canalizing functions).

7 Theory: the main result

We will introduce notation needed to state the main theorem. Let us fix a positive integer $\ell$. For a binary string $\alpha\in S:=\{0,1\}^{\ell}$, we define:

• •

$|\alpha|$ denotes the number of ones;

• •

$\bar{\alpha}$ denotes component-wise negation;

• •

$s(\alpha)$ denotes a cyclic shift to the right.

For binary strings $\alpha,\beta\in\{0,1\}^{\ell}$ we define

[TABLE]

Then we define a $2^{\ell}\times 2^{\ell}$ matrix $G_{\ell}$ by

[TABLE]

where we interpret numbers $1\leqslant a,b\leqslant 2^{\ell}$ as binary sequences of length $\ell$.

Theorem 7.1**.**

Let $A_{\ell}$ be the limit of the expected number of attractors of length $\ell$ in a random Boolean network of canalizing depth one (see Definition 4.6) when the number of variables $n$ goes to infinity. Then

[TABLE]

where $P_{G_{\ell}}$ is the characteristic polynomial of matrix $G_{\ell}$ introduced above. In particular, we have

[TABLE]

Remark 7.1**.**

The plots below show that the result of Theorem 7.1 agrees with our simulations

Remark 7.2**.**

As explained in Remark A.1, Theorem 7.1 stills holds if we replace a random Boolean network of canalizing depth one with a random Boolean network defined by canalizing functions.

Example 7.1**.**

Let $\ell=2$. Then, for example, we have $f(0,2)=f(0,1)=\frac{1}{2}$ and $g(0,1)=g(3,1)=\frac{1}{4}$. In total, we have

[TABLE]

Remark 7.3**.**

Theorem 7.1 and Corollary A.1 imply that $A_{\ell}>\frac{1}{\ell}$ for every $\ell>1$. On the other hand, a direct computation shows that the expected number of attractors of length $\ell$ in a random Boolean network (without any canalization requirements) is $\frac{1}{\ell}$. This is consistent with our observations from Section 6.1.

Remark 7.4**.**

A Sage script for computing numbers $A_{\ell}$ is available at https://github.com/MathTauAthogen/Canalizing-Depth-Dynamics/blob/master/core/theory.sage.

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, Theorem 7.1, 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, Dimitrova, and Macauley [17].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 this 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.

Acknowledgments

The authors are grateful to Claus Kadelka, Christian Krattenthaler, and Doron Zeilberger for helpful discussions and to the referees for valuable suggestions and comments. GP was partially supported by NSF grants CCF-1564132, CCF-1563942, DMS-1760448, DMS-1853482, and DMS-1853650 by PSC-CUNY grants #69827-0047 and #60098-0048. RL was partially supported by Grants NIH 1U01EB024501-01 and NSF CBET-1750183. EP, GP, and WQ are grateful to the New York Math Circle where their collaboration started.

Data Availability Statement

Python/sage code and the results of simulations used to support the findings of this study have been deposited at https://github.com/MathTauAthogen/Canalizing-Depth-Dynamics.

Appendix A: Proofs

Notation A.1**.**

We fix a positive integer $\ell$.

• •

For every $1\leqslant i<j\leqslant\ell$, we define a subset $S_{i,j}\subset S=\{0,1\}^{\ell}$ by

[TABLE]

• •

For every $1\leqslant i<j\leqslant\ell$, let $G_{\ell;i,j}$ be the submatrix of $G_{\ell}$ with rows and columns having indices from $S_{i,j}$.

Lemma A.1**.**

For every $\ell$, we have

$G_{\ell}^{T}$* is stochastic (see [21, § 8.5]), and $G_{\ell}$ has exactly one eigenvalue being equal to $1$.* 2. 2.

for every $1\leqslant i<j\leqslant\ell$, there exists a $2^{\ell-1}\times 2^{\ell-1}$-matrix $C_{\ell;i,j}$ with nonnegative entries such that $\frac{2^{\ell+2}}{2^{\ell+2}-1}(G_{\ell;i,j}+C_{\ell;i,j})^{T}$ is stochastic, and has exactly one of the eigenvalues being equal to $1$.

Proof.

We will first show that $G_{\ell}^{T}$ is stochastic and irreducible (see [21, § 3.11]).

By definition, showing that $G_{\ell}^{T}$ is stochastic is equivalent to proving that, for every $\beta\in S:=\{0,1\}^{\ell}$,

[TABLE]

Since shift just permutes binary strings, this sum is equal to $\sum\limits_{\beta\in S}g(\alpha,\beta)$. For a fixed $\beta$ and $k\leqslant|\beta|$, the number of $\alpha\in S$ such $\alpha\vee\beta=\beta$ and $|\alpha|=k$ is equal to $\binom{|\beta|}{k}$. Thus

[TABLE]

To prove irreducibility, we observe that, if $\bm{0}\in S$ denotes a zero binary string, then $g(\alpha,\bm{0})\neq 0$ and $g(\bm{0},\alpha)\neq 0$ for every $\alpha\in S$. Then [21, § 3.11, Exercise 12a] implies that $G_{\ell}^{T}$ is irreducible.

Since $G_{\ell}^{T}$ is stochastic, its largest eigenvalue is equal to $1$ [21, § 8.5, p. 156]. Since $G_{\ell}^{T}$ is irreducible, [21, Theorem 8.2] implies that $1$ is a simple eigenvalue.

To prove the second part of the lemma, we fix $1\leqslant i<j\leqslant\ell$. We will show that for every $\beta\in S_{i,j}$

[TABLE]

Indeed, let $\gamma$ be a binary string with all zeroes and one at the $i$-th position. Then, since $g(\gamma,s(\beta))\geqslant\frac{1}{2^{|\beta|+2}}\geqslant\frac{1}{2^{\ell+2}}$, we have

[TABLE]

Inequality (3) implies that there exists a matrix $C_{\ell;i,j}$ with nonnegative entries such that $\frac{2^{\ell+2}}{2^{\ell+2}-1}(G_{\ell;i,j}+C_{\ell;i,j})^{T}$ is stochastic.

Since $\bm{0}\in S_{i,j}$, the same argument as in the proof of the first part of the lemma shows that $\frac{2^{\ell+2}}{2^{\ell+2}-1}(G_{\ell;i,j}+C_{\ell;i,j})^{T}$ is stochastic, and has exactly one of the eigenvalues being equal to $1$. ∎

Corollary A.1**.**

Let $P_{\ell}(t)$ be the charactersitic polynomial of $G_{\ell}$. Then, for every $\ell>1$, $|P_{\ell}^{\prime}(1)|<1$.

Notation A.2**.**

Fix a positive integer $n$. For vectors $\bm{a}=(a_{1},\ldots,a_{n})\in\mathbb{Z}_{\geqslant 0}^{n}$ and $\bm{b}=(b_{1},\ldots,b_{n})\in\mathbb{Z}_{\geqslant 0}^{n}$, we denote

[TABLE]

Lemma A.2**.**

Let $A$ be an $s\times s$ stochastic matrix with only one of the eigenvalues being one. We set

[TABLE]

Let $P_{A}(t)$ be the characteristic polynomial of $A$. Then $\lim\limits_{n\to\infty}C(A)_{n}=\frac{1}{P_{A}^{\prime}(1)}$.

Proof.

We recall that the Lambert $W$ function [2] is the principal branch of the inverse of $xe^{x}$. We will use the notation $y(z)=-W(-z)$ from [7] so that $y(z)=ze^{y(z)}$. Function $y(z)$ has a singularity of the square-root type at $z=1/e$ and has the following expansion around this point (see [7, p. 107])

[TABLE]

From this, we obtain

[TABLE]

The main result of [1] implies that, for every complex $s\times s$ matrix $A$, we have

[TABLE]

Since $A^{T}$ is stochastic, we have $\sum\limits_{i=1}^{n}a_{i,j}=1$, so

[TABLE]

If we perform a substitution $x=y(z)$ and use the definition of the Lambert $W$ function, we obtain

[TABLE]

From this, we obtain

[TABLE]

$F(z)$ can be rewritten as

[TABLE]

Finding the asymptotic behavior of the Taylor coefficients of $F(z)$ would yield an asymptotic for $C(A)_{n}$. We will do this using singularity analysis [6, Chapter VI] (similarly to [7, Theorem 2]). Since $|y(z)|<1$ for $|z|<1/e$ (see [2, Fig. 1]) and all roots of $P_{A}$ lie in the unit circle due to the stochasticity of $A$, $\frac{1}{e}$ is the singularity of $F(z)$ with the smallest absolute value. Due to Lemma A.1, $P_{A}(t)=(1-t)Q_{A}(t)$, where $Q_{A}(1)\neq 0$. Using (5), we obtain the following expansion of $F(z)$ around $1/e$:

[TABLE]

Singularity analysis [6, Corollary VI.1] implies that

[TABLE]

Using Stirling’s formula, we get

[TABLE]

Using $P_{A}^{\prime}=-Q_{A}^{\prime}+(1-t)Q_{\ell}$, we deduce $P_{A}^{\prime}(1)=-Q_{A}(1)$, and this finishes the proof. ∎

Lemma A.3**.**

On the set of all Boolean networks with $n$ states consider two probability distributions:

1. (A)

all the networks with canalizing depth one have the same probability, all others have probability zero; 2. (B)

the probability assigned to each network is proportional to the product of the number of canalizing variables of the functions defining this network.

We fix a positive integer $\ell$. By $A_{\ell,n}$ and $B_{\ell,n}$ we denote the average number of attractors of length $\ell$ in a random Boolean network with $n$ states with respect to distributions (A) and (B), respectively. Then

[TABLE]

Example A.1**.**

We will illustrate the distribution (B) by an example. Consider three following networks with two states:

[TABLE]

Since the canalizing depth of $x_{1}+x_{2}$ is zero, $P_{B}(\mathbf{f}_{1})$, the probability of $\mathbf{f}_{1}$ with respect to $B$, is zero. Since the canalizing depths of $x_{1}x_{2}$ and $x_{1}$ are $2$ and $1$, respectively, the ratio $P_{B}(\mathbf{f}_{2})/P_{B}(\mathbf{f}_{3})$ is equal to $\frac{2\cdot 1}{2\cdot 2}=1/2$.

Proof.

Let $F_{n}$ and $F^{\ast}_{n}$ be the number of Boolean functions in $n$ variables with of canalizing depth exactly one and more than one, respectively. We will use the following bounds:

$F^{\ast}_{n}\leqslant n^{2}\cdot 4\cdot 4\cdot 2^{2^{n-2}}$. We look term-by-term. There are at most $n^{2}$ ways to choose first and second canalizing variables. There are at most $4$ choices for the canalizing outputs and at most $4$ choices for canalizing values for these two variables. There are at most $2^{2^{n-2}}$ core functions, since that is all possible functions, which may or may not be canalizing. Since redundant arrangements of canalizing variables are not accounted for, this must overcount. 2. 2.

$F_{n}\geqslant 2^{2^{n-1}}-(n-1)\cdot 2\cdot 2\cdot 2^{2^{n-2}}$. This is a lower bound for the number of non-canalizing core function in $n-1$ variables because $(n-1)\cdot 2\cdot 2\cdot 2^{2^{n-2}}$ is an upper bound on the number of canalizing functions in $n-1$ variables (obtained in the same way as the bound above).

We also introduce

[TABLE]

For $X$ being (A) or (B) and positive integer $n$, let $P_{X,n}$ denote the probability (it is always the same) of choosing a network from distribution $X$ with all functions being of depth exactly one. Let $P_{n}^{\ast}$ be the maximal probability of choosing a network from (B) with at least one function being of depth more than one, respectively. By $S_{n}$ and $S_{n}^{\ast}$ we denote the total number of attractors of length $\ell$ in networks with all functions being of depth exactly one and with at least one function being of depth more than one, respectively.

The statement of the lemma is equivalent to the statement that

[TABLE]

Using the notation introduced above, we can bound $A_{\ell,n}-B_{\ell,n}$ as

[TABLE]

We set $U_{n}:=S_{n}(P_{n,A}-P_{n,B})$ and $V_{n}:=P_{n}^{\ast}S_{n}^{\ast}$. Then (12) would follow from $\lim\limits_{n\to\infty}U_{n}=0$ and $\lim\limits_{n\to\infty}V_{n}=0$, so we will prove these two equalities.

Since any network has at most $2^{n}$ attractors of length $\ell$, $S_{n}\leqslant 2^{n}F_{n}^{n}$. Since the total sum of the products of canalizing depths over all Boolean networks does not exceed $(F_{n}+nF_{n}^{\ast})^{n}$, we have $P_{n,B}\geqslant\frac{1}{(F_{n}+nF_{n}^{\ast})^{n}}$. Since $P_{n,A}=\frac{1}{F_{n}^{n}}$, we have

[TABLE]

(11) implies that $nR_{n}<1$ for large enough $n$. Hence, for large enough $n$, we have

[TABLE]

By similar arguments, $P_{n}^{*}\leqslant\frac{n^{n}}{F_{n}^{n}}$ and $S_{n}^{*}\leqslant 2^{n}n(F_{n}+F_{n}^{\ast})^{n-1}F_{n}^{\ast}$ so:

[TABLE]

Since $R_{n}<1$ for large enough $n$, using (11), we have

[TABLE]

Remark A.1**.**

The proof of Lemma A.3 will be valid if we replace distribution (B) with any other distribution $(C)$ such that, for every Boolean network $\mathbf{f}=(f_{1},\ldots,f_{n})$,

• •

if at least one of $f_{i}$’s is non-canalizing, $P_{C}(\mathbf{f})=0$;

• •

there exists a constant $P_{n,C}$ such that, if the canalizing depth of every $f_{i}$ is one, then $P_{C}(\mathbf{f})=P_{n,C}$;

• •

we have $\frac{P_{C}(\mathbf{f})}{P_{n,C}}\leqslant\frac{P_{B}(\mathbf{f})}{P_{n,B}}$ (using notation from the proof of Lemma A.3).

The above properties hold, for example, for the following distribution

1. (C)

all the networks defined by canalizing functions have the same probability, all others have probability zero.

Using this distribution instead of (A), we see that Theorem 7.1 holds also for a random Boolean network defined by canalizing functions.

Lemma A.4**.**

We will use Notation A.1 and notation from Lemma A.2. Then, for every positive integers $\ell$ and $n$, we have

[TABLE]

Proof.

We fix $n$. Consider a tuple $\bm{X}=(X_{1},\ldots,X_{\ell})$ of $\ell$ distinct elements of $\{0,1\}^{n}$. For $1\leqslant i\leqslant n$, we denote $\bm{X}_{i}:=(X_{1,i},\ldots,X_{n,i})$. For $\alpha\in S$, let

[TABLE]

Then $\sum\limits_{\alpha\in S}n_{\alpha}=n$. First, we will show that

[TABLE]

where $\mathbf{n}=(n_{0},n_{1},\ldots,n_{2^{\ell}-1})$.

To prove (15), we will use that the functions $f_{i}$ ($i=1,\ldots,n$) in the network are chosen independently to decompose the left-hand side as

[TABLE]

where we use notation $X_{n+1}=X_{1}$ and the probability of each Boolean function to be chosen is assumed to be proportional to the number of its canalizing variables. We show that, for every $1\leqslant i\leqslant n$,

[TABLE]

Then (15) would follow from multiplying (16) for all $i$. To prove (16), without loss of generality, we consider $i=1$. Consider a set

[TABLE]

with a uniform probability distribution $P_{\Omega}$. Observe that for a function $f$ with canalizing variables $x_{k_{1}}$, $\ldots$, $x_{k_{s}}$, we have

[TABLE]

If we can show that, for every $1\leqslant k\leqslant n$,

[TABLE]

then (16) would follow by summing up (17) over all $k$ and using the law of total probability.

We consider one of the canalizing variables of $f$, say, $x_{k}$. Let $c$ be the canalizing value of $x_{k_{1}}$, and let $v$ be the value taken by $f$ when $x_{k_{1}}=c$. Then $(c,v)\in\{0,1\}^{2}$, and all these four cases have the same probability due to the symmetry. As $g(\alpha,s(\beta))=\frac{1}{4}(h(\alpha,\beta)+h(\bar{\alpha},\beta)+h(\alpha,\bar{\beta})+h(\bar{\alpha},\bar{\beta}))$, it is sufficient to show that

[TABLE]

and then sum for all $(c,v)\in\{0,1\}^{2}$.

To prove (18), consider any $j$, say $j=1$. There are then 4 cases for the values of $X_{1,k}$ and $X_{2,1}$.

$X_{1,k}=1$ and $X_{2,1}$ is [math] or $1$. With probability $\frac{1}{2}$, we have $f(X_{1})=X_{2,1}$. This is true due to symmetry, as for any $f_{1}$ which takes on the value $w$ at $X_{1}$, we can produce another function $g$ that is equal to [math] if $X_{1,k}=0$ and $\bar{f}_{1}$ if $X_{1,k}=1$. Then $g(X_{1})=\bar{w}$. 2. 2.

$X_{1,k}=0$ and $X_{2,1}=1$. Since $X_{1,k}=c$, the probability of $f(X_{1})=X_{2,1}\neq v=0$ is zero. 3. 3.

$X_{1,k}=X_{2,1}=0$. Since $X_{1,k}=c$ and $X_{2,1}=v$, the canalization property implies that $f(X_{1})=X_{2,1}$ with probability one.

The only case in which $\bm{X}_{1}\vee s(\bm{X}_{k})\neq s(\bm{X}_{k})$ is where there is at least one $j$ such that case 2 is realized. In this case, the probability in the left-hand side of (18) will be zero. Otherwise, each occurrence of case 1 will multiply the total probability by $\frac{1}{2}$ and each occurrence of case 3 will multiply the total probability by $1$. Thus, we show that the left-hand side of (18) is indeed equal to $h(\bm{X}_{1},s(\bm{X}_{k}))$. This finishes the proof of (15).

To finish the proof of the lemma, we set

[TABLE]

Summing (15) over all $\ell$-tuples $(X_{1},\ldots,X_{\ell})$ of distinct elements of $\{0,1\}^{n}$, we obtain (see (4))

[TABLE]

On the other hand, if $\mathbf{n}$ is supported on one some $S_{i,j}$, then $G_{\ell}\mathbf{n}=G_{\ell;i,j}\mathbf{n}|_{S_{i,j}}$, where $\mathbf{n}|_{S_{i,j}}$ denotes the restriction of $\mathbf{n}$ on the coordinates from $S_{i,j}$. This implies that

[TABLE]

This finishes the proof of the lemma. ∎

Proof of Theorem 7.1.

We fix positive integer $\ell$. In the notation of Lemma A.3, we have $A_{\ell}=\lim\limits_{n\to\infty}A_{\ell,n}$. Lemma A.3 implies that $A_{\ell}=\lim\limits_{n\to\infty}B_{\ell,n}$. We fix any $1\leqslant i<j\leqslant\ell$, and let $C_{\ell;i,j}$ be the matrix given by Lemma A.1. We set $M:=\frac{2^{\ell+2}}{2^{\ell+2}-1}(G_{\ell;i,j}+C_{\ell;i,j})$. Then

[TABLE]

Lemma A.2 implies that $\lim\limits_{n\to\infty}C(M)_{n}$ is finite, thus we have that $\lim\limits_{n\to\infty}C(G_{\ell;i,j})_{n}=0$. We finish the proof of the theorem by considering the limit of (14) and applying Lemma A.2 to $G_{\ell}$. ∎

Appendix B: Complexity analysis

Proposition B.1**.**

Complexity of Algorithm 3 is $\mathrm{O}(k^{3})$.

Proof.

First we show the complexity of a single run the algorithm, i.e., not taking into account the recursive call, is $\mathrm{O}(k^{2})$.

First we show the complexity of a single run the algorithm, i.e., not taking into account the recursive call, is $\mathrm{O}(k^{2})$.

First we show the complexity of a single run the algorithm, i.e., not taking into account the recursive call, is $\mathrm{O}(k^{2})$.

First we show the complexity of a single run the algorithm, i.e., not taking into account the recursive call, is $\mathrm{O}(k^{2})$.

First we show the complexity of a single run the algorithm, i.e., not taking into account the recursive call, is $\mathrm{O}(k^{2})$. Since the first $k$ rows of the Pascal’s triangle can be precomputed in $\mathrm{O}(k^{2})$, the complexity of step 1 is also $\mathrm{O}(k^{2})$. Similarly, the complexity of step 3 is $\mathrm{O}(k^{2})$. It remains to observe that step 2 takes $\mathrm{O}(1)$ and step 4 takes $\mathrm{O}(k^{2})$ (indeed, selecting a subset of size $j$ amounts to selecting and removing $j$ indices). In total, we obtain $\mathrm{O}(k^{2})$.

The depth of the recursion calls is at most $k$. Since the complexity of each single call is $O(k^{2})$, so the total complexity is $O(k^{3})$. ∎

Lemma B.1**.**

The average complexity of the algorithm for generating a function in $n>0$ variables which is either $1$ or noncanalizing described in Remark 5.2 is $\mathrm{O}(n2^{n})$.

Proof.

[25, p. 116] implies that the proportion of functions which are canalizing in $n$ variables is bounded from above by $\frac{4n}{2^{2^{n-1}}}$. Note that [25] considers constant functions canalizing which we do not. Thus the probability $P_{n}$ of choosing a function which is either $1$ or noncanalizing is bounded from above by

[TABLE]

This bound is less than $\frac{3}{4}$ for all values of $n$ except 1 and 2, but we can compute directly that $P_{1}=\frac{3}{4}$ and $P_{2}=\frac{13}{16}$. Therefore, the number of times the generation of a function needs to be repeated averages to $\frac{1}{1-P_{n}}$, which does not exceed $4$, so the average complexity of the whole procedure is the same as of a single generation step.

The complexity of a single step consists of generating a random function (which is $\mathrm{O}(2^{n})$) and checking whether it is canalizing or not. We perform this check by running linearly through the table for each variable, so the complexity is $\mathrm{O}(n2^{n})$ time. Thus, the total complexity is indeed $\mathrm{O}(n2^{n})$. ∎

Lemma B.2**.**

There is a constant $c<1$ such that the probability that a function generated in steps 1 and 2 of Algorithm 2 does not satisfy one of the of the conditions (E1) or (E2) is bounded by $c$ for every $n$.

Proof.

Notice that

[TABLE]

We will show that there is a constant $c<1$ such that $P(\text{\ref{e1} is false}\mid r\neq 1)$ and $P(\text{\ref{e2} is false}\mid r=1)$ do not exceed $c$.

• •

$P(\text{\ref{e1} is false}\mid r\neq 1)$. The probability of having $k_{r}=1$ (the only possible $k_{r}<2$) is just the proportion of ordered partitions with a single element at the end. We can construct all of these by picking an element and then picking a partition of the remaining elements, so this creates $k\cdot p_{k-1}$ possibilities. Thus the probability this occurring is $\dfrac{kp_{k-1}}{p_{k}}$. [8, Eq. (5)] implies that this approaches $\ln(2)<1$ as $n$ goes to infinity. Thus, there exists such $c$.

• •

$P(\text{\ref{e2} is false}\mid r=1)$. The probability of ever picking $b=1$ is just $\frac{1}{2}$, so we can take $c=\frac{1}{2}$.

∎

Proposition B.2**.**

Complexity of Algorithm 2 is $\mathrm{O}(n2^{n})$.

Proof.

Lemma B.2 implies that the average number of reruns in step 3 is constant. Thus, the complexity of the algorithm is the same as of a single run.

Proposition B.1 and Lemma B.1 imply that the complexity of step 1 is $\mathrm{O}(k^{3}+(n-k)2^{n-k})$. Step 2 generates a truth table for the function. There are $2^{n}$ input-output pairs, and computing the function takes at most $k$ steps, so this is $\mathrm{O}(k2^{n})$. In step 3, the conditions (E1) or (E2) are verified in $\mathrm{O}(2^{n})$ time.

Summing everything, we obtain $\mathrm{O}(k^{3}+(n-k)2^{n-k}+k2^{n})=\mathrm{O}(n2^{n})$ ∎