# SAT-Based Method for Finding Attractors in Asynchronous Multi-Valued Networks

Takehide Soh, Morgan Magnin, Daniel Le Berre, Mutsunori Banbara, Naoyuki

Tamura

► To cite this version:

Takehide Soh, Morgan Magnin, Daniel Le Berre, Mutsunori Banbara, Naoyuki Tamura. SAT-Based Method for Finding Attractors in Asynchronous Multi-Valued Networks. 14th International Conference on Bioinformatics Models, Methods and Algorithms (BIOINFORMATICS 2023), Feb 2023, Lisbon, Portugal. hal-03964870

HAL Id: hal-03964870

https://hal.science/hal-03964870

Submitted on 13 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# SAT-based Method for Finding Attractors in Asynchronous Multi-valued Networks

Takehide Soh<fnref n="1" /> IDª, Morgan Magnin<fnref n="2" /> IDᵇ, Daniel Le Berre<fnref n="3" /> IDᶜ, Mutsunori Banbara<fnref n="4" /> ID⁵, and
Naoyuki Tamura<fnref n="1" /> IDᵉ.

**Keywords:** Automata Network, Attractor, Constraint Programming, SAT.

**Abstract:** In this paper, we propose a SAT-based method for finding attractors of bounded size in asynchronous automata networks. The automata network is a multi-valued mathematical model which has been studied for the qualitative modeling of biological regulatory networks. An attractor is a minimal set of states in automata networks that cannot be escaped and thus loops indefinitely. Attractors are crucial to validate the initial design of a biological model and predict possible asymptotic behaviors, e.g., how cells may result through maturation in differentiated cell types. Developing an efficient computational method to find attractors is thus an important research topic. Our contribution is a translation of the problem of finding attractors of automata networks into a sequence of propositional satisfiability (SAT) problems. We also propose to add two optional constraints to improve the computation time of attractors. Experiments are carried out using 30 automata networks, 8 coming from real biological case studies and 22 crafted ones with controlled attractor size. The experimental results show that our method scales better than the state-of-the-art ASP method when the size of the attractors increases.

## 1 INTRODUCTION

**Background.** Understanding the mechanisms involved in biological regulation is a fundamental issue in analyzing living systems. The formal study of the dynamics of biological systems raises many problems, e.g., identification of attractors, bifurcations, and reachability, that are combinatorial by essence. Making these issues scalable requires designing efficient methods that rely on efficient programming frameworks. Considering this viewpoint, the long-term behavior of a regulatory network's dynamics is of specific interest (Wuensche, 1998). Such outline has been interpreted as distinct responses of the organism, such as differentiating into distinct cell

types in multicellular organisms (Huang et al., 2005). Moreover, when refining a model of a living system, one way to remove inconsistencies or to predict missing information in biological models consists in comparing the attractors of the model with the experimentally observed long-term behavior. This explains why our property of interest in this paper is the computation of attractors, an attractor being a minimal set of states that cannot be escaped and thus loops indefinitely.

**Automata Network.** Various kinds of mathematical models have been proposed for the qualitative modeling of biological regulatory networks (BRNs). These models include neural networks, differential equations, Petri Nets, Boolean Networks (BN), probabilistic Boolean networks, and other multi-valued models such as synchronous/asynchronous Automata Networks (AN).

In this paper, we focus on a subclass of automata networks called Asynchronous Automata Networks (AAN) (Folschette et al., 2015; Paulevé, 2016a),

<footnote n="1">¹Kobe University, Information Infrastructure and Digital Transformation Initiatives Headquarters, 1-1, Rokko-dai, Nada, Kobe, Hyogo 657-8501 Japan</footnote>

<footnote n="2">²Nantes Université, École Centrale Nantes, CNRS, LS2N, UMR 6004, F-44000 Nantes, France</footnote>

<footnote n="3">³Univ. Artois, CNRS, Centre de Recherche en Informatique de Lens (CRIL), F-62300 Lens, France</footnote>

<footnote n="4">⁴Nagoya University, Graduate School of Informatics, Furo-cho, Chikusa-ku, Nagoya, 464-8601, Japan
soh@lion.kobe-u.ac.jp, morgan.magnin@ec-nantes.fr, leberre@cril-lab.fr, banbara@nagoya-u.jp, tamura@kobe-u.ac.jp</footnote>

<footnote>a ID: https://orcid.org/0000-0001-5897-9192</footnote>

<footnote>b ID: https://orcid.org/0000-0001-5443-0506</footnote>

<footnote>c ID: https://orcid.org/0000-0003-3221-9923</footnote>

<footnote>d ID: https://orcid.org/0000-0002-5388-727X</footnote>

<footnote>e ID: https://orcid.org/0000-0002-5466-1010</footnote>
which is convenient to model BRNs. In Automata Networks, biological components (e.g., genes) are represented as abstracted in the form of *automata*. The different local states (that are not restricted to Boolean values) of each automaton correspond to the different discrete qualitative levels of the components represented by the automaton. Interaction between biological components is modeled in an atomic way by local transitions on the automaton, where each one is conditioned by a set of required local states in different automata and can modify the local state of a unique automaton. In other words, AANs allow to have multiple requirements for a local transition to occur, but *not* to synchronize several local transitions in different automata. In this sense, they are considered *Asynchronous Automata Networks*. Synchronous semantics could also be defined on such a model, but this is out of the scope of the current paper. While building a model, the choice of an appropriate semantics is of crucial importance with regard to dynamical properties. This discussion is a matter of research in itself, with many papers balancing the merits of each semantics. For some applications, like the biological ones, asynchronous semantics is said to capture more realistic behaviors in the sense that, at a given time, a single gene can change its expression level. For a more comprehensive discussion on this aspect, (Garg et al., 2008) considers the differences and respective merits of synchronous and asynchronous semantics to model biological networks and identify attractors. Note also that, depending on the chosen semantics (asynchronous or synchronous), AANs encompass the well-known Boolean frameworks of René Thomas (Thomas, 1973) and Stuart Kauffman (Kauffman, 1969).

**Studies on Automata Networks.** Automata networks are a very concise and convenient framework to model BRNs. Their connection with other formal models used in the Systems Biology community (like logical networks, which they generalize) has been formally proven in the past (Paulevé et al., 2012; Folschette et al., 2015). In a more practical perspective, respective translations between model frameworks have previously been implemented and integrated altogether in the CoLoMoTo initiative (Naldi et al., 2018) through the tool implementing (Paulevé, 2017). As such, the existing literature defining logical models of gene regulatory networks or signaling networks, e.g., (Grieco et al., 2013; Abou-Jaoudé et al., 2015; Sahin et al., 2009) (among many others) can benefit from results obtained on automata networks (Levy et al., 2018). Many approaches have thus been developed to analyze dynamical properties on automata networks (Paulevé et al., 2013; Paulevé, 2016b; Fitime et al., 2017; Chai et al., 2020).

In (Rougny et al., 2021), the authors offer a representative example of the benefits of automata networks to capture the behavior of complex biological systems like the circadian clock while recognizing that such works would greatly benefit the analysis of long-term behaviors, i.e. attractors, untractable when too many nodes and interactions are involved.

**Computing Attractors.** In recent years, the design of efficient approaches to identify attractors in biological regulator networks mainly relied on the ASP framework (Gebser et al., 2012). In (Ben Abdallah et al., 2017), the authors proposed a first ASP encoding to enumerate attractors up to a given length. Their approach is iterative and without any assumption on the topology of the interaction graph. Contrastingly, the authors of (Khaled and Benhamou, 2020) enumerate attractors in asynchronous Boolean networks that have the form of a circuit using ASP by introducing a new resolution semantics that does not use the usual negation by failure. Instead, they rely on a weak version of the negation by failure that allows them to provide an enumeration approach preventing the simulation of the underlying networks. This approach targets networks with a given structure (all nodes must have an incoming edge and some results have been obtained only for cyclic interaction graphs). In this paper, we do not assume any specific structure for the networks.

**Using SAT instead of ASP.** If ASP solvers became very popular thanks to their versatile modeling capability, SAT solvers became in two decades de facto NP-complete oracles for solving a wide range of problems, either decision or optimization ones, on a wide range of complexity (from NP-complete to PSPACE-complete) (Biere et al., 2021) thanks to the availability of numerous, robust and increasingly more efficient solvers. The use of SAT solvers requires translating the original problem into a set of clauses, with Boolean variables. Such translation is most often ad-hoc, but can also use an intermediate modeling language for which a generic translation to SAT exists. In this paper, we use constraints with integer variables as an intermediate modeling language. We propose a constraints-based method for finding attractors of bounded size in asynchronous automata networks implemented by translation into a sequence of SAT problems. Our contribution is the translation of the problem of finding attractors of automata networks into the enumeration of the solutions of a set of constraints. Notably, we redefine attractors by using two predicates about local playability and global playability. Moreover, we propose to improve the computation time of our approach thanks
<figure><img src="image_2.png" /><figcaption>Figure 1: Automata Network Example: grey is used to represent one possible state of the model, that is (a<sub>1</sub>, b<sub>0</sub>, c<sub>2</sub>).</figcaption></figure>

to the addition of two optional constraints. The non-
Boolean constraints are finally translated to a Boolean
domain to complete the translation into a SAT prob-
lem, which allows using a wide variety of very effi-
cient SAT solvers. Experiments are carried out using
30 automata networks. The SAT approach is com-
pared favorably to the existing ASP approach on in-
stances containing larger-sized attractors.

2 PRELIMINARIES

**2.1 Automata Networks**

**Definition 1.** An automata network is a triple
(Σ, S, T) where Σ, S and T are defined as follows:

• Σ is a set of automata. We use a symbol α to denote an automaton α ∈ Σ.

• S = Π_{α∈Σ} S_α is the finite set of **global states**,
where S_α = {α_0, α_1, ...} is a finite set of local
states of an automaton α ∈ Σ and α_ν denotes an
automaton α is in the state v.

* T is a finite set of transitions. A transition can change the state of only one automaton α and is represented in the form of (α<sub>u</sub>, α<sub>v</sub>, C). It changes the state of the automaton α from u to v using a condition C ⊂ ⋃<sub>α∈Σ</sub> S<sub>α</sub>. In addition, we define the transitions changing the state of an automaton α as T<sub>α</sub> ⊆ {(α<sub>u</sub>, α<sub>v</sub>, C) | α<sub>u</sub>, α<sub>v</sub> ∈ S<sub>α</sub>, u ≠ v, C ⊂ ⋃<sub>α∈Σ</sub> S<sub>α</sub>}.

Figure 1 shows an example of an automata net-
work. It contains three automata named a, b, c. We
focus on the automaton a to explain the example.
The automaton a can be in one of the two states 0 or 1. The automaton a has three transitions T_a = { (a_0, a_1, {c_0}) , (a_0, a_1, {b_1, c_2}) , (a_1, a_0, {b_0}) } . The transition (a_0, a_1, {b_1, c_2}) means that the transition changes the state of the automaton a from 0 to 1 when

the two conditions $b_1$ – automaton b is in state 1 – and
$c_2$ – automaton c is in state 2 – hold.

## 2.2 Playability in Asynchronous Automata Network

This section explains how a transition $t \in T$ on a sin-
gle automaton is used to define the transitions among
global states $g \in S$ in an automata network $(\Sigma, S, T)$.

**Asynchronous Update.** The state of an automaton can be updated either synchronously or asynchronously. The synchronous update is the simplest method and all automata are updated simultaneously. The asynchronous update is a more realistic model and all automata are not necessarily updated. A difference between the two update methods appears in their state transition graphs (STGs). The synchronous update is deterministic and thus the outdegree of each node in STG is at most one. The asynchronous update is nondeterministic and thus the outdegree of each node in STG is not limited. In this paper, we focus on the asynchronous update method. In particular, we allow only one automaton to update, called *generalized asynchronous*, which is known to contain many other update schemes as its special case (Gershenson, 2002).

In the following, we treat a global state g as a set
rather than a tuple if it is clear from the context.

**Definition 2.** (Local Playability of Transition). A transition $(\alpha_u, \alpha_v, C) \in T$ is *locally playable* in a global state $g \in S$ if $\{\alpha_u\} \cup C \subseteq g$ holds.

Let $g_i, g_j$ be two global states in $S$ such that $g_i \neq g_j$ and $t = (\alpha_u, \alpha_v, C) \in T_\alpha$ is a transition changing the state of $\alpha \in \Sigma$.

**Definition 3.** (Global Playability of Transition). A transition $t = (\alpha_u, \alpha_v, C) \in T$ is *globally playable* from a global state $g_i \in S$ to a global state $g_j \in S$ such that the transition $t$ is *locally playable* in the global state $g_i, \alpha_v \in g_j$, and $|g_i \setminus g_j| = 1$.
<figure><img src="image_3.png" /><figcaption>Figure 2: State Transition Graph Drawn from the Example Automata Network: red, blue, and grey are used to represent attractors of size 1, 2 and 4, respectively.</figcaption></figure>

Note that the condition $|g_i \setminus g_j| = 1$ guarantees the difference between $g_i$ and $g_j$ is only one automaton, i.e., the automaton $\alpha$ whose state is changed from $u$ to $v$ according to the generalized asynchronous update.

## 2.3 Attractors

We define two predicates about the playability of a transition $t$ with respect to global states which facilitate the understanding of the definition of trap domains and our proposed constraint model described in Section 3.2.

$$P(t, g_i) \stackrel{def}{\iff} \{\alpha_u\} \cup C \subseteq g_i$$

$$P(t, g_i, g_j) \stackrel{def}{\iff} P(t, g_i) \wedge \alpha_v \in g_j \wedge |g_i \setminus g_j| = 1$$

Then, trap domains and attractors are defined as follows.

**Definition 4.** (Trap Domain). $G \subseteq S$ is a trap domain if and only if $\forall t \in T. \forall g_i \in G. (P(t, g_i) \rightarrow \exists g_j \in G. (P(t, g_i, g_j))$

**Definition 5.** (Attractor (or Minimal Trap Domain)). $G \subseteq S$ is an attractor if and only if there is no $G' \subset G$ such that $G'$ is a trap domain.

That is, attractors are trap domains that are inclusion-minimal. The size of an attractor $G$ is the cardinality of $G$ denoted as $|G|$.

Figure 2 shows the state transition graph of the automata network of Figure 1. All attractors in this network are:

* Size 1: $\{(a_0, b_1, c_1)\}, \{(a_1, b_1, c_0)\}$

* Size 2: $\{(a_1, b_1, c_2), (a_1, b_1, c_1)\}$

* Size 4: $\{(a_0, b_0, c_0), (a_0, b_0, c_1), (a_0, b_0, c_2), (a_1, b_0, c_0)\}$

In the state transition graph, attractors can be seen as strongly connected components (SCCs).

**Definition 6.** Bounded Attractor Enumeration Problem (BAE). BAE is the problem of finding all attractors sized less than or equal to a given size.

* **Input:** An automata network $(\Sigma, S, T)$ and a bound $k$.

* **Output:** All attractors of the given automata network whose sizes are less than or equal to $k$.

Given a BAE whose inputs are the automata network of Figure 1 and $k=1$, the output is 2 attractors. In the case of $k=2$, the output is 3 attractors. In the case of $k=4$, the output is 4 attractors.

# 3 A SAT-BASED MODEL

The proposed SAT-based approach is based on two parts. The first one is our constraint model for constructing trap domains. The second one is the incremental computation of trap domains. We ensure the minimality of trap domains by computing them starting from the smallest one.

## 3.1 Propositional logic recap

We consider propositional formulas in conjunctive normal form (CNF). A literal is either a propositional variable or its negation. A clause is a disjunction of literals and a CNF formula is a conjunction of clauses. An assignment is a mapping from a set of propositional variables to the Boolean values true or false. A propositional formula is satisfiable (SAT) if there
**Algorithm 1:** Algorithm for Computing Attractors

**Data:** BAE: Automata network $(\Sigma, S, T)$ and a bound $k$

**Result:** All attractors of size $\le k$

1 **Function** Main((Σ, S, T) and a bound k):

2 $\Omega := \{\};$ // set of trap domains

3 **for** i = 1 to k **do**

4 $\Psi_i := \text{Construct } \Psi_i^L \wedge \Psi_i^G \wedge \Psi_i^T \wedge \bigwedge_{\omega \in \Omega} \text{BlockClauses}(\omega, i);$

5 **while** $\Psi_i$ has a solution $\omega$ **do**

6 $\Omega := \Omega \cup \{f(\omega)\};$ // $f(\omega)$ returns a trap domain from the solution $\omega$.

7 $\Psi_i = \Psi_i \wedge \text{BlockClauses}(\omega, i);$

8 **return** $\Omega$;

9 **Function** BlockClauses(Solution $\omega$, Current Bound i):

10 generate clauses $\Psi_i^\omega$ avoiding computation of trap domains containing a trap domain $f(\omega);$

11 **return** $\Psi_i^\omega$

exists an assignment that satisfies it and is unsatisfiable (UNSAT) otherwise. In addition to $\lor, \land, \neg$, we use Boolean operators $\leftrightarrow, \rightarrow$ which are easily translated in CNF formulas. A SAT solver is a program that computes a solution (if any) from a CNF formula.

Note that the constraint model in the following section is partially not propositional, i.e., we use some constraints over integer variables for clarity to the reader. Those non-Boolean constraints are surrounded by symbols $\llcorner$. They will be translated into Boolean constraints using the encodings in Table 1 in Section 5.

## 3.2 Constraint Model for Trap Domain

**Variables.** For each $t \in T$ and $1 \le i \le k$, we introduce a propositional variable $p_{t,i}$ that is true when a transition $t \in T$ is locally playable from a global state $g_i \in S$. For each $t \in T$ and $1 \le i \ne j \le k$, we introduce a propositional variable $p_{t,i,j}$ that is true when a transition $t \in T$ is globally playable from a global state $g_i \in S$ to a global state $g_j \in S$. For each $\alpha \in \Sigma$ and $1 \le i \le k$, we introduce an integer variable $x_{\alpha,i}$ that is $x_{\alpha,i} = u$ iff an automaton $\alpha$ is in a local state $u \in S_\alpha$ at a global state $g_i \in S$, i.e., $\alpha_u \in g_i$.

**Constraints.** We introduce the following constraints to compute trap domains using the variables explained. A symbol $\bullet t$ denotes the condition of the local playability of a transition $t$, that is, $\bullet t = \{\alpha_u\} \cup C$ when $t = (\alpha_u, \alpha_v, C)$. The first constraint is the definition of the local playability of transitions. We denote this constraint as $\Psi_k^L$.

$$p_{t,i} \leftrightarrow \left( \bigwedge_{\alpha_u \in \bullet t} \llcorner x_{\alpha,i} = u \rrcorner \right) \quad (t \in T, 1 \le i \le k)$$

The second constraint is the definition of the global playability of transitions. We denote this constraint as $\Psi_k^G$.

$$p_{t,i,j} \leftrightarrow \left( \begin{array}{l} p_{t,i} \wedge \llcorner x_{\alpha,j} = v \rrcorner \wedge \bigwedge_{\alpha' \in \Sigma \setminus \alpha} \llcorner x_{\alpha',i} = x_{\alpha',j} \rrcorner \\ (t = (\alpha_u, \alpha_v, C) \in T, 1 \le i \ne j \le k) \end{array} \right)$$

The third constraint is the definition of the trap domain. We denote this constraint as $\Psi_k^T$.

$$p_{t,i} \rightarrow \left( \bigvee_{\substack{j \ne i \\ 1 \le j \le k}} p_{t,i,j} \right) \quad (t \in T, 1 \le i \le k)$$

**Solutions.** Let $\omega$ be any solution of $\Psi_k^L \wedge \Psi_k^G \wedge \Psi_k^T$, which are assignments to propositional or integer variables $t_i, t_{i,j}$, and $x_{\alpha,i}$ that satisfies it. In a solution $\omega$, $x_{\alpha,i} = u$ means $\alpha_u \in g_i$. Let $f$ be a mapping from a solution to global states. For instance, suppose that a solution $\omega$ contains $x_{a,0} = 1, x_{b,0} = 1, x_{c,0} = 2$ and $x_{a,1} = 1, x_{b,1} = 1, x_{c,1} = 1$. Then, global states $f(\omega)$ are $\{g_1, g_2\}$, where $g_1 = (a_1, b_1, c_2)$ and $g_2 = (a_1, b_1, c_1)$.

**Proposition 1.** Global states $f(\omega)$ from solutions $\omega$ of $\Psi_k^L \wedge \Psi_k^G \wedge \Psi_k^T$ are trap domains.

(Proof) The two propositional variables $p_{t,i}$ in $\Psi_k^L$ and $p_{t,i,j}$ in $\Psi_k^G$ correspond to the definition of the two predicates $P(t, g_i)$ and $P(t, g_i, g_j)$, respectively. The constraint $\Psi_k^T$ is the same as the definition of the trap domain: $\forall t \in T. \forall g_i \in G. (P(t, g_i) \rightarrow \exists g_j (P(t, g_i, g_j)).$ Thus, global states $f(\omega)$ are trap domains. $\square$

Note that we do not need $\leftarrow$ of $\Psi_k^G$ to solve BAE since $p_{t,i,j}$ appears only positively in $\Psi_k^T$.
## 3.3 Computing Attractors

Algorithm 1 shows an algorithm to compute attractors using the constraint model from the previous section. Input is a BAE instance, i.e., an automata network $(\Sigma, S, T)$ and a bound $k$. Output is all attractors whose sizes are less than or equal to $k$. Line 2 initializes a set $\Omega$ that contains attractors. Lines 3 to 7 incrementally compute attractors from size 1 to size $k$. Line 5 launches a SAT solver to compute an attractor of size $i$. Lines 9 to 11 discard the attractors found so far. Specifically, the blocking clauses $\Psi_i^\omega$ are as follows.

$$ \bigwedge_{1 \le j \le i} \bigvee_{\alpha \in \Sigma} \llbracket x_{\alpha, j} \ne u \rrbracket \quad (g \in f(w)) $$

**Proposition 2.** Algorithm 1 outputs all attractors sized less than or equal to $k$.

(Proof) We use mathematical induction on $k$. i) in the case of $k = 1$, Algorithm 1 outputs all trap domains sized 1. Obviously, those trap domains are minimal. Thus, the proposition holds at $k = 1$. ii) if the proposition holds at $k$, then Algorithm 1 outputs all trap domains sized $k+1$ and those trap domains are minimal since all trap domains less than or equal to $k$ are enumerated and blocked. Therefore, by i) and ii), the proposition holds. $\square$

# 4 Tuning Our Model

To represent the next constraints, we need to introduce new propositional variables. A propositional variable $p_{t,i} \to$ is true iff a propositional variable $p_{t,i,j}$ is true for some $j$. A propositional variable $p_{t,j} \leftarrow$ is true iff a propositional variable $p_{t,i,j}$ is true for some $i$. A propositional variable $p_t$ is true iff propositional variables $p_{t,i} \to$ or $p_{t,i} \leftarrow$ are true for some $i$.

$$ p_{t,i} \leftrightarrow \bigvee_{1 \le j \le k} p_{t,i,j} \quad (t \in T, 1 \le i \le k) $$

$$ p_{t,j} \leftrightarrow \bigvee_{1 \le i \le k} p_{t,i,j} \quad (t \in T, 1 \le j \le k) $$

$$ p_t \leftrightarrow \bigvee_{1 \le i \le k} p_{t,i} \leftrightarrow p_{t,i} \quad (t \in T) $$

Note that the propositional variable $p_{t,i} \to$ is redundant because it is the same as the propositional variable $p_{t,i}$ in Section3.2. We rewrite here for clarity to the reader.

## 4.1 For attractors of size greater than 1

When we search attractors sized $\ge 2$, we can add the following constraints since each global state in attractors should have at least one indegree and one outdegree in STGs.

$$ \bigvee_{t \in T} p_{t,i} \quad (1 \le i \le k) $$

$$ \bigvee_{t \in T} p_{t,j} \quad (1 \le j \le k) $$

These constraints are useful when there are a huge number of attractors sized 1 since we do not need to block each of them individually.

## 4.2 With optional cycle constraints

Let $T'_\alpha$ be a set of transitions, which changes the state of an automaton $\alpha$, i.e., $T'_\alpha \subseteq T_\alpha$. We call $T'_\alpha$ are cycles of transitions in an automaton $\alpha$ when $T'_\alpha$ are (not simple, and possibly multiple) cycles.

An important property of attractors sized $\ge 2$ is that these consist of (not simple, and possibly multiple) cycles of transitions in each automata $T_\alpha$. Consider the attractor $\{(a_1, b_1, c_2), (a_1, b_1, c_1)\}$ in the example, it consists of the cycle of transitions of the automata $c$, i.e., $(c_1, c_2, \{a_1, b_1\})$ and $(c_2, c_1, \{\})$. In general, we can say the following proposition.

**Proposition 3.** (Transitions consisting attractors). Let $T'_\alpha$ be a set of transitions such that $T'_\alpha \subseteq T_\alpha$. Given an attractor $G$, we denote all transitions globally playable in $G$ as $\cup_{\alpha \in \Sigma} T'_\alpha$. Then, for each $\alpha$, $T'_\alpha$ are cycles of transitions or an empty set.

(Proof) We use a proof by contradiction. Suppose that the set of transitions $T'_\alpha$ are not cycles of transitions. Then, there is at least one state that cannot be reached by using transitions of $T'_\alpha$. Since any attractors are strongly connected components in its state transition graph, it is a contradiction. $\square$

For instance, simple cycles of transitions in Figure 1 are listed as follows.

$$ \gamma_1 = \{(a_0, a_1, \{c_0\}), (a_1, a_0, \{b_0\})\} $$

$$ \gamma_2 = \{(a_0, a_1, \{b_1, c_2\}), (a_1, a_0, \{b_0\})\} $$

$$ \gamma_3 = \{(c_0, c_1, \{a_0\}), (c_1, c_0, \{a_0, b_0\})\} $$

$$ \gamma_4 = \{(c_1, c_2, \{a_1, b_1\}), (c_2, c_1, \{\})\} $$

$$ \gamma_5 = \{(c_1, c_2, \{a_0, b_0\}), (c_2, c_1, \{\})\} $$

Then, all sets of transitions globally playable in attractors are combinations from the following list.

$$ T'_a \in \{\}, \gamma_1, \gamma_2, \gamma_1 \cup \gamma_2 \} $$

$$ T'_b \in \{\} $$

$$ T'_c \in \{\}, \gamma_3, \gamma_4, \gamma_5, \gamma_3 \cup \gamma_4, \gamma_3 \cup \gamma_5, \gamma_4 \cup \gamma_5, \gamma_3 \cup \gamma_4 \cup \gamma_5 \} $$

For instance, all transitions globally playable in the attractors sized $\ge 2$ in the example are as follows.

* $\{(a_1, b_1, c_2), (a_1, b_1, c_1)\}$
  - $T'_a = \{\}$ and $T'_c = \gamma_4$.

* $\{(a_0, b_0, c_0), (a_0, b_0, c_1), (a_0, b_0, c_2), (a_1, b_0, c_0)\}$
$$-T'_{a} = \gamma_{1} \text{ and } T'_{c} = \gamma_{3} \cup \gamma_{4}.$$

Using Proposition 3, we can add the following optional constraints. Let $\Gamma$ be the set of simple cycles in a given automata network. We introduce a propositional variable $p_{\gamma}$ that is true if a cycle $\gamma \in \Gamma$ is contained in a solution. We add the following constraints.

$$p_{\gamma} \rightarrow \bigwedge_{t \in \gamma} p_t \quad (\gamma \in \Gamma)$$

$$p_t \rightarrow \bigvee_{\substack{\gamma \in \Gamma \\ t \in \gamma}} p_{\gamma} \quad (t \in T)$$

The first constraint ensures that if a cycle $\gamma$ is contained in a solution then all its transitions also belong to the solution. The second constraint ensures that if a transition $t$ is playable then one of the cycles $\gamma$ that contain the transition $t$ belongs to the solution.

## 4.3 With optional ordering constraints

Solutions of the basic constraints discussed in Section 3.2 contain many symmetries. In the case of an attractor $\{(a_1, b_1, c_2), (a_1, b_1, c_1)\}$ from the example in Figure 2, two solutions can be computed. The first one $\omega_1$ is:

$$\bullet \quad x_{a,0} = 1, x_{b,0} = 1, x_{c,0} = 2 \text{ (i.e., } (a_1, b_1, c_2))$$

$$\bullet \quad x_{a,1} = 1, x_{b,1} = 1, x_{c,1} = 1 \text{ (i.e., } (a_1, b_1, c_1))$$

The second one $\omega_2$ is:

$$\bullet \quad x_{a,0} = 1, x_{b,0} = 1, x_{c,0} = 1 \text{ (i.e., } (a_1, b_1, c_1))$$

$$\bullet \quad x_{a,1} = 1, x_{b,1} = 1, x_{c,1} = 2 \text{ (i.e., } (a_1, b_1, c_2))$$

Those two solutions are identical in the sense that both solutions represent the same attractor. The transitions are simply not in the same order.

Using the current constraints, the solver could find any of them. The other one will be later discarded by the blocking clauses.

It is possible to force the order of the transitions, i.e. to represent each attractor only once. This is a standard technique in constraints programming called "symmetry breaking" (Crawford et al., 1996).

Given a sequence of global states $g_1, g_2, \dots, g_k$, we add constraints forcing $g_1 \prec g_2 \wedge \dots \wedge g_{k-1} \prec g_k$, where $g_i \prec g_{i+1}$ means the global state $g_i$ is lexically strictly smaller than $g_{i+1}$. Let $\Sigma$ be $\{a^1, a^2, \dots, a^n\}$. Then the constraints $g_i \prec g_{i+1}$ are defined as follows. Note that $\prec$ is $<$ if $j = n$ and $\le$ otherwise.

$$[[x_{a^j, i} \le x_{a^j, i+1}]] \quad (j = 1)$$

$$\left( \bigwedge_{1 \le m \le j-1} [[x_{a^m, i} = x_{a^m, i+1}]] \right) \rightarrow [[x_{a^j, i} \prec x_{a^j, i+1}]] \quad (2 \le j \le n)$$

<table><caption>Table 1: Encoding Constraints</caption><thead><tr><th>Constraints</th><th>Clauses</th></tr></thead><tbody><tr><td>[[x = v]]</td><td>$P_{x=v}$</td></tr><tr><td>[[x = x<sup>'</sup>]]</td><td>$\bigwedge_{\substack{d \in D(x) \\ d' \in D(x') \\ d=d'}} P_{x=d} \leftrightarrow P_{x'=d'}$</td></tr><tr><td>[[x < x<sup>'</sup>]]</td><td>$\bigwedge_{\substack{d \in D(x) \\ d' \in D(x') \\ \neg(d < d')}} \neg P_{x=d} \vee \neg P_{x'=d'}$</td></tr><tr><td>[[x \le x<sup>'</sup>]]</td><td>$\bigwedge_{\substack{d \in D(x) \\ d' \in D(x') \\ \neg(d \le d')}} \neg P_{x=d} \vee \neg P_{x'=d'}$</td></tr></tbody></table>

By this constraint, only $\omega_2$ can be a solution. It is known that symmetry breaking constraints are often effective for unsatisfiable instances since it reduces search spaces.

# 5 Encoding Constraints in CNF

The constraints in Section 3.2 and 4 are almost propositional constraints but we need to encode some integer constraints surrounded by $[[}$]. There have been many encoding methods that encode integer constraints into propositional clauses: direct encoding (de Kleer, 1989; Walsh, 2000), log encoding (Iwama and Miyazaki, 1994; Gelder, 2008), order encoding (Crawford and Baker, 1994; Tamura et al., 2009), and hybrid encoding (Soh et al., 2017). The benchmarks used in the next section contain automata with at most 3 states, and biological systems rarely have more than 4 states (Folschette et al., 2015). Thus, we use the direct encoding because it is simple and efficient when there are few integer or domain values to encode.

Let $D$ be a mapping from a variable to its domain values. Let $p_{x=d}$ be a propositional variable that is true when $x = d$ holds, where $d \in D(x)$. We encode integer variables' axiom, a variable is assigned exactly one value, as follows.

$$\bigvee_{d \in D(x)} p_{x=d} \quad (x \in X)$$

$$\bigwedge_{\substack{d, d' \in D(x) \\ d \neq d'}} \neg p_{x=d} \vee \neg p_{x'=d'} \quad (x \in X)$$

Using the propositional variable $p_{x=d}$, integer constraints are encoded as summarized in Table 1. Note that the translation of symmetry breaking constraints in Section 4.3 can be done in linear size by using Tseitin translation (Tseitin, 1968).
<table><caption>Table 2: Number of Constraints and Clauses</caption><thead><tr><th>Name</th><th>#Constraints</th><th>#Clauses</th></tr></thead><tbody><tr><td>Ψ<sub>k</sub><sup>L</sup></td><td>O(|T||k|)</td><td>O(|T||k||●t|)</td></tr><tr><td>Ψ<sub>k</sub><sup>G</sup></td><td>O(|T||k|<sup>2</sup>||Σ|)</td><td>O(|T||k|<sup>2</sup>||Σ||S<sub>α</sub>|<sup>2</sup>)</td></tr><tr><td>Ψ<sub>k</sub><sup>T</sup></td><td>O(|T||k|)</td><td>O(|T||k|)</td></tr><tr><td>Cycle</td><td>O(|T|)</td><td>O(|T|)</td></tr><tr><td>Sym. Break.</td><td>O(|k||Σ|)</td><td>O(|k||Σ||S<sub>α</sub>|<sup>2</sup>)</td></tr></tbody></table>

Table 2 shows the order of the number of constraints and clauses. We here assume the followings to ease the discussion: $|S_α|$ and $|\bullet t|$ are constant values; $|T| > |\Gamma|; |\Gamma| \gg |\gamma|$. The table shows us that the encoding is dominated by the definition of globally playable transitions when $k$ becomes large.

# 6 Experiments

To check the performance of the proposed SAT-based method, we carried out an experimental comparison against the state-of-the-art ASP-based method from (Ben Abdallah et al., 2017). Our approach is implemented using Scala 2.12.4 to generate the CNF from the automata network description. We use two SAT solvers in our implementation. For $k = 1$, we use BDD_MINISAT_ALL (Toda and Soh, 2016). This solver is good at problems that are easy but have many solutions. For $k > 1$, we use CaDiCaL (Biere et al., 2020). This solver is good at problems that are hard to solve.

## 6.1 Experimental Conditions

**Existing Benchmark.** 30 automata networks are used. The first 8 networks are all instances used in the previous study (Ben Abdallah et al., 2017), which are inspired from real organisms and found in the literature: Example (Ben Abdallah et al., 2017), Lambda phage (Thieffry and Thomas, 1995), Trp- reg (Simão et al., 2005), Fission-yeast (Davidich and Bornholdt, 2008), Mamm (Fauré et al., 2006), Tcr- sig (Klamt et al., 2006), FGF (Mbodj et al., 2013), and T-helper (Abou-Jaoudé et al., 2015). These net- works contain relatively small-sized attractors.

**Crafted Benchmark.** We further contribute to an artificial benchmark set named star to check the scalability w.r.t. the number of automata, the number of transitions, and the size of cyclic attractors. This benchmark is available on the repository<fnref n="1" />.

Given a parameter $N > 1$, starN is an automata network ($\Sigma, S, T$) defined as follows.

$$ \bullet \Sigma = \{a^i | 1 \le i \le N\} $$

$$ \bullet S = \{a_0^1, a_1^1\} \times \{a_0^2, a_1^2\} \times \cdots \times \{a_0^N, a_1^N\} $$

<figure><img src="image_4.png" /><figcaption>Figure 3: The State Transition Graph of star03</figcaption></figure>

$$ \bullet T = \{(a_0^i, a_1^i, \{a_0^j | 1 \le j \le N, i \ne j\}) | 1 \le i \le N\} \cup \{(a_1^i, a_0^i, \{\}) | 1 \le i \le N\}. $$

Each starN automata network contains $N$ automata that have two states, $2N$ transitions, and only one attractor of size $N+1$.

star03 is given below.

$$ \bullet \Sigma = \{a^1, a^2, a^3\} $$

$$ \bullet S = \{a_0^1, a_1^1\} \times \{a_0^2, a_1^2\} \times \{a_0^3, a_1^3\} $$

$$ \bullet T = T_{a1} \cup T_{a2} \cup T_{a3} $$

$$ \bullet T_{a1} = \{(a_0^1, a_1^1, \{a_0^2, a_0^3\}), (a_1^1, a_0^1, \{\})\} $$

$$ \bullet T_{a2} = \{(a_0^2, a_1^2, \{a_0^1, a_0^3\}), (a_1^2, a_0^2, \{\})\} $$

$$ \bullet T_{a3} = \{(a_0^3, a_1^3, \{a_0^1, a_0^2\}), (a_1^3, a_0^3, \{\})\} $$

Figure 3 shows the attractor of star03 in gray.

**Bound k.** For all benchmarks, the attractors are known. In this experiment, we set a bound $k$ for each automata network so that solvers can compute the largest attractors.

**Environment.** We execute all experiments on the machine that equips 3GHz CPU and 16GB RAM. The time limit is 3 hours.

## 6.2 Results

Table 3 shows comparisons between the proposed method and the state-of-the-art existing ASP-based method (Ben Abdallah et al., 2017). The first column denotes the name of the automata networks. The second column denotes the number of automata included. The third column denotes the number of transitions included. The fourth column denotes the sizes of the attractors. The number of attractors for each size is displayed between the parenthesis. The fifth column denotes the bound $k$ given to the solvers. The sixth column denotes the maximum number of states of automata denoted as $S_{max}$. The seventh to tenth columns denote the CPU time of the proposed method. “Basic” denotes the basic constraint model described in Section 3.2 with at least one constraint in Section 4.1. “+C” denotes the basic constraint model with the cycle constraint discussed in Section 4.2. “+S” denotes the basic constraint model with the symmetry breaking constraint discussed in

<footnote n="1">¹https://doi.org/10.5281/zenodo.7460388</footnote>
<table><caption>Table 3: Comparisons between Proposed SAT-based Methods and the Existing Method</caption><thead><tr><th rowspan="2">Instance</th><th colspan="5">Statictics</th><th colspan="4">Proposed</th><th rowspan="2">Existing ASP</th></tr><tr><th>|Σ|</th><th>|T|</th><th>Attractors</th><th>k</th><th>|S<sub>max</sub>|</th><th>Basic</th><th>+C</th><th>+S</th><th>+C+S</th></tr></thead><tbody><tr><td>Example</td><td>4</td><td>12</td><td>1(3):2(1):4(1)</td><td>4</td><td>3</td><td>1.6</td><td>1.6</td><td>1.7</td><td>1.7</td><td><strong>0.2</strong></td></tr><tr><td>Lambda phage</td><td>4</td><td>46</td><td>1(1):2(1)</td><td>2</td><td>4</td><td>1.7</td><td>1.8</td><td>1.8</td><td>1.7</td><td><strong>0.1</strong></td></tr><tr><td>Trp-reg</td><td>4</td><td>14</td><td>1(2):4(1)</td><td>4</td><td>3</td><td>1.7</td><td>1.7</td><td>1.7</td><td>1.7</td><td><strong>0.2</strong></td></tr><tr><td>Fission-yeast</td><td>9</td><td>43</td><td>1(1)</td><td>1</td><td>3</td><td>1.4</td><td>1.4</td><td>1.5</td><td>1.4</td><td><strong>&lt;0.1</strong></td></tr><tr><td>Mamm.</td><td>10</td><td>34</td><td>1(1)</td><td>1</td><td>2</td><td>1.3</td><td>1.3</td><td>1.3</td><td>1.3</td><td><strong>&lt;0.1</strong></td></tr><tr><td>Tcrsig</td><td>40</td><td>85</td><td>1(8)</td><td>1</td><td>2</td><td>1.6</td><td>1.7</td><td>1.8</td><td>1.8</td><td><strong>&lt;0.1</strong></td></tr><tr><td>FGF</td><td>59</td><td>102</td><td>1(1536)</td><td>1</td><td>3</td><td>1.8</td><td>1.8</td><td>1.9</td><td>1.8</td><td><strong>&lt;0.1</strong></td></tr><tr><td>T-helper</td><td>101</td><td>316</td><td>1(5875504)</td><td>1</td><td>3</td><td>88.6</td><td><strong>86.6</strong></td><td>88.9</td><td>88.4</td><td>148.3</td></tr><tr><td>star02</td><td>2</td><td>4</td><td>3(1)</td><td>3</td><td>2</td><td>1.3</td><td>1.3</td><td>1.3</td><td>1.3</td><td><strong>0.1</strong></td></tr><tr><td>star03</td><td>3</td><td>6</td><td>4(1)</td><td>4</td><td>2</td><td>1.4</td><td>1.5</td><td>1.3</td><td>1.5</td><td><strong>0.3</strong></td></tr><tr><td>star04</td><td>4</td><td>8</td><td>5(1)</td><td>5</td><td>2</td><td>1.7</td><td>1.6</td><td>1.6</td><td>1.6</td><td><strong>0.6</strong></td></tr><tr><td>star05</td><td>5</td><td>10</td><td>6(1)</td><td>6</td><td>2</td><td>2.0</td><td>1.9</td><td>1.7</td><td>1.8</td><td><strong>1.6</strong></td></tr><tr><td>star06</td><td>6</td><td>12</td><td>7(1)</td><td>7</td><td>2</td><td>2.6</td><td>2.2</td><td><strong>2.1</strong></td><td><strong>2.1</strong></td><td>3.6</td></tr><tr><td>star07</td><td>7</td><td>14</td><td>8(1)</td><td>8</td><td>2</td><td>6.5</td><td>2.8</td><td>2.4</td><td><strong>2.3</strong></td><td>10.7</td></tr><tr><td>star08</td><td>8</td><td>16</td><td>9(1)</td><td>9</td><td>2</td><td>46.5</td><td>5.9</td><td>2.9</td><td><strong>2.7</strong></td><td>27.1</td></tr><tr><td>star09</td><td>9</td><td>18</td><td>10(1)</td><td>10</td><td>2</td><td>538.6</td><td>36.7</td><td>3.4</td><td><strong>3.2</strong></td><td>65.2</td></tr><tr><td>star10</td><td>10</td><td>20</td><td>11(1)</td><td>11</td><td>2</td><td>6598.8</td><td>442.7</td><td>4.1</td><td><strong>4.0</strong></td><td>148.0</td></tr><tr><td>star11</td><td>11</td><td>22</td><td>12(1)</td><td>12</td><td>2</td><td>T.O.</td><td>8787.0</td><td><strong>4.6</strong></td><td><strong>4.6</strong></td><td>1088.2</td></tr><tr><td>star12</td><td>12</td><td>24</td><td>13(1)</td><td>13</td><td>2</td><td>T.O.</td><td>T.O.</td><td>5.9</td><td><strong>4.8</strong></td><td>9959.3</td></tr><tr><td>star13</td><td>13</td><td>26</td><td>14(1)</td><td>14</td><td>2</td><td>T.O.</td><td>T.O.</td><td>7.8</td><td><strong>5.2</strong></td><td>T.O.</td></tr><tr><td>star14</td><td>14</td><td>28</td><td>15(1)</td><td>15</td><td>2</td><td>T.O.</td><td>T.O.</td><td>12.1</td><td><strong>5.9</strong></td><td>T.O.</td></tr><tr><td>star15</td><td>15</td><td>30</td><td>16(1)</td><td>16</td><td>2</td><td>T.O.</td><td>T.O.</td><td>19.8</td><td><strong>7.6</strong></td><td>T.O.</td></tr><tr><td>star16</td><td>16</td><td>32</td><td>17(1)</td><td>17</td><td>2</td><td>T.O.</td><td>T.O.</td><td>24.5</td><td><strong>9.7</strong></td><td>T.O.</td></tr><tr><td>star17</td><td>17</td><td>34</td><td>18(1)</td><td>18</td><td>2</td><td>T.O.</td><td>T.O.</td><td>38.8</td><td><strong>11.5</strong></td><td>T.O.</td></tr><tr><td>star18</td><td>18</td><td>36</td><td>19(1)</td><td>19</td><td>2</td><td>T.O.</td><td>T.O.</td><td>57.7</td><td><strong>14.5</strong></td><td>T.O.</td></tr><tr><td>star19</td><td>19</td><td>38</td><td>20(1)</td><td>20</td><td>2</td><td>T.O.</td><td>T.O.</td><td>87.7</td><td><strong>19.3</strong></td><td>T.O.</td></tr><tr><td>star20</td><td>20</td><td>40</td><td>21(1)</td><td>21</td><td>2</td><td>T.O.</td><td>T.O.</td><td>146.3</td><td><strong>28.5</strong></td><td>T.O.</td></tr><tr><td>star30</td><td>30</td><td>60</td><td>31(1)</td><td>31</td><td>2</td><td>T.O.</td><td>T.O.</td><td>6250.3</td><td><strong>614.4</strong></td><td>T.O.</td></tr><tr><td>star40</td><td>40</td><td>80</td><td>41(1)</td><td>41</td><td>2</td><td>T.O.</td><td>T.O.</td><td>T.O.</td><td><strong>7296.7</strong></td><td>T.O.</td></tr></tbody></table>

Section 4.3. “+C+S” denotes all constraints described in Section 3.2 and Section 4. The eleventh column denotes the CPU time of the existing ASP-based method.

In the biologically inspired benchmarks, all methods solved them within a few seconds. The main reason is that they contain small-sized attractors and thus *k* is also small. The difference between the proposed methods and the ASP-based method mainly comes from the implementation language: C++ and Scala. The proposed method is implemented on Scala which is running on a Java Virtual Machine (JVM) for the modeling part and thus there is a small disadvantage on this point. One exception is T-helper which contains a large number of attractors. On this problem, our approach takes advantage of the SAT solver *BDD_MINISAT_ALL* dedicated to the fast enumeration of solutions.

On *star* benchmarks, the difference between each method is more obvious. The CPU time of the ASP-based method increases exponentially, and it cannot

solve *star13* within 3 hours. On the other hand, although “Basic” cannot solve *star11*, the two optional constraints successfully improve “Basic”. The best results are obtained by combining those constraints.

# 7 CONCLUSIONS

This paper describes a SAT-based method for finding attractors of bounded size in asynchronous automata networks. Attractors are crucial to validate the initial design of a biological model and predict possible asymptotic behaviors, e.g., how cells may result through maturation in differentiated cell types. Given a bound *k*, we propose a constraint model of trap domains and compute attractors by incrementally computing them from 1 to *k*. We further propose two optional kinds of constraints to improve the computation time of our approach. The first one denotes cycle constraints which utilize a property of transitions in au-
tomata networks. The second one denotes symmetry breaking constraints which reduce redundant search spaces contained in the initial constraint model. Experimental evaluations are carried out over 30 automata networks. While the proposed SAT-based approach and the state-of-the-art ASP one could hardly be discriminated on the few existing biologically inspired benchmarks, their behavior is quite different on the crafted benchmarks we contribute. While the performance of the initial SAT-based model is not so effective, adding optional constraints allows it to scale much better than the ASP approach on benchmarks with controlled attractor size. Future work is listed as follows. Extending the comparisons to methods of other networks like Boolean network (Mori and Akutsu, 2022; Trinh et al., 2022; Rozum et al., 2021; Inoue, 2011) is necessary. Extending the proposed method to more complex networks like Petri net is interesting. Supplementary materials of this paper is available on the repository <fnref n="2" />.

## ACKNOWLEDGEMENTS

This work was financially supported by the “PHC Sakura” program (43009SC, JPJSBP120193213), implemented by the French Ministry for Europe and Foreign Affairs, the French Ministry of Higher Education, Research and Innovation and the Japan Society for Promotion of Science. This work was also supported by JSPS KAKENHI Grant Numbers JP20K11748 and JP20H05794.

## REFERENCES

Abou-Jaoudé, W., Monteiro, P. T., Naldi, A., Grandclaudon, M., Soumelis, V., Chaouiya, C., and Thieffry, D. (2015). Model checking to assess t-helper cell plasticity. *Frontiers in bioengineering and biotechnology*, 2:86.

Ben Abdallah, E., Folschette, M., Roux, O., and Magnin, M. (2017). Asp-based method for the enumeration of attractors in non-deterministic synchronous and asynchronous multi-valued networks. *Algorithms for Molecular Biology*, 12(1):1–23.

Biere, A., Fazekas, K., Fleury, M., and Heisinger, M. (2020). CaDiCaL, Kissat, Paracooba, Plingeling and Treengeling entering the SAT Competition 2020. In Balyo, T., Froleyks, N., Heule, M., Iser,

M., Järvisalo, M., and Suda, M., editors, *Proc. of SAT Competition 2020 – Solver and Benchmark Descriptions*, volume B-2020-1 of *Department of Computer Science Report Series B*, pages 51–53. University of Helsinki.

Biere, A., Heule, M., van Maaren, H., and Walsh, T., editors (2021). *Handbook of Satisfiability-Second Edition*, volume 336 of *Frontiers in Artificial Intelligence and Applications*. IOS Press.

Chai, X., Ribeiro, T., Magnin, M., Roux, O., and Inoue, K. (2020). Static analysis and stochastic search for reachability problem. *Electronic Notes in Theoretical Computer Science*, 350:139–158.

Crawford, J. M. and Baker, A. B. (1994). Experimental results on the application of satisfiability algorithms to scheduling problems. In *Proceedings of the 12th National Conference on Artificial Intelligence (AAAI 1994)*, pages 1092–1097.

Crawford, J. M., Ginsberg, M. L., Luks, E. M., and Roy, A. (1996). Symmetry-breaking predicates for search problems. In Aiello, L. C., Doyle, J., and Shapiro, S. C., editors, *Proceedings of the Fifth International Conference on Principles of Knowledge Representation and Reasoning (KR’96), Cambridge, Massachusetts, USA, November 5-8, 1996*, pages 148–159. Morgan Kaufmann.

Davidich, M. I. and Bornholdt, S. (2008). Boolean network model predicts cell cycle sequence of fission yeast. *PloS one*, 3(2):e1672.

de Kleer, J. (1989). A comparison of ATMS and CSP techniques. In *Proceedings of the 11th International Joint Conference on Artificial Intelligence (IJCAI 1989)*, pages 290–296.

Fauré, A., Naldi, A., Chaouiya, C., and Thieffry, D. (2006). Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle. *Bioinformatics*, 22(14):e124–e131.

Fitime, L. F., Roux, O., Guziolowski, C., and Paulevé, L. (2017). Identification of bifurcation transitions in biological regulatory networks using Answer-Set Programming. *Algorithms for Molecular Biology*, 12(1):19.

Folschette, M., Paulevé, L., Magnin, M., and Roux, O. (2015). Sufficient conditions for reachability in automata networks with priorities. *Theoretical Computer Science*, 608:66–83.

Garg, A., Di Cara, A., Xenarios, I., Mendoza, L., and De Micheli, G. (2008). Synchronous versus asynchronous modeling of gene regulatory networks. *Bioinformatics*, 24(17):1917–1925.

<footnote n="2">²https://github.com/TakehideSoh/bioinformatics2023-supplemental-material</footnote>
Gebser, M., Kaminski, R., Kaufmann, B., and Schaub, T. (2012). *Answer Set Solving in Practice*. Synthesis Lectures on Artificial Intelligence and Machine Learning. Morgan & Claypool Publishers.

Gelder, A. V. (2008). Another look at graph coloring via propositional satisfiability. *Discrete Applied Mathematics*, 156(2):230–243.

Gershenson, C. (2002). Classification of random boolean networks. In *Proceedings of the Eighth International Conference on Artificial Life*, pages 1–8.

Grieco, L., Calzone, L., Bernard-Pierrot, I., Radvanyi, F., Kahn-Perles, B., and Thieffry, D. (2013). Integrative modelling of the influence of mapk network on cancer cell fate decision. *PLoS computational biology*, 9(10):e1003286.

Huang, S., Eichler, G., Bar-Yam, Y., and Ingber, D. E. (2005). Cell fates as high-dimensional attractor states of a complex gene regulatory network. *Physical review letters*, 94(12):128701.

Inoue, K. (2011). Logic programming for boolean networks. In *IJCAI 2011, Proceedings of the 22nd International Joint Conference on Artificial Intelligence, Barcelona, Catalonia, Spain, July 16-22, 2011*, pages 924–930.

Iwama, K. and Miyazaki, S. (1994). SAT-variable complexity of hard combinatorial problems. In *Proceedings of the IFIP 13th World Computer Congress*, pages 253–258.

Kauffman, S. A. (1969). Metabolic stability and epigenesis in randomly constructed genetic nets. *Journal of Theoretical Biology*, 22(3):437–467.

Khaled, T. and Benhamou, B. (2020). An asp-based approach for boolean networks representation and attractor detection. In *LPAR*, pages 317–333.

Klamt, S., Saez-Rodriguez, J., Lindquist, J., Simeoni, L., and Gilles, E. (2006). A methodology for the structural and functional analysis of signaling and regulatory networks. *BMC Bioinformatics*, 7(1):56.

Levy, N., Naldi, A., Hernandez, C., Stoll, G., Thieffry, D., Zinovyev, A., Calzone, L., and Paulevé, L. (2018). Prediction of Mutations to Control Pathways Enabling Tumour Cell Invasion with the CoLoMoTo Interactive Notebook (Tutorial). *Frontiers in Physiology*, 9:787.

Mbodj, A., Junion, G., Brun, C., Furlong, E. E., and Thieffry, D. (2013). Logical modelling of drosophila signalling pathways. *Molecular BioSystems*, 9(9):2248–2258.

Mori, T. and Akutsu, T. (2022). Mini review attractor detection and enumeration algorithms for

boolean networks. *Computational and Structural Biotechnology Journal*.

Naldi, A., Hernandez, C., Levy, N., Stoll, G., Monteiro, P. T., Chaouiya, C., Helikar, T., Zinovyev, A., Calzone, L., Cohen-Boulakia, S., et al. (2018). The colomoto interactive notebook: accessible and reproducible computational analyses for qualitative biological networks. *Frontiers in physiology*, 9:680.

Paulevé, L. (2016a). Goal-oriented reduction of automata networks. In *International Conference on Computational Methods in Systems Biology*, volume 9859 of *Lecture Notes in Bioinformatics*, pages 252–272. Springer.

Paulevé, L. (2016b). Goal-oriented reduction of automata networks. In *International Conference on Computational Methods in Systems Biology*, pages 252–272. Springer.

Paulevé, L. (2017). Pint: a static analyzer for transient dynamics of qualitative networks with IPython interface. In *CMSB 2017 - 15th conference on Computational Methods for Systems Biology*, volume 10545 of *Lecture Notes in Computer Science*, pages 309–316. Springer International Publishing.

Paulevé, L., Andrieux, G., and Koepp, H. (2013). Under-approximating cut sets for reachability in large scale automata networks. In *International Conference on Computer Aided Verification*, pages 69–84. Springer.

Paulevé, L., Magnin, M., and Roux, O. (2012). From the Process Hitting to Petri Nets and Back. Technical Report hal-00744807, ETH Zürich.

Rougny, A., Paulevé, L., Teboul, M., and Delaunay, F. (2021). A detailed map of coupled circadian clock and cell cycle with qualitative dynamics validation. *BMC bioinformatics*, 22(1):1–24.

Rozum, J. C., Deritei, D., Park, K. H., Gómez Tejeda Zañudo, J., and Albert, R. (2021). pystablemotifs: Python library for attractor identification and control in Boolean networks. *Bioinformatics*, 38(5):1465–1466.

Sahin, Ö., Fröhlich, H., Löbke, C., Korf, U., Burmester, S., Majety, M., Mattern, J., Schupp, I., Chaouiya, C., Thieffry, D., et al. (2009). Modeling erbb receptor-regulated g1/s transition to find novel targets for de novo trastuzumab resistance. *BMC systems biology*, 3(1):1–20.

Simão, E., Remy, E., Thieffry, D., and Chaouiya, C. (2005). Qualitative modelling of regulated metabolic pathways: application to the tryptophan biosynthesis in e. coli. *Bioinformatics*, 21(suppl 2):ii190–ii196.
Soh, T., Banbara, M., and Tamura, N. (2017). Proposal and evaluation of hybrid encoding of CSP to SAT integrating order and log encodings. *International Journal on Artificial Intelligence Tools*, 26(1):1–29.

Tamura, N., Taga, A., Kitagawa, S., and Banbara, M. (2009). Compiling finite linear CSP into SAT. *Constraints*, 14(2):254–272.

Thieffry, D. and Thomas, R. (1995). Dynamical behaviour of biological regulatory networks—ii. immunity control in bacteriophage lambda. *Bulletin of Mathematical Biology*, 57(2):277–297.

Thomas, R. (1973). Boolean formalization of genetic control circuits. *Journal of Theoretical Biology*, 42(3):563–585.

Toda, T. and Soh, T. (2016). Implementing efficient all solutions SAT solvers. *ACM J. Exp. Algorithms*, 21(1):1.12:1–1.12:44.

Trinh, V.-G., Hiraishi, K., and Benhamou, B. (2022). Computing attractors of large-scale asynchronous boolean networks using minimal trap spaces. In *Proceedings of the 13th ACM International Conference on Bioinformatics, Computational Biology and Health Informatics*, pages 1–10.

Tseitin, G. S. (1968). On the complexity of derivations in the propositional calculus. *Studies in Mathematics and Mathematical Logic Part II*, pages 115–125.

Walsh, T. (2000). SAT v CSP. In *Proceedings of the 6th International Conference on Principles and Practice of Constraint Programming (CP 2000)*, pages 441–456.

Wuensche, A. (1998). Genomic regulation modeled as a network with basins of attraction. In *Pacific Symposium on Biocomputing*, volume 3, pages 89–102.