Ordered dynamics in biased and cooperative Boolean networks

This paper contributes to the theoretical analysis of the qualitative behavior of two types of Boolean networks: biased and cooperative ones. A Boolean network is biased if at least a speciﬁed fraction of its regulatory functions returns one Boolean value more often than the other and is cooperative if there are no negative interactions between the variables. We prove nontrivial upper bounds on the maximum length of periodic orbits in such networks under the assumption that the maximum number of inputs and outputs per node is a ﬁxed constant r . For the case of n -dimensional networks with r = 2 in which only AND and OR are allowed, we ﬁnd an upper bound of 10 n /4 , which is asymptotically optimal in view of previously published counterexamples. The theoretical results are supplemented by simulations of the generic behavior of cooperative networks which indicate that for large indegrees, trajectories tend to converge rapidly towards a steady state or a small periodic orbit. The latter starkly contrasts with the behavior of random arbitrary Boolean networks. MSC: 05A15; 06A07; 34C12; 39A33; 94C10; 92C42


Introduction
Mathematical models of biological problems have become increasingly sophisticated with the advent of new quantitative techniques. Models involving dozens of proteins at the cell level abound in the literature [, ], mostly created using continuous methods. An interesting counterpart to continuous systems are Boolean networks, characterized by discrete time dynamics and the existence of only two possible states ( and ) for each variable at any given time []. In the face of parameter uncertainty, modeling using Boolean networks can produce qualitative predictions at the right level of detail for many applications.
This paper continues a line of theoretical research exploring the qualitative behavior of so-called cooperative Boolean networks. This concept corresponds to systems that have exclusively positive interactions among their variables, and in particular are monotone, which means that they have only positive feedback loops. Such networks have been proposed as important tools for gaining insight into the dynamics of gene regulatory and other biological networks []. In the case of systems with at most two inputs for each variable, a cooperative Boolean network can only have constant, AND, OR and COPY regulatory functions. The continuous counterpart of this definition has been studied extensively [], and under mild conditions the generic trajectory of a continuous cooperative system is guaranteed to converge towards the set of equilibria. http://www.advancesindifferenceequations.com/content/2013/1/313 Much of the literature on the dynamics of Boolean networks has focused on the study of so-called random Boolean networks or RBNs, where an n-dimensional Boolean network is randomly drawn from the class of all such networks with a maximum of r inputs per node. For r > , the trajectory of a randomly chosen initial state tends to reach a periodic orbit whose size scales exponentially in n. Such long orbits are considered a hallmark of chaotic dynamics. Another hallmark of chaotic dynamics is sensitivity to initial conditions. This can be conceptualized in a variety of ways, for example, in terms of p-instability, as defined at the end of Section . A third hallmark of chaos are relatively few or no eventually frozen nodes, that is, nodes whose state will change only finitely often during the trajectory. These hallmarks usually, but not always, go together. It is a matter of judgement which hallmarks are deemed necessary for considering the dynamics truly chaotic; see [] for examples and further discussion on this point. In terms of these hallmarks, the dynamics of RBNs tends to become, on average, more chaotic as r is increased or if the average bias of the regulatory functions, defined as the difference between . and the proportion of input vectors on which it takes the Boolean value , decreases. For surveys of these and related results, see [-].
In Section  we briefly explore the analogue for cooperative Boolean networks. That is, for r ∈ {, , , } we simulate the trajectories of randomly chosen initial states in Boolean networks that are randomly chosen from the class of all cooperative Boolean networks with at most r inputs per node. We find that for r >  trajectories tend to reach a steady state or a very small periodic orbit after a few steps, and this pattern appears to be stronger for larger r and n. We also give an argument in support of a conjecture that if there is no restriction on the number of inputs, for sufficiently large n, the vast majority of the trajectories will reach either of the steady-state vectors (, . . . , ) or (, . . . , ). Thus, on average, random cooperative Boolean networks tend to have much shorter periodic orbits than corresponding RBNs, and in contrast to the latter, their dynamics appears to become less chaotic, at least in terms of lengths of periodic orbits that are reached from random initial conditions and the proportion of eventually frozen nodes, as the number r of allowed inputs per node increases.
These observations about average behavior naturally lead to the question whether there exist provable and nontrivial upper bounds on the length of periodic orbits, or at least on the median length of periodic orbits, that will be reached from randomly drawn initial conditions. Since the state space of an n-dimensional Boolean network has size  n , a bound might be considered nontrivial if it does not exceed c n for some constant c < .
If only the maximal number of inputs and outputs per node is restricted, even to , then no such bound can be derived; nice examples are already given in [, ]. On the other hand, very stringent upper bounds can be derived for some very special classes of networks, such as the XOR Boolean networks [, ], or the class of networks that are allowed to use only AND-operators or OR-operators, but not both [].
For cooperative Boolean networks, the situation is more subtle. The maximum possible size of any periodic orbit in such networks is known to be n n , which by Stirling's formula scales like  n+ √ π n , faster than c n for any c <  [, , ]. It has previously been shown that if in addition to requiring cooperativity we restrict the maximum number r of inputs and the number of outputs per node, then nontrivial upper bounds on the (maximal or even median) lengths of periodic orbits still do not exist, even for r =  (see Theorem  below, which summarizes some results from [, -]). However, for c >  n/ , the constructions http://www.advancesindifferenceequations.com/content/2013/1/313 used in the proof of this theorem give systems in which the vast majority of variables take just one input, that is, their regulatory function is COPY. Note that while COPY is an unbiased function, both AND and OR are biased. For example, binary AND returns  for only one quarter of its inputs. In other words, for every  < α < , most constructions used for the proof of Theorem  give Boolean systems in which a proportion of more than (α) of all regulatory functions are unbiased.
In contrast, suitable assumptions on the amount of bias do imply nontrivial upper bounds on the possible lengths of periodic orbits. Our main result, Theorem , shows that upper bounds of the form c n for some c <  on the lengths of periodic orbits exist under the assumption of fixed bounds on the minimal proportion of biased regulatory functions and the maximal number of inputs and outputs per node. In particular, this result applies to cooperative Boolean networks with a fixed minimal proportion of binary AND or binary OR regulatory functions. By Theorem (ii), an upper bound on the number of inputs alone is not sufficient in this type of result. For r = , if the proportion of biased regulatory functions is , then we can take c =  / , which is optimal in view of Theorem (iii).
Even more stringent bounds apply in some cases if we additionally require a certain amount of sensitivity to initial conditions. Theorem  gives upper bounds for the case when all regulatory functions take exactly two inputs and are biased. It implies upper bounds on the lengths of the form c n for any c > √  in systems of this form that are punstable for p sufficiently close to . By Theorem . of [], the latter bound is asymptotically optimal. Note that in contrast to Theorem , no upper bound on the number of outputs per node is required in Theorem .
Preliminary versions of the proofs given in Sections - were posted in the preprint []. One purpose of this paper is to make streamlined and more easily readable versions of these arguments available in peer-reviewed journal form. The material in Section  is entirely new.

Definitions
An n-dimensional Boolean network (g, ) consists of a state space = {, } n , together with a time evolution map g : → . Each individual component g k : → {, } is called the kth regulatory function, and a trajectory of the system is a function s : N → such that s(t + ) = g(s(t)) for every t ∈ N. Since the state space is finite, each trajectory must eventually become periodic. The set of states in the periodic part of the trajectory will be called a periodic orbit. Periodic orbits that are reached from a given initial condition are often called attractors or limit cycles of the trajectory in the literature. The basin of attraction of a given periodic orbit is the set of all initial conditions from which the orbit is reached.
Define the partial order ≤ on by x ≤ y if x i ≤ y i for every i. A regulatory function g k is cooperative if x ≤ y implies g k (x) ≤ g k (y). A Boolean network is cooperative if each of its regulatory functions is cooperative, i.e., if x ≤ y implies g(x) ≤ g(y). An equivalent definition can be obtained by requiring that each regulatory function can be written as a composition of Boolean functions that use only the operators AND and OR with one or more inputs. Since the negation operator is not needed, cooperative Boolean networks are exactly the ones without negative interactions.
We define a (b, r)-Boolean network as a system in which each regulatory function g k actually depends on at most r inputs, and each variable has at most b outputs (i.e., it affects the http://www.advancesindifferenceequations.com/content/2013/1/313 dynamics of at most b other variables). If r = , we call the system quadratic; a (, )-system is called bi-quadratic. A regulatory function that depends on only one variable is called monic; a non-monic quadratic regulatory function is called strictly quadratic. A Boolean network with only quadratic regulatory functions will be called a strictly quadratic network. Notice that in a strictly quadratic, bi-quadratic network every variable must have exactly two inputs and two outputs.
An n-dimensional Boolean network is said to be c-chaotic, for c < , if it has a periodic orbit of length at least c n . If a periodic orbit of this length is reached from a proportion of at least p initial conditions, then the system is called p-c-chaotic.
The bias of a regulatory function g k is the fraction of input vectors for which the function outputs . If = ., we say that g k is unbiased; if = ., we say that g k is biased. For  < α < , we say that a Boolean system is α-biased if a proportion of at least α of its regulatory functions consists of biased functions. Note that this last definition neither distinguishes between bias towards  or  nor all by itself implies how strong the bias is.
A Boolean system is p-unstable if a random single-bit flip in a randomly chosen initial state moves the trajectory into the basin of attraction of a different periodic orbit with probability at least p. Let us call a Boolean system p--unstable if a random single-bit flip in a randomly chosen initial state results in a different state after one updating step with probability at least p (even if the new and old trajectories might merge after several updating steps). Thus every p-unstable Boolean system must be p--unstable, but not vice versa.

Bounds on the length of periodic orbits
For easier reference, let us quote the following result.
Theorem  Let c, c  , p be constants with  < p <  < c <  and  < c  <  / . Then, for all sufficiently large n, there exist n-dimensional cooperative Boolean networks that are, respectively: (i) bi-quadratic and p-c-chaotic, (ii) strictly quadratic and p-c-chaotic, (iii) strictly quadratic, bi-quadratic and p-c  -chaotic. Theorem  Let  ≥ α > , and let b, r be positive integers. Then there exists a positive constant c(α, b, r) <  such that for every c > c(α, b, r) and sufficiently large n, there is no cchaotic, n-dimensional, α-biased (b, r)-Boolean system. In particular, for every c > c(α, , r) and sufficiently large n, every c-chaotic, n-dimensional cooperative (, r)-Boolean system uses COPY for more than (α)n of its regulatory functions.
Proof We will prove only the first part of the theorem; the second part then follows from the observation that the only nonconstant cooperative unbiased Boolean function that takes at most two inputs is the COPY function. http://www.advancesindifferenceequations.com/content/2013/1/313 Fix α, b, r as in the assumptions, and let ( , g) be an n-dimensional (b, r)-Boolean system. Assume that at least αn of the regulatory functions g k are biased.
Let S = (s : ∈ [L]) be a sequence of (not necessarily pairwise distinct) elements of = {, } n . If the elements of S happen to be pairwise distinct, then we will speak of S being a subset of . Let i ∈ [n] and consider the ratio More generally, let v ∈ [n] and σ : with i  < · · · < i v , we define ratios ξ σ I (S) as follows: One can think of ζ * I (S) ≥  as a measure of the bias of the sequence S on the set I. Notice that ζ * I (S) =  implies ζ * J (S) =  for all nonempty subsets J of I, and in particular that ζ i (S) = . for all i ∈ I.
If S is randomly chosen, then for i, I as above the expected values are E(ζ i (S)) = . and E(ξ σ I (S)) =  -|I| . Thus if S is sufficiently large, one should be able to choose a subset S * ⊆ S with ζ i (S * ) = ., ζ * I (S * ) =  and |S * | |S| arbitrarily close to . The following lemma shows that periodic orbits S in ( , g) are far from random in this respect.
Lemma  There exists τ >  such that for all ( , g) as above, all biased g k , all periodic orbits S of ( , g) and all subsets S * with |S * | ≥ (τ )|S|, it holds that ζ * I (S * ) =  for I = {k} ∪ dom(g k ).
Proof Let ε >  be such that for each of the finitely many biased Boolean functions with Choose g k as in the assumption and assume wlog that ≥ . + ε; the proof in the case when ≤ .ε is symmetric. Let J be the domain of g k , and let I = {k} ∪ J.
Consider a subset S * of a periodic orbit S with |S * | ≥ (τ )|S|. As we mentioned in the discussion that follows (), if ζ * J (S * ) = , then also ζ * I (S * ) =  and there is nothing to prove. So assume that ζ * J (S * ) = . Then all possible input vectors for g k appear in exactly equal proportion in S * , and it follows that |{s ∈ S * : g k (s) = }| ≥ |S * | ≥ (τ )|S| . The following inequalities show that in this case ζ k (S * ) > ., which implies ζ * I (S * ) = .
Lemma  There exists a constant β >  such that for τ as in Lemma , no periodic orbit of ( , g) can be β-τ -balanced.
A recursive construction allows us to find αn/(r + b(r -) + ) pairwise disjoint sets I k corresponding to biased functions g k . That is, for β := α r+b(r-)+ , there exists a family D of pairwise disjoint subsets I ⊂ [n] of size  ≤ |I| ≤ r +  each, with |D| ≥ βn, such that each set in D is of the form I k for some biased g k . By Lemma , this family witnesses that no periodic orbit of ( , g) can be β-τ -balanced. Now Theorem  becomes a consequence of the following result.
Proof Let β, τ be as in the assumptions and assume throughout this argument that n is a sufficiently large positive integer. Let ε =  -r- τ . Consider an arbitrary family D of pairwise disjoint subsets I of [n] with |D| ≥ βn and  ≤ |I| ≤ r +  for each I ∈ D.
The proof proceeds in three stages. First we will define, for  < c < , random variables η c . Next we will show that if there exists any subset S + ⊆ {, } n of size c n that is not β-τbalanced for D, then In the third stage we will complete the proof by showing that for c sufficiently close to , inequality () fails.
Stage : Defining η c Fix c with  < c <  and consider the space of all sequences S = (s : ∈ [L]) of randomly and independently (with replacement) chosen states in {, } n of length L = c n .
Consider [v] . Note that v ∈ [r + ]. We will treat ξ σ I = ξ σ I (S) and ζ * I = ζ * I (S) as random variables and usually suppress their dependence on S in our notation. We define Stage : Nonbalanced sets imply () Consider the following procedure for randomly drawing a subset S † of {, } n of size c n : First draw a sequence S as in Stage . Then the set S red = {s : ∈ [L]} associated to the sequence S has size c n -N , where N = N(S) is a random variable that represents the number of s in S that duplicate an earlier term in the sequence. Form S † = S † (S) by http://www.advancesindifferenceequations.com/content/2013/1/313 adding a randomly chosen subset of {, } n \S red of size N to S red . By symmetry, this gives the uniform distribution, and moreover, for any given M with  ≤ M < c n , the conditional probability of any potential outcome S + is still uniform, that is, The expected value of N can be estimated as In particular, we can assume E(N) < ε  c n . Now it follows from () and Markov's inequality Now assume that there is any subset S + ⊂ {, } n of size c n that is not β-τ -balanced for D, and suppose that S exists for which N(S) ≤ εc n and S † (S) = S + . Consider . These sets are pairwise disjoint, and their union S * satisfies ζ * I (S * ) = . Moreover, |S + \S * | ≤  r+ εc n = τ c n . This would contradict the assumption that S + was not β-τ -balanced for D. Thus we must have and () follows from () and ().
Stage : Inequality () fails for sufficiently large c . Note that ξ σ I is the mean of c n independent random variables that take values in {, }. Thus Hoeffding's inequality [] applies, and we get where λ = e -ε  < , which in turn implies that P(ζ * I ≥ ε) ≤  v λ c n . Since the random variables ζ * I for I ∈ D are independent, it follows that Thus for c > λ β and n sufficiently large, estimates () and () imply as required.
4 More stringent bounds for the case b = r = 2, α = 1 Let c(α, b, r) denote the smallest real number for which the conclusion of Theorem  holds in the class of cooperative Boolean networks. Theorem  states that c(α, b, r) > , but gives no estimates of the magnitude of the difference. Some lower bounds for the differences were extracted from the proof of the theorem in Appendix C of preprint []; they appear to substantially underestimate the actual difference. For example, the proof of Theorem  shows only that c(, , ) ≤ .. On the other hand, Corollary  of the preprint [] gives an upper bound of c(α, , ) ≤  (-α)/ , which is more stringent if α is sufficiently close to . In particular, since for c > c(, , ) no strictly quadratic, bi-quadratic Boolean network can be c-chaotic, this result implies together with Theorem (iii) that the latter estimate is optimal in this case and c(, , ) =  / ≈ .. We will prove this result here only for the special case of cooperative systems, which makes the argument more transparent. The proof for the general case given in [] uses essentially the same ideas.
Theorem  Consider a strictly quadratic, bi-quadratic, cooperative Boolean system ( , f ) of dimension n with a periodic orbit of size c n . Then c ≤  / .
Proof Let ( , f ) be as in the assumption, with (f  , . . . , f n ). Call a subset I ⊆ [n] closed if each f i for i ∈ I takes inputs only from I, and call I minimal closed if no proper subset of I is closed. Given the constraints on the indegree and outdegree of each node, it follows that all in-and outdegrees are equal to two, and that all outputs of elements from a minimal closed set I are also inside I. Therefore [n] is the union of pairwise disjoint minimal closed sets I , for ∈ [L], with  ≤ |I | ≤ n for all . Let us call a vector g = (g  , . . . , g k ) of Boolean functions on [k] a Boolean k-block if the corresponding Boolean system ({, } k , g) is strictly quadratic and bi-quadratic and [k] is minimal closed in this system. Let us define R(k) for k ≥  as the maximal size of the range of any Boolean k-block. Let (k) = (R(k)) /k , and let = max{ (k) :  ≤ k}.
Returning to ( , f ), we can see that L = |I | = n and for any C in the range of g, in particular, for any C that is a periodic orbit of the system, we must have Now the theorem is a consequence of the following lemma.
Proof For k =  we have () = √  <  / , so let us assume k > . Consider a Boolean k-block g = (g  , . . . , g k ). It is not hard to see that after a suitable renumbering of the input http://www.advancesindifferenceequations.com/content/2013/1/313 variables, we can assume wlog that g i takes inputs s i , s i+ for all i ∈ [k -] and g k takes inputs s  , s k . For ease of notation, let us identify k +  with . We call j ∈ [k] an o-node if g j = s j ∨ s j+ and an a-node if g j = s j ∧ s j+ . The proof of the lemma now splits into two cases.
Case : The sets of o-nodes and a-nodes are both nonempty In this case we can group the variables in [k] into one or more consecutive intervals I m = {j : i m < j ≤ i m+ } of length |I m | ≥  in such a way that for each m, either i m+ will be an o-node while all nodes j with i m < j ≤ i m+ are a-nodes or i m+ will be an a-node while all nodes j with i m < j ≤ i m+ are o-nodes. Wlog assume the latter; the proof in the other case is symmetric. Then the restriction of a vector in the range of g to the variables in I m cannot have an isolated  in any of the positions i m + , . . . , i m+ , that is, if g(s) j =  for some j with i m < j ≤ i m+ , then we must also have g(s) j =  for some j ∈ I m with |jj | = .
Let ψ(n) be the number of Boolean vectors of length n that do not contain isolated s, except possibly in their first position. The related function ϕ(n) that counts the number of Boolean vectors of length n that do not contain any isolated s whatsoever was studied in []. Any vector that has exactly one isolated , in its first position, must start with the sequence , and it follows that for n ≥ , It was reported in [] that ϕ(n) = cγ n + z  z n  +z z n  , where γ ≈ . is the real root of x  -x  + x - = , c ≈ ., and z  , z  are complex numbers with |z  | < . and |z  | < . Thus which together with () gives a crude upper bound Since ψ() =  <  / and γ <  / ≈ ., it follows that ψ(n) ≤  n/ for all n ≥ . Now we conclude that the range R of g satisfies Case : The set of o-nodes is empty or the set of a-nodes is empty Wlog assume that all nodes are o-nodes; the proof in the other case is symmetric. In this case, any vector in the range of g cannot have isolated s, except in the case where it takes the value  both on its first and last positions (so that these s are not isolated with respect to a circular topology on the indices). Let τ (n) denote the number of such vectors of length n. A vector in our target set could have no isolated s whatsoever, could be obtained by adding a  at the beginning of a sequence that has at most one isolated , in its starting position, and then appending the sequence  at the end, or reversing such a vector. This double-counts the vectors in the target set that start with  and end with http://www.advancesindifferenceequations.com/content/2013/1/313 . Thus in view of () we obtain τ (n) = ϕ(n) + ψ(n -)ϕ(n -) = ϕ(n) + ϕ(n -)ϕ(n -) + ϕ(n -).

p-Unstable strictly quadratic systems
In this section we derive even more stringent bounds on the maximal length of periodic orbits in quadratic systems that are α-biased (with sufficiently large α) than in the previous section, under the assumption that the system shows sufficiently large sensitivity to initial conditions. Let us first put this result into the context of previously constructed counterexamples.
For every positive integer n, there exists a -unstable, strictly quadratic, bi-quadratic cooperative Boolean system of dimension n (Proposition  of []). Moreover, for all positive constants p <  < c <  and all sufficiently large n, there exist n-dimensional cooperative Boolean networks that are simultaneously bi-quadratic, c-chaotic, p-unstable (Theorem  of []) and even p-c-chaotic (Theorem . of []). As long as c < √ , the analogue of the latter result holds for strictly quadratic, bi-quadratic cooperative Boolean networks (Theorem . of []). By choosing α =  and considering the limit as p → in the following theorem, we see that the upper bound on c in the latter counterexamples is optimal.
Theorem  Let  ≤ α ≤ , let c be a constant such that  √ . < c < , and let p > .α  + ln(.c)  ln . . Then no α-biased quadratic Boolean system can simultaneously be c-chaotic and p--unstable.
Proof Let c be as in the assumption and assume ( , g) is a c-chaotic strictly quadratic cooperative Boolean system of dimension n. A pair (j, j ) with j, j ∈ [n] will be called dominating if there are i, i ∈ [n] such that {s i , s i } is the set of input variables for both g j and g j , and at least one of the functions g j , g j is biased. Let I be the set of all s i with outdegree  that act as an input of a dominating pair. Let (j, j ) ∈ J, and let i, i be the corresponding indices of the inputs. Assume wlog that g j is biased and takes the value  on input vectors x, y, z ∈ {, } {i,i } . Then g j takes the same value on at least two of the vectors x, y, z, and it follows that the range of g j × g j has size at most . Since we assumed that there exists a periodic orbit of length at least c n , we must have  |J|/  n-|J| ≥ c n , and the lemma follows by taking logarithms. Now let I  denote the set of variables s i with outdegree zero, I  denote the set of variables s i with outdegree  that act as an input to a biased regulatory function, and I  the set of variables s i outside of I with outdegree  that act as an input to two biased regulatory functions.

Lemma 
Consider a random initial state s() and the state s * () obtained by flipping the value of the variable s i ().
Clearly, if s i ∈ I  , then s() = s * (), since the variable s i is not used at all to calculate the next state.
If s i ∈ I  , then there is exactly one s j for which s i acts as an input. Since g j is assumed biased, it must take another input i , and there exists a value k ∈ {, } such that the value of g j does not depend on the input s i as long as s i = k . The latter happens with probability .; thus with probability ≥ ., we will have s() = s * (). Now consider the case when s i ∈ I  . Then there are j = j such that s i acts as an input to both g j and g j , and g j , g j are both biased. Since i / ∈ I by definition of I  , the other input variables s i , s i of g j , g j are distinct and different from s i . Since g j , g j are both biased, there exist k , k ∈ {, } such that neither the value of g j nor the value of g j depends on the input s i as long as s i = k and s i = k . The latter happens with probability .; thus with probability ≥ ., we will have s() = s * (). Let us first assume for simplicity that ( , g) is strictly quadratic and bi-quadratic. Since the sum of all indegrees is equal to the sum of all outdegrees, the sets I  and I  are empty in this case. The set of variables that act as an input to at least one unbiased regulatory function has size less than (α)n, thus we must have |I  | > αnn -|I|. Then, by Lemma , Since with probability . a single-bit flip of a node in I  has no effect after one time step, inequality () implies that and the theorem follows for this special case. Now let us turn to the general case when ( , g) is only assumed quadratic. For each k, let I * k denote the set of variables with outdegree k. Then I *  = I  , but I  as defined above may be a proper subset of I *  ; similarly for I  . Since the sum of all outdegrees must equal the sum of all indegrees and the system was assumed to be quadratic, we must have We can deduce from () that which in turn implies As we mentioned above, I  = I *  . Moreover, since fewer than (α) of all regulatory functions are unbiased, we have Now it follows from () and () that From the above and () we conclude that and the theorem follows.

Generic behavior of cooperative networks
While our work so far has investigated the most extreme possible behavior (in terms of long periodic orbits), in this section we focus on the average behavior of cooperative Boolean networks. We first report some observations from simulation studies of cooperative Boolean networks that were randomly chosen from the class of all such networks with a specified maximum number r of inputs per node. Subsequently, we will formulate a conjecture on the expected behavior when the restriction on the number of inputs is removed and adduce some evidence for it. In order to select a random Boolean network, it is important to understand how to select random cooperative regulatory functions g : {, } r → {, } with at most r inputs. Such functions are in bijective correspondence with antichains aka Sperner families of P = {, } r , that is, subsets of P such that no two elements are comparable under the partial order ≤. For a given antichain A, one can define the function and one can verify that the assignment A → g A is a bijection. r = , as the corresponding Dedekind number is already ,,. We refrained from using ad hoc sampling methods for r ≥  as they do not guarantee a uniform distribution, and the inherent biases might invalidate the results. Given n, one antichain A k in {, } r was selected at random and independently for each variable. Moreover, the sets D k of inputs to each regulatory function were chosen randomly and independently in order to create a random Boolean network. The regulatory functions g k were defined as in (), with A = A k and x representing the restriction of the input vector to the domain D k of g k after a suitable permutation of the variables. A random initial condition was selected, and the trajectory was calculated for all t ≤ . This procedure was repeated for n =  and n = ,, as well as r = , , , , with  repetitions for each setting. The results of this simulation are displayed in Figure , using the norm |s(t)| = i s i (t) at each time step.
Notice that for r =  the norms of the trajectories tend to remain near n  . As r increases, the norms tend to converge towards  or n, indicating convergence of the trajectory towards (, . . . , ) or (, . . . , ) or to a small periodic orbit. The speed of convergence towards these steady states is similar for n =  and n = ,, but when r is increased from  to , we observe slightly faster convergence of the norm.
In order to estimate the distribution of the number of periodic orbits, their length, and the size of their basins of attraction, for each of r = , , ,  we sampled  networks. For each network we uniformly sampled  initial conditions, ran the system until it converged towards a periodic orbit, and measured the length of the orbit. We kept track of the number of trajectories that converged towards a given periodic orbit. We classified periodic orbits into four types as follows: the first type comprises the constant vectors (, . . . , ) and (, . . . , ), which are always steady states in the system, i.e., periodic orbits of length ; the second type comprises all other steady states; the third type comprises orbits of length  and ; and the fourth type comprises all orbits of length larger than .
In Figure (a)-(b) we display an estimated histogram for the total number of observed periodic orbits that were reached in each of the sampled networks, for r = , ,  (a) and r =  (b). Notice that this information is necessarily incomplete, since many periodic or-http://www.advancesindifferenceequations.com/content/2013/1/313 bits may not have been found. Nevertheless, one can clearly see that for r =  the average number of periodic orbits is significantly reduced. Figure (c) classifies networks by their indegree r and the type of a periodic orbit, and it displays the average number of periodic orbits across all sampled networks. Figure (d) complements this information by describing the average size of the basin of attraction of a given periodic orbit, measured as a percentage of the state space. We also measured the size distribution of the basins of attraction. We partitioned the interval [, ] into five equal intervals I  , . . . , I  , and for each such interval, we counted how many periodic orbits of a given type have a basin of attraction with estimated size within that interval. The results for each value of r are displayed in Figure (e) through (h).
Notice that even though there might be several long periodic orbits in a given network, the basins of attraction of such periodic orbits tend to be very small. Periodic orbits of length >  are reached from a significant fraction of initial conditions for r ∈ {, , }, but for r =  this fraction is already much smaller, and for r =  it becomes negligible. Overall, we found that most trajectories for all r that we investigated converge towards a steady state. For r ∈ {, }, almost all of these are different from the constant vectors (, . . . , ) and (, . . . , ), and for r =  most trajectories converge towards these two vectors. The largest basins of attraction on average comprise about half of the state space; for instance, if r =  the average steady state of type  attracts % of the trajectories. Since there are two such steady states, they combine to attract an average of % of all trajectories.
Close inspection of a subset of our simulations confirmed that the most complex behavior occurs for r = , with many trajectories reaching periodic orbits of the fourth type, that is, of length > , that appear to have very small basins of attraction. On the other hand, some networks for r =  appear to have a single steady state that attracts almost all initial conditions. http://www.advancesindifferenceequations.com/content/2013/1/313 The above simulations indicate that for sufficiently large r, random cooperative networks tend to have highly ordered dynamics, and that the amount of chaos, at least in terms of the lengths of observed periodic orbits and the percentage of eventually frozen nodes, is very low. The latter is exactly the opposite of what is known for RBNs without the assumption of cooperativity, and it calls for a theoretical explanation.
For fixed r, consider the family A(r) of all antichains in {, } r with the uniform distribution. Furthermore, for x ∈ {, } r , let p(r, x) denote the probability that x ≤ a for some a ∈ A, when A is randomly chosen from A(r). By symmetry, p(r, x) depends only on r and |x|, so we can extend it to a function defined on pairs of positive integers, with p(r, k) = p(r, x) for |x| = k. See Figure   While there exists a vast literature on antichains (see [] for an encyclopedic review), the function p(r, k) appears not to have been explicitly investigated. The literature indicates that the vast majority of antichains tend to mostly contain elements of size very close to r  ; in particular, the proofs of asymptotic formulas for the number of antichains on {, } r that were derived in [] and [] point in this direction. Thus it appears that if k(r) grows a little slower than r  , we should have lim r→∞ p(r, k(r)) = , and if k(r) grows a little faster than r  , we should have lim r→∞ p(r, k(r)) = . Let us call a function κ : Z + → Z + fast if p(r, r  +κ(r)) = o(r - ) and p(r, r  -κ(r)) = -o(r - ).

Conjecture  There exists a fast function κ such that κ(r) = o( √ r).
This conjecture appears to be highly plausible [], but may not immediately follow from known results. To see why it seems plausible, consider the related function p * (r, k) that is defined like p(r, k), but for antichains A that are randomly drawn from the family of all antichains on {, } r that consist of elements with identical norms. Then for k > r  the probability p * (r, k) is bounded from above by the probability that a randomly chosen antichain in the restricted class consists of elements with norm for some ≥ k. Note that for k > r  we have k ≤ r r/ r-k k+ . This in turn implies that if k > r+c √ r  for some c > , then (   ) http://www.advancesindifferenceequations.com/content/2013/1/313 and the analogue of Conjecture  for p * (r, k) with k > r follows since r r/ grows much faster than √ r. For k < r we get the same order-of-magnitude estimates, and the analogous result follows from the additional observation that if A is a randomly chosen antichain that consists of elements with norm > k, then for a given x with |x| = k, the probability that A contains an a with x ≤ a is equal to  - -( k ) .
It may be possible to extract a proof of the conjecture by carefully analyzing the arguments in [, ], but this is not easy and would go beyond the scope of this paper.
Our next result shows that Conjecture  has a striking consequence for the dynamics of random cooperative networks without restrictions on the number of inputs r per variable. The explicit calculations and simulations that we were able to perform for small values of r (see Figures  and ) add some credibility to this consequence of the conjecture.
Lemma  Assume that Conjecture  holds. Then the probability that the trajectory of a randomly chosen initial condition in a randomly chosen n-dimensional cooperative Boolean network reaches either the steady state vector (, . . . , ) or the steady state vector (, . . . , ) after one updating step approaches  as n → ∞.
Proof For fixed n, choosing a random cooperative network can be conceptualized as the same procedure as we have been following in our simulations, but with r = n. Thus for each k ∈ [n] we randomly and independently choose an antichain A k ∈ A(n), and then we define the regulatory function g k = g A k , as in (). Now fix a probability  < q < , and let κ be as in Conjecture . Consider a randomly chosen initial condition s(). Then, for sufficiently large n, the probability that ||s()|n  | < κ(n) will be smaller than -q  . Moreover, if n  -|s()| ≥ κ(n), then the probability that s() = (, . . . , ), i.e., that at least one of the functions takes the value g A k (s()) = , can be assumed to be < -q  ; similarly, if |s()| -n  ≥ κ(n), then the probability that s() = (, . . . , ), i.e., that at least one of the functions takes the value g A k (s()) = , can be assumed to be < -q  . This implies that with probability ≥ q, one of these two steady states will be reached after one updating step. Now consider what happens if we draw random cooperative networks and initial conditions as described at the beginning of this section for r moderately large but fixed. If |s()| = n + α √ n, then the expected value of |x| for the input vectors of the functions g A k is r( + α n ). Thus if α > , we should see on average more input vectors with norm significantly larger than r  than input vectors with norm significantly smaller than r  and vice versa. Thus if p(r, k) decays sufficiently fast to  as |k -r  | increases, |s()| -r  should be expected to have a larger absolute value and the same sign as |s()| -r  , and the effect should amplify during subsequent iterations, until a steady state or a periodic orbit with all norms close to  or n is reached. This seems to be exactly what we observe in our simulations for r ∈ {, }.

Conclusions
Very long periodic orbits are a hallmark of chaotic dynamics in Boolean networks. It has long been known that such features as few inputs per node and prevalence of significantly biased regulatory functions reduce, on average, the propensity towards chaotic behavior http://www.advancesindifferenceequations.com/content/2013/1/313 in such networks, in particular, the prevalence of very long periodic orbits. However, nontrivial upper bounds on the lengths of periodic orbits had previously been proved only for some very restricted classes of Boolean networks.
The main result of this paper, Theorem , shows the existence of such bounds if both the number of inputs and the number of outputs per node are bounded by a constant and a fixed minimal fraction of all regulatory functions is assumed to be biased. Since there are only finitely many different Boolean functions on a given set of inputs, the assumptions of this theorem imply a lower bound on the amount of bias we see in each regulatory function. However, the analogue of the theorem fails if it is only assumed that a fixed fraction of regulatory functions shows a specified minimal amount of bias, even when the number of inputs (but not outputs) of each regulatory function is bounded. This follows from a previously published result that is reproduced as Theorem (ii) here.
If both the number of inputs and that of outputs are bounded by  and all regulatory functions are biased, then  n/ is an upper bound on the length of periodic orbits in ndimensional networks. This result is proved in full generality in the preprint []; here we illustrated the main ideas of the proof by deriving the special case for cooperative networks (Theorem ). By previously published results, the upper bound is optimal. It can be further improved for systems that show p--instability for sufficiently large p <  (Theorem ). For p → this gives an upper bound of  n/ on the length of periodic orbits in n-dimensional strictly quadratic cooperative Boolean networks. Again the latter bound is optimal by a result of [], even for systems that have the stronger property of being p-unstable. Both pinstability and p--instability are formal counterparts of the general notion of sensitivity to initial conditions, which is considered to be another hallmark of chaotic dynamics. While sensitive dependence for Boolean systems can be formalized in a variety of ways (see []), p--instability for at least some p >  is implied by all meaningful formal definitions.
Cooperativity is the absence of negative interactions, in particular, negative feedback loops. Previous empirical studies had already shown that reducing the amount of negative feedback tends to reduce, on average, the propensity towards chaotic dynamics in Boolean networks []. Simulation studies reported in Section  above indicate that cooperative Boolean networks that are randomly drawn from the class of all such networks with at most r inputs per node tend to have trajectories that quickly converge to a steady state or a small periodic orbit. The effect becomes more pronounced as r increases. For r ∈ {, }, almost all sampled trajectories converged to steady states, and for r = , most of these were constant Boolean vectors. Thus the dynamics of random cooperative Boolean networks appears to become less chaotic, at least in terms of long periodic orbits and the proportion of eventually frozen nodes, as the number of inputs per node increases; exactly the opposite of what has been observed for RBNs without an assumption of cooperativity. We stated a conjecture about the distribution of sizes of sets in randomly chosen antichains that appears to conform with, but not immediately follow from, known results about these combinatorial objects. The conjecture implies that if a cooperative Boolean network and an initial condition are randomly chosen (without any restrictions on the number of inputs per node), the probability that a steady state is reached after only one step approaches  as the dimension of the network increases without bound.
These results open up several new avenues of further research. Since Theorem  assumes a fixed bound on the number of outputs, it does not directly apply to the study of RBNs from the class for which an upper bound r on the number of inputs is the only http://www.advancesindifferenceequations.com/content/2013/1/313 restriction. With probability approaching  as n → ∞, the fraction of biased regulatory functions will be very close to the fraction of biased Boolean functions among all those with r inputs, so the vast majority of such networks will be α-biased for suitable α. However, the number of outputs per node will roughly follow a Poisson distribution with parameter λ = r, with only very few nodes having a (moderately) large number of outputs. The known counterexamples do not preclude a generalization of Theorem  to this weaker assumption on the distribution of outputs, and it becomes a natural question whether one can prove or disprove such a version of the theorem. Similarly, it might be of interest to find generalizations of Theorems  and  that deal with the case r =  to corresponding networks with r > .
In Section  we defined a function p(r, k) that expresses the size distributions of sets in antichains aka Sperner families. As we argued at the end of that section, a characterization of the asymptotic behavior of the function p(r, k) should allow one to derive analytic bounds on the probability of trajectories in random cooperative Boolean networks with up to r inputs per node eventually reaching a steady state, and upper bounds on the expected lengths of the transients. In our opinion, the function p(r, k) deserves more systematic study by combinatorists. A characterization of the asymptotic behavior of this function would allow precise quantification of the well-known observation that for large r the overwhelming majority of Sperner families on {, } r contains predominantly vectors with close to r/ coordinates equal to . Of particular interest from our point of view would be a proof of Conjecture , which hypothesizes specific bounds on asymptotic behavior of this function.