Local negative circuits and cyclic attractors in Boolean networks with at most five components

We consider the following question on the relationship between the asymptotic behaviours of asynchronous dynamics of Boolean networks and their regulatory structures: does the presence of a cyclic attractor imply the existence of a local negative circuit in the regulatory graph? When the number of model components $n$ verifies $n \geq 6$, the answer is known to be negative. We show that the question can be translated into a Boolean satisfiability problem on $n \cdot 2^n$ variables. A Boolean formula expressing the absence of local negative circuits and a necessary condition for the existence of cyclic attractors is found unsatisfiable for $n \leq 5$. In other words, for Boolean networks with up to $5$ components, the presence of a cyclic attractor requires the existence of a local negative circuit.


Introduction
Boolean networks are used to model the dynamics resulting from the interactions between n regulatory components that can assume only two values, 0 and 1, and are therefore naturally described as maps from {0, 1} n to itself. Any such map uniquely identifies an asynchronous dynamics, which requires at most one component to change at each step. A regulatory graph defined by a Boolean network is a graph with one node for each regulatory component, and directed, signed edges that represent regulatory interactions. A regulation from a component to another might be observable only at certain states. Therefore, for each state of the system, a local regulatory graph is defined by considering only the regulations that can be observed at that state.
Since the explicit construction and analysis of asynchronous dynamics is generally impractical, the capability of regulatory structures to inform about the network dynamics has been often investigated. In particular, relationships have been established between the presence of circuits in regulatory graphs and the asymptotic asynchronous behaviours of Boolean networks. In absence of regulatory circuits, the dynamics always reaches a unique fixed point [13], whereas local positive circuits are required for multistationarity [5,9] and negative circuits for oscillations [5,7]. Here we consider the following question (for studies addressing related issues, see for example [5,7,8,10,12]): Question 1. Does the presence of a cyclic attractor imply the existence of a negative circuit in a local regulatory graph?
A counterexample for the multilevel case, i.e., where the discrete variables can take their values in a broader range than {0, 1}, was presented by Richard [7]. Recently, a number of counterexamples have been identified for the Boolean setting. Ruet [12] exhibited a procedure to create counterexamples in the Boolean case, for every n ≥ 7, n being the number of variables; these are maps admitting an antipodal attractive cycle and no local negative circuits in the regulatory graph. Tonello [17] and Fauré and Kaji [3] identified different Boolean versions of Richard's discrete example, that provide counterexamples to Question 1 for n = 6. A map with an antipodal attractive cycle and no local regulatory circuits also exists for n = 6 (we present such a map in appendix A).
Question 1 remains open for n ≤ 5. Even for such a small number of components, the range of possible dynamical behaviours is vast, and connections between the network regulatory structure and its associated dynamics are not immediate. However, answers to problems such as the one described in Question 1 clarify general rules and can provide guidance, for instance, to gene network modellers seeking to capture a certain dynamical behaviour.
In this work, we describe how Question 1 can be translated into a Boolean satisfiability problem (SAT). To this end, for a fixed number n of regulatory components, we consider n · 2 n Boolean variables, representing the values taken by the n components of the Boolean map on the 2 n states in {0, 1} n . We then describe how the features referred to in Question 1 can be encoded as Boolean expressions on the n · 2 n variables. More precisely, we define a Boolean formula that encodes both the absence of local negative circuits and a necessary condition for the presence of a cyclic attractor. In addition, we reduce the search space by exploiting symmetries of regulatory networks, so that, for small n, the problem can be analysed by a satisfiability solver in a few hours. The solver finds the formula unsatisfiable for n ≤ 5, and provides further examples for n = 6.
The relevant definitions and background are introduced in section 2, whereas section 3 is dedicated to recasting Question 1 as a Boolean satisfiability problem. We discuss our results in section 4.

Background
In this section, we fix some notations and introduce the main definitions. We denote by B the set {0, 1}, and consider n ∈ N. The elements of B n are also called states. The state x ∈ B n with x i = 0, i = 1, . . . , n will be denoted 0. Given x ∈ B n and a set of indices I ⊆ {1, . . . , n}, we denote byx I the state that satisfiesx I i = 1 − x i for i ∈ I, andx I i = x i for i / ∈ I. If I = {i} for some i, we simply writex i forx {i} . Given two states x, y ∈ B n , d(x, y) denotes the Hamming distance between x and y. We call n-dimensional hypercube graph the directed graph on B n with an edge from x ∈ B n to y ∈ B n whenever d(x, y) = 1.
A Boolean network is defined by a map f : B n → B n . The dynamical system defined by f is also referred to as the synchronous dynamics. The asynchronous state transition graph or asynchronous dynamics AD f defined by f is a graph on B n with an edge from x ∈ B n tox i for all i ∈ {1, . . . , n} such that f i (x) = x i . We write (x, y) for the edge (transition) from x to y.
A non-empty subset D ⊆ B n is trap domain for AD f if, for every edge (x, y), x ∈ D implies y ∈ D. The minimal trap domains with respect to the inclusion are called attractors for the dynamics of the network. Attractors that consist of a single state are called fixed points or stable states; the other attractors are referred to as cyclic attractors.
Boolean networks are used to model the interactions between regulatory components. The interactions are derived from a Boolean map f as follows. For each state x ∈ B n , we define the local regulatory graph G f (x) of f at x ∈ B n as a labelled directed graph with {1, . . . , n} as set of nodes.
The graph G f (x) contains an edge from node j to node i, also called interaction between j and i, when ; the edge is represented as j → i and is labelled with ). The label s is also called the sign of the interaction, and accounts for the regulatory effect of j upon i at the state x.
The global regulatory graph G f of f is the multi-directed labelled graph on {1, . . . , n} that contains an edge j → i of sign s if there exists a state for which the local regulatory graph contains an interaction j → i of sign s. In the global regulatory graph parallel edges are permitted to account for different regulatory effects that can be observed at different states. The sign of a path i 1 → i 2 → · · · → i k in a regulatory graph is defined as the product of the signs of its edges. A circuit in a regulatory graph is a path i 1 → i 2 → · · · → i k with i 1 = i k and such that the indices i 1 , . . . , i k−1 are all distinct. We recall a useful result which can be found in [8,Remark 1] and [11,Lemma 5.2]. Lemma 1. Let C be a circuit of G f (x) with set of vertices I. If the cardinality of {i ∈ I| f i (x) = x i } is even (resp., odd), then C is a positive (resp. negative) circuit.

Regulatory circuits and asymptotic behaviours
Following R. Thomas early conjectures [16], asymptotic properties of the asynchronous state transition graph have been connected to the existence and the signs of regulatory circuits.
Shih and Dong [13] established that, if no local regulatory circuit exists, then the map admits a unique fixed point. The result was extended to the multilevel setting by Richard [6].
The presence of multiple attractors was shown to require the existence of a local positive circuit [9]. The existence of a cyclic attractor requires instead the (global) regulatory graph to include a negative circuit. This was proved in [5] for the case of an attractive cycle (a cycle in the asynchronous dynamics that is an attractor), and in the general case of a cyclic attractor in [7].
Cyclic attractors are compatible, however, with the absence of local negative circuits. This was first shown in [7] in the multilevel case. Boolean networks with a cyclic attractor and no local negative circuits were presented in [12], with a method to create maps with antipodal attractive cycles and no local negative circuits, for n ≥ 7. Tonello [17] and Fauré and Kaji [3] exhibited maps with cyclic attractors and no local negative circuits, for n = 6. Maps with antipodal attractive cycles and no local negative circuits also exist for n = 6; a procedure that extends the one in [12] is presented, for completeness, in appendix A.
In this work we consider Question 1 in the remaining cases (n ≤ 5). We show that the problem can be approached as a Boolean satisfiability problem, and find that all maps from B n to itself with a cyclic attractor define a local negative circuit.

Automorphisms of the n-hypercube
In this section we present some relationships between Boolean networks and symmetries of the hypercube; these will be used to translate Question 1 into a Boolean expression (see section 3.3).
We first introduce some additional notations. Given I ⊆ {1, . . . , n}, ψ I denotes the map defined by ψ I (x) =x I for all x ∈ B n . We call S n the group of permutations of {1, . . . , n}; S n acts on B n by permuting the coordinates: . We consider here the maps of the form U = ψ I • σ for some I ⊆ {1, . . . , n} and some σ ∈ S n . These are all the automorphisms of the n-hypercube (see for instance [14,12]).
Given the maps U = ψ I • σ and f : The following proposition relates the asynchronous state transition graphs and regulatory graphs of f and f U , asserting that they have the same structures. In addition, albeit the signs of the interactions of the regulatory graphs can differ, the signs of the regulatory circuits are the same. An example illustrating this property is given in fig. 1.
Standard arrows j → i denote interactions with positive sign, and arrows with a vertical tip j ⊣ i represent negative interactions. The asynchronous state transition graphs have the same "shape": the map in (b) can be obtained from the map in (a) by swapping the two components, and changing 0 with 1 for the second component. In other The regulatory graphs of the two maps also have the same edges. The positive interactions on the left correspond to negative interactions on the right; however, the sign of the loop is negative in both regulatory graphs, and the sign of the circuit involving the two components is positive in both graphs.
(i) The state transition graphs AD f and AD f U are isomorphic.
(ii) For each x ∈ B n , the graphs G f (x) and G f U (U (x)), seen as unlabelled directed graphs, are isomorphic. In addition, corresponding circuits have the same signs.
the graph isomorphism is given by U . This follows from the observation that and , and therefore f U Given a circuit C in G f (x) with support on L ⊆ {1, . . . , n}, σ(L) is therefore the support of a circuit C U in G f U (U (x)). In addition, from point (i), we have that the sets {i ∈ L|f i (x) = x i } and {i ∈ σ(L)|f U i (U (x)) = U (x) i } have the same cardinality. We conclude by observing that, by Lemma 1, the circuit C is positive (resp. negative) if and only if the cardinality of {i ∈ L|f i (x) = x i } is even (resp. odd), hence if and only if the cardinality of {i ∈ σ(L)|f U i (U (x)) = U (x) i } is even (resp. odd), if and only if C U is positive (resp. negative).
It follows from the proposition that a property relating the asymptotic behaviour of the asynchronous dynamics and the regulatory circuits holds for a map if and only if it holds for any of its conjugated maps under symmetry. We will use this fact when writing Question 1 as a Boolean satisfiability problem in the next section.

Recasting Question 1 as a Boolean satisfiability problem
For each n, Question 1 requires that we determine (or exclude the existence of) a map f from B n = {0, 1} n to itself. We therefore consider as variables of the problem n · 2 n Boolean variables that we denote as We first describe how the absence of negative circuits in the local regulatory graph G f (x) can be translated into a set of expressions on the variables (2).

Imposing the absence of local negative circuits
To express the sign condition on the circuits, we consider each local graph as a complete graph on the nodes {1, . . . , n}. Then, we consider every possible circuit on this graph, and we impose that each circuit has a non-negative sign. For small values of n, this requirement leads to a satisfiability problem that is computationally manageable. The number of elementary circuits of length k in a complete graph on n nodes is given by n k (k − 1)!. Hence we have to consider, for instance, 89 circuits for n = 5, and 415 circuits for n = 6. Let C n denote the set of all possible circuits on the complete graph on {1, . . . , n}.
Given a state x ∈ B n , if an interaction exists in G f (x) from j to i, then its sign is given by the The following Boolean expression asserts that the interaction from j to i is positive: and the following Boolean expression asserts that the interaction is negative: We can now write a formula expressing that, given a state x, a circuit c is negative in G f (x), that is to say, the circuit c contains an odd number of negative interactions, the remaining interactions being positive. We write m for the length of the circuit, and c − and c + for the interactions in c with negative or positive sign, respectively. We obtain the following formula: The absence of local negative circuits in the regulatory graph is therefore specified by the formula which we can write in CNF form as

A simpler question: absence of fixed points
Before considering Question 1 in its generality, we describe how a special case of the question can be easily translated into a Boolean satisfiability problem. The question is the following: Question 2. Does the absence of fixed points imply the existence of a local negative circuit in the regulatory graph?
The absence of local negative circuits being formulated as in section 3.1, we now need to formulate the absence of fixed points. To express that a state x ∈ B n is not a fixed point for f we can write the following formula: The formula expressing the absence of fixed points for f can be written as: Since the state 0 is not fixed, there exists an index i such that f i (0) = 1. Consider a permutation σ ∈ S n that sends i to 1. The map g = σ • f • σ −1 satisfies g 1 (0) = 1; in addition, by Proposition 1, g and f have local circuits with the same signs. We can therefore assume that the first coordinate of f (0) is 1. The formula corresponding to Question 2 is therefore: The unsatisfiability of this problem is thus determined, for n = 5, in minutes, by the satisfiability solvers we considered (see section 3.4). The solvers also identify other examples of maps with no fixed points and no local negative circuits in the regulatory graph, for n = 6. The existence of a cyclic attractor is less straightforward to express; we describe our approach in the next section.

A necessary condition for the existence of a cyclic attractor
In this section we consider Question 1 in its generality. We need therefore to assert that the asynchronous state transition graph of f admits a cyclic attractor. The approach is based on the following observation. Proof. If AD f admits a cyclic attractor, then the conclusion is true for any state x in the cyclic attractor. Conversely, suppose that x is a state with the described property, and call R the set of points reachable from x in the asynchronous state transition graph. Then the minimal trap domain contained in R does not contain any fixed point, hence it must contain a cyclic attractor for AD f . Proposition 2 translates the existence of a cyclic attractor into a condition on the paths in the asynchronous state transition graph. It is, however, computationally problematic to impose that, if AD f contains a path of any length from x to y, then y is not a fixed point. We therefore consider the following condition instead.

Condition 1.
There exists a state x ∈ B n such that, for each y ∈ B n , if there is an acyclic path in AD f from x to y of length at most k, then y is not a fixed point.
It is clear from Proposition 2 that, for each k ≥ 0, Condition 1 is a necessary condition for the existence of a cyclic attractor. Our strategy is therefore to impose the absence of local negative circuits, as well as Condition 1 for increasing values of k, until we find that the problem is unsatisfiable.
In order to express Condition 1, we need to encode the existence of a given path in the asynchronous state transition graph. Given a pair of states (x, y) such that d(x, y) = 1, if x j = y j we can require that the edge (x, y) is in AD f by imposing Given a sequence of states π = (x 0 , x 1 , . . . , x k ) such that d(x i , x i+1 ) = 1, i = 0, . . . , k − 1, we can require that the sequence defines a path in AD f by imposing k constraints of the form in (8): Given a state x ∈ B n , let P k (x) denote the set of acyclic paths in the n-dimensional hypercube graph that start from x and have length less or equal to k. If π is a path in AD f , we write t(π) for the last node of the path. We express Condition 1 for a state x ∈ B n , using (5), as follows: Condition 1 requires the existence of a state x ∈ B n that verifies (10). Suppose that a map f satisfies (10) for some x ∈ B n , and that its local regulatory graphs do not admit any negative circuit. Consider j such that f j (x) = x j , and consider a permutation σ ∈ S n that swaps j and 1. Define I = {i ∈ {1, . . . , n}|σ(x) i = 0}. Then, by Proposition 1, the map f U with U = ψ I • σ satisfies (10) for x = 0, and its local regulatory graphs do not admit any negative circuit. In addition, f 1 (0) = 1. We have therefore that, to exclude the existence of maps with cyclic attractors and no local negative circuits, it is sufficient to consider expression (10) for x = 0, and assume f 1 (0) = 1. By combining (10) with (4), we have, for fixed k, the Boolean formula which we can use to answer Question 1. Notice that Q 1 is a generalisation of (7), where fewer points are required to be non-fixed. Using (10) and (9), (11) is easily written in CNF form.

Results
We created CNF files in DIMACS CNF format, a standard input format accepted by most SAT solvers. The files start with a line that begins with p cnf followed by the number of variables and the number of clauses. One line for each clause then follows. Each clause is expressed by listing the indices of the variables involved in the clause separated by spaces, using negative numbers for negated variables. A zero is added at the end of each clause line. The files were created with a Python script (source code available at github.com/etonello/regulatory-network-sat).
. The path of length 4 leaving the origin reaches a fixed point; none of the negative circuits admitted by regulatory graph are local.
Using the satisfiability solver Lingeling [1], we found that, if k is set to 2, 4, 6, 11 respectively, for n = 2, 3, 4, 5, the problem described by (11) is unsatisfiable. This means that, for n ≤ 5, all maps that admit a cyclic attractor must have a local negative circuit.
The lengths k = 2, 4, 6, 11 are the minimum lengths that lead to the unsatisfiability of the formula in (11). In other words, there exists at least one map in dimension 2 (respectively 3, 4 and 5) such that the paths of length at most 1 (respectively 3, 5 and 10) do not reach a fixed point, and the associated regulatory graph does not admit a local negative circuit. Examples of such maps are given in fig. 2, for n = 2 and n = 3. Figure 3 illustrates instead the idea of the result obtained for n = 2 and n = 3, for two special cases of asynchronous state transition graphs admitting a unique path leaving the origin: since this path reaches 3 (respectively 5) different states, the regulatory graph must admit a local negative circuit, somewhere in the state space.
The CNF file for n = 5 and k = 11 on the 160 variables consists of about 2.6 million clauses (the number of clauses for each constraint is given in table 1). The satisfiability solver Lingeling [1] was used to determine the unsatisfiability and to generate a proof, expressed in the standard DRAT notation [19]. For n = 5 and k = 11, the file for the proof is about 1GB in size, and was verified using the SAT checking tool chain GRAT [4]. The CNF file and the proof of unsatisfiability generated for n = 5, k = 11 are available as Supplementary Materials.
No local negative circuit of dimension 1 or 2 exists; however, since the only path leaving the origin has length 5, the regulatory graph must admit a local negative circuit involving all three variables. The unique attractor for the asynchronous state transition graph is a fixed point.

Conclusion
In this work we have considered the question of whether a regulatory network whose asynchronous state transition graph contains a cyclic attractor must admit a local negative circuit. For n ≥ 6, only the existence of a negative circuit in the global regulatory structure is guaranteed [7]. We have written the question as a Boolean satisfiability problem, and SAT solvers found the problem unsatisfiable for n ≤ 5. Behaviours of gene regulatory networks have been previously investigated using SAT (see, for instance [15,2,18]). Here we demonstrated that Boolean satisfiability problems can be utilised not only to examine the behaviour of a given network, but also to explore the existence of maps with desired properties, specifically, properties of the associated regulatory structure. We actually verified that, in absence of local negative circuits, Condition 1, that is implied by the existence of a cyclic attractor, cannot be satisfied, for k sufficiently large. Condition 1 requires that, for at least one state in the state space, paths of lengths at most k leaving that state cannot reach a fixed point. We found that Condition 1 with k = 2, 4, 6, 11 is sufficient for the existence of a local negative circuit in the regulatory graph, for dimensions n = 2, 3, 4, 5, respectively. The absence of local negative circuits is instead compatible with Condition 1 for k ≤ 1, 3, 5 and 10, in dimensions n = 2, 3, 4, 5, respectively.
It is natural to ask whether a relation can be established between the values identified for k via the satisfiability problems and specific properties of the n-hypercube. Such an understanding could help in clarifying the change in behaviours between n = 5 and n = 6. These points remain open for further research.