A SAT-based Method for Counting All Singleton Attractors in Boolean Networks

1 Introduction

Boolean network (BN) [Kauffman, 1969b; Kauffman, 1993] is a mathematical model representing state transition systems. In a BN, each node is assigned a Boolean value of 0 (inactive) or 1 (active). The state of each node is then determined by a Boolean function based on the states of its input nodes in the previous discrete time step. The analysis of both static and dynamic properties of transition systems using BN has recently gained increasing importance in various fields. In particular, in biological systems, its application includes gene regulatory networks [Kauffman, 1969a; Daizhan Cheng, 2011], neural networks [Anthony, 2001], cell cycle control networks [M¨unzner et al., 2019], and signal transduction networks [Mori et al., 2015]. With the expanding range of such applications, repositories like GINsim [Naldi et al., 2018], CellCollective [Helikar et al., 2012] and Biomodels [Malik-Sheriff et al., 2020], which describe biological networks in a formal format like BN, are being published and maintained.

Attractor is a steady state that a BN eventually reaches as time progresses. It can be either a singleton attractor (consisting of a single state) or a periodic attractor (cycling through multiple states). Attractors may also represent homeostasis, ensuring that the network can return to a steady state after a perturbation. Given its importance in understanding the system, this has led to extensive research on attractor detection [Mori and Akutsu, 2022]. In particular, analysis of singleton attractors (a.k.a. fixed points) plays a crucial role in understanding complex biological models since many biological systems stabilize to an observable phenotype rather than exhibit cyclic behavior. For instance, in gene regulatory networks, singleton attractors represent different cell types or states [Kauffman, 1993; Abou-Jaoud´e et al., 2016], thus helping understand how cells differentiate.

Attractor computation is thus an important subject. There are mainly two major approaches: BDD- and ASP-based methods [Garg et al., 2008; Su and Pang, 2021; Klarner et al., 2017; Benes et al., 2022; Trinh et al., 2023]. Those studies focused on enumerating attractors. However, recently published BNs are large and contain a huge number of attractors. For instance, the BN named "#124 NSP9-PROTEIN" consists of 252 nodes and has over 1020 singleton attractors. For such BNs, enumerating all attractors is neither realistic nor particularly useful. Indeed, there are BNs for which all attractors cannot be enumerated using state-of-the-art tools.

This limitation suggests counting attractors instead of enumerating them, as the number of attractors is still useful for assessing effects of intentional gene activation/inhibition or providing insights to network modelers. Thus, recent enumeration tools provide an option that returns only the number of attractors to reduce I/O costs, enabling retrieval of attractor count without displaying them. However, even with this feature, there remain BNs for which the number of attractors is still unknown. It is also challenging to thoroughly revise existing BDD- or ASP-based methods to adopt pure counting approaches. In the case of BDDs, a BDD is constructed as a compressed representation of all attractors, making enumeration and counting essentially equivalent. For ASP, it would be feasible if enumeration of stable models could be replaced with counting. However, the latest ASP encodings [Trinh et al., 2023] include disjunctive rules, whereas the latest stable model counter [Kabir et al., 2024] does not support disjunctive rules, making this approach difficult to implement.

Boolean Satisfiability Testing (SAT) problem [Biere et al., 2021] is the problem of determining if there exists an assignment of variables that satisfies a given Boolean formula. SAT

has widespread applications in areas such as model checking, scheduling, automated planning, constraint satisfaction, and more [Biere et al., 2021]. Due to its significance, SAT has garnered extensive research attention, leading to the development of various efficient SAT solvers, whose input is formulas in Conjunctive Normal Form (CNF). Advances in SAT technologies have extended their applicability to a range of challenging and practical combinatorial problems. Furthermore, SAT has several variants, such as AllSAT and #SAT whose goals are to enumerate or count all models of a given formula. In particular, for SAT solvers and #SAT solvers, there are international competitions every year 1 , and these competitions witness the improvement of the performance of SAT and #SAT solvers.

In this paper, we propose a SAT-based method for counting all singleton attractors of a given BN. To fully carry the power of the state-of-the-art solvers, the difficulty is twofold: the Boolean formulation of singleton attractors such that models of the formula and attractors have one-to-one correspondences; the translation of arbitrary formulas into CNF formulas to allow the use of standard SAT or #SAT solvers. To tackle this, we develop a Boolean formulation of singleton attractors and study three CNF translation methods, including one using the classical Tseitin translation [Tseitin, 1968], one using a novel small implicant translation, and the last one using the hybridization of the other two approaches. We evaluated the performance of the proposed method with the stateof-the-art tools for BNs by comparing them with the 624 instances used in the literature [Trinh et al., 2023] and 19 additional instances obtained from the BBM repository updates. The contributions of this paper are summarized as follows:

• A Boolean formulation (not in CNF) for singleton attractors in BNs whose models and attractors have oneto-one correspondence.
• The evaluation of three translation methods from this formulation to CNF formulas, including Tseitin, dedicated small implicants, and their hybridization.
• Comprehensive comparisons with state-of-the-art tools using 243 real-world and 400 artificial BNs.
• Determination of the number of singleton attractors for 9 real-world BNs and 73 large, crafted BNs, both of which have never been revealed before.

2 Preliminaries

2.1 Boolean Logic Recap

We start with the definition of a general form of Boolean formulas. Hereafter, unless otherwise specified, Boolean formula will be referred to as formula. A Boolean variable v is a formula. If P is a formula, ¬P is also a formula. If P and Q are formulas, then P ∆ Q are formulas, where ∆ can be Boolean operators ∨, ∧,↔, and →. A disjunction of literals is called a clause. A conjunction of literals is called a term. A formula is in Conjunctive Normal Form (CNF) if it is a conjunction of one or more clauses. A formula is in Disjunctive Normal Form (DNF) if it is a disjunction of one or more terms. A formula is in Negation Normal Form (NNF) if it contains only the logical connectives ∧, ∨, and ¬, and the ¬ operator applies only to literals. A partial assignment is a mapping from a subset of Boolean variables to the Boolean values 0 (false) or 1 (true). A complete assignment (or simply assignment) is a mapping from the set of Boolean variables to the Boolean values 0 (false) or 1 (true). We use the symbol α to denote assignments and extend it to formulas, defining a mapping that takes a formula and returns its truth value. A Boolean formula is satisfiable (SAT) if there exists an assignment that satisfies it and is unsatisfiable (UNSAT) otherwise. A partial model (or complete model) is a partial assignment (or complete assignment) that satisfies a given formula. A SAT solver is a program that computes a model (if any) from a CNF formula. A model counter/enumerator is a SAT solver specialized for counting/enumerating all models from a CNF formula.

1 https://satcompetition.github.io/, https://mccompetition.org/

2.2 Boolean Network and Attractor

Definition 1. A Boolean network (V, F) is defined as follows. V = {v1, v2, . . . , vn} is a set of Boolean variables, also referred to as nodes. F = (f1, f2, . . . , fn) is a vector of Boolean formulas, where fi is called an update function for the node vi. The update function fi : {0, 1}k(i) → {0, 1} is a Boolean formula over Boolean variables Vi = {vi,1, vi,2, . . . , vi,k(i)} ⊆ V , where k is a function that returns the number of variables contained in fi. We call the set of variables Vi as the set of input variables of vi.

When a BN represents a gene regulatory network, each vi corresponds to each gene and values 0 and 1 correspond to the inactive and active states of genes. Suppose that the states of a BN evolve over discrete time points {0, 1, 2, . . .}. Then, the state of a BN at a discrete time point t ∈ {0, 1, 2, . . .} is given by a vector →v(t) = (vt1 , vt2 , . . . , vtn ), where vt i represents the state of the node vi at time t. The value of the node vi at time t + 1 is determined by the update function fi applied to the values of its input nodes at time t as follows. Here we introduce f t i to simplify the description of fi.

Several update schemes were studied: synchronous (all nodes updated simultaneously), asynchronous (one node at a time), and generalized semantics (any subset of nodes updated simultaneously). These different semantics preserve singleton attractors [Paulev´e and Sen´e, 2022]. As such, the discussion proceeds with asynchronous BNs without loss of generality. We define the relation of asynchronous transition T between two states →v(t) and →v(t + 1) as follows.

This denotes that the value of at most one variable vi can change between two states in the asynchronous transitions.

Definition 2. A State Transition Graph (STG) (S, T) of a BN (V, F) is defined as follows where S = {0, 1}n is the set of states, and T is the transition relation given by the BN. The state transition graph is a directed graph consisting of the set of states S as nodes, and the arcs interpreted from the transition T .

Definition 3. A Singleton Attractor is a state that cannot transition to any other state but itself. In other words, a singleton attractor is a state that does not have outgoing arcs to any other node in STG.

Figure 1 illustrates a BN example and its STG, where 110 (red colored) is the only singleton attractor. Note that selftransition-loops are omitted in this Figure. Then, our target counting problem is defined as follows:

Definition 4. Given BN (V, F ), a Singleton Attractor Counting Problem (SACP) is the problem of counting the number of all singleton attractors of (V, F).

3 RelatedWork and Existing Tools

There have been various tools that can compute attractors in BNs. Boolsim [Garg et al., 2008] uses algorithms based on reduced ordered binary decision diagrams (ROBDD) and is capable of computing singleton and periodic attractors. Cabean [Su and Pang, 2021] is also a ROBDD-based method. PyBoolNet [Klarner et al., 2017] uses the popular model checker NuSMV and ASP solvers. AEON [Benes et al., 2022] uses ROBDD and is scalable to compute more than a trillion singleton attractors. fASP [Trinh et al., 2023] uses an ASP-based method to enumerate singleton attractors.

Automata network (AN) [Folschette et al., 2015] is another mathematical model for state transition systems and is a subclass of safe Petri nets [Paulev´e et al., 2012]. AN-ASP [Abdallah et al., 2017] is an ASP-based method to calculate the attractors in AN. SAF [Soh et al., 2023a] is a tool employing a SAT-based method [Soh et al., 2023b] to calculate attractors in AN. While a BN can be represented as an AN via transformations using tools like BioLQM [Naldi, 2018], the transition rules in automata networks are different, and it is not possible to specify complex logical expressions in the firing conditions.

At the time of writing, no methods or tools focus solely on counting attractors, likely due to technical challenges (Section 1) and the adequacy of enumeration for medium-sized BNs. However, larger and more complex BNs have made enumeration insufficient, motivating our research. ASPbased methods, the state of the art for attractor enumeration, could support counting if a stable model counter were available. In 2024, sharpASP [Kabir et al., 2024] was introduced, but its lack of support for disjunctive programs limits applicability to methods like fASP. If this is addressed, ASP-based methods could become strong candidates.

4 ProposedMethods

One approach to constructing our proposed SAT-based singleton attractor counter is to utilize existing SAT-based methods. As mentioned in the previous section, there is a SATbased method for ANs [Soh et al., 2023a], and by applying the BN-to-AN conversion tool BioLQM along with the latest model counter before and after this method, a SATbased counter can be realized. However, preliminary experiments showed that this approach did not even outperform the best enumeration-based methods (see Section 5). Other approaches have applied Boolean formulas to AND/OR BNs (i.e., BNs consisting of AND/OR nodes), assuming that update functions are given in DNF. While they are theoretically capable of handling arbitrary functions, they do not provide concrete procedures for doing so [Tamura and Akutsu, 2009; Inoue, 2011]. In this section, we thus propose a constraint model for counting all singleton attractors in BNs with arbitrary update function and present methods to translate this model into CNF compactly.

4.1 Boolean Formulation of Singleton Attractors

We start with the Boolean formulation that must transition to any other state but itself. Suppose that vi is the node whose value changes in the transition between t and t+1. If its value changes from 0 (at t) to 1 (at t + 1) then ¬vt i ∧ f t i holds at time t. Similarly, if its value changes from 1 (at t) to 0 (at t + 1) then vt i ∧ ¬f t i holds at time t. By negating those two constraints, we obtain the constraint for singleton attractors at time t as follows.

Since singleton attractors can be detected at any time point t, we arbitrarily consider the time t = 0. Thus, for simplicity, we omit the notation of t . Finally, we obtain the Boolean formula Ψ as follows.

Proposition 1. Given BN (V, F), any complete model of the formula Ψ can be interpreted as a singleton attractor of (V, F ) and vice versa.

Proof. (⇒) Let α be a model of Ψ consisting of a Boolean value vector →s = (s1, s2, . . . sn) such that α(vi) = si. It does not satisfy the two transition conditions (¬vi ∧ fi) and (vi ∧¬fi). Then, the state →s cannot transition to any other state but itself. Thus, →s is an attractor of (V, F). (⇐) If a value vector →s is an attractor of (V, F). Then, the state →s does not transition to any other state but itself. A model α consisting of →s does not satisfy the two transition conditions (¬vi ∧ fi) and (vi ∧ ¬fi). Thus, α is a model of Ψ.

To utilize the power of SAT techniques, we need to translate the arbitrary Boolean functions fi into a CNF formula. In the following sections, we evaluate three methods to translate Ψ into a CNF formula.

4.2 Direct Translation

The most common way to translate an arbitrary Boolean function in CNF is to use the Tseitin translation [Tseitin, 1968]. Direct translation uses bidirectional Tseitin translation to convert the arbitrary formula Ψ into a CNF formula ΨD. It is performed by recursively replacing any "complex" subformula by a new variable and appending the definition of the new variable to the formula. The process terminates when the formula is either a conjunction or a disjunction of literals. Those definitions can be simply expressed by sets of clauses without adding new variables. We denote the Tseitin translation by TT. Using the Tseitin translation, Direct translation is given as follows.

It is known that bidirectional Tseitin translation preserves the number of complete models, i.e. the original formula and the translated CNF have exactly the same number of complete models. Thus, the following proposition holds.

Proposition 2. Given BN (V, F), any complete model of the formula ΨD restricted to the original variables of Ψ is a complete model of Ψ and any complete model of Ψ can be extended to a complete model of ΨD.

Example 1. Suppose that f1 = (v2 ∨ v3 ∨ v4) ∧ (v2 ∨ ¬v5) is an update function of a given BN. By Direct translation described above, we obtain the following Boolean formula, where l1, p1 and p2 are newly introduced auxiliary Boolean variables.

This Direct translation can be implemented without any complex computation. However, the disadvantage of this translation is that sometimes the number of auxiliary variables and companion clauses becomes comparatively large, which can have a negative impact on the solver runtime as we will see in the experimental section.

4.3 Indirect Translation

In this approach, we will use again Tseitin traslation, but to feed the Boolean function to a model enumerator to obtain a DNF representation without adding new variables.

Let Mbe a function that translates a Boolean formula into the set of all its (complete) models. Then, for a formula ϕ, we can obtain a logically equivalent formula in disjunctive normal form (DNF) by using a DNF translation function DNF as follows.

Then, using this translation, we can translate the singleton attractor formula Ψ into a CNF formula as follows, since the negation of a DNF is a CNF and the conjunction of two CNF is still a CNF:

Then, the following proposition holds.

Proposition 3. Given BN (V, F ), any complete model of formula ΨI can be interpreted as a complete model of Ψ and vice versa.

To perform this conversion using an off-the-shelf model enumerator, we use again the Tseitin translation. The models of TT(fi) contain both auxiliary variables from the Tseitin translation and the original variables of fi. We can just filter out the auxiliary variables to get all the models of fi from all the models of TT (fi).

However, the number of literals in each clause corresponds to the number of input variables for the function fi, and the number of complete models may be exponential in the number of input variables of fi. This can lead to ΨI with a large number of lengthy clauses, which may prevent the translation into DNF to complete or produce a formula too large for the SAT solvers or model counters. An alternative approach is to compute the prime implicants of the formula, i.e. subset-minimal sets of literals which satisfy the formula, instead of its models. We employ the translation-based method proposed in [Jabbour et al., 2014]: each model of the translated CNF ψPI corresponds to a prime implicant of the original CNF ψ. Algorithm 1 shows the final procedure PModelEnum. We call PModelEnum(ψ,∞) for implementing the aforementioned functionM.

The prime implicant of a Tseitin translated formula is not necessarily a prime implicant of the original formula. Take, for instance, the formula f = (v1 ∧ v2) ∨ (v1 ∧ ¬v2). It will be translated using the Tseitin translation as TT(f) = (p1 ∨ p2) ∧ (p1 ↔ v1 ∧ v2) ∧ (p2 ↔ v1 ∧ ¬v2). The single prime implicant of f is v1. The prime implicants of TT(f) are p1 ∧ ¬p2 ∧ v1 ∧ v2 and ¬p1 ∧ p2 ∧ v1 ∧ ¬v2, which gives v1 ∧ v2 and v1 ∧ ¬v2 when restricted to the original variables, which is exactly the set of models, i.e. M(f). So computing the prime implicants of the Tseitin translated formula is a pragmatic effort to reduce the size of the DNF encoding. Note that an approach to compute exactly the PI of non clausal Boolean formula also exists [Previti et al., 2015]. We implemented PModelEnum in Scala and use Sat4j [Le Berre and Parrain, 2010] as ModelEnumerator to perform all those calls to model enumerators in memory: for most benchmarks, the procedure takes less than 10 seconds.

Since we are interested in computing the prime implicants of our Boolean function fi, we can simplify the Tseitin translation by only using an implication from the auxiliary variables to the original variables (instead of a bi-implication) and apply a negation normal form (NNF) [Darwiche and Marquis, 2002] to avoid negated auxiliary variables. This will reduce the number of clauses, not the number of variables, and will also affect the number of prime implicants found. Take for instance f = (v1∧v2)∨(v3∧v4), which has 2 prime implicants (v1 ∧v2 and v3 ∧v4). The one-direction Tseitin encoding will produce 1T T (f ) = (p1∨p2)∧(p1 → v1∧v2)∧(p2 → v3∧v4)

Algorithm 1 PModelEnum

Input: formula ϕ, cutoff co

Output: if enumeration succeeded or not, models M

• 1: Ω := one-direction Tseitin translation of NNF(ϕ)
• 2: ΩPI := PI translation [Jabbour et al., 2014] of Ω
• 3: (isOK,M) :=ModelEnumerator(ΩPI, co)
• 4: M' := {decode(m) | m ∈M}
• 5: return (isOK,M ′ )

which has two prime implicants, p1 ∧ ¬p2 ∧ v1 ∧ v2 and ¬p1 ∧ p2 ∧ v3 ∧ v4 while TT(f) has four prime implicants, p1 ∧ ¬p2 ∧ v1 ∧ v2 ∧ ¬v3, p1 ∧ ¬p2 ∧ v1 ∧ v2 ∧ ¬v4, ¬p1 ∧ p2 ∧ ¬v1 ∧ v3 ∧ v4 and ¬p1 ∧ p2 ∧ ¬v2 ∧ v3 ∧ v4.

Example 2. Suppose that f1 = (v2 ∨v3 ∨v4)∧(v2 ∨¬v5) is an update function. Using Algorithm 1, the following prime implicants of ¬v1 ∧ f1 and v1 ∧ ¬f1 are obtained:

Then, we obtain the following CNF formula as ΨI .

The upper three clauses are for ¬v1 ∧ f1, and the lower two clauses are for v1 ∧ ¬f1.

So an alternative way to obtain a DNF of ϕ is as follows:

While not guaranteed that DNF P (ϕ) represents the complete set of prime implicants, the impact of model compression is significant. For instance, we confirmed that for a BN rule, the number of clauses was reduced from 1,073,741,823 to 30 within a reasonable translation time, compared to the naive model-enumeration-based translation DNF (ϕ).

4.4 Proofs for Indirect translation

In this section, we will prove DNF(ϕ) ↔ DNF P (ϕ). To provide the proof, we start with the two properties of the Tseitin transformation in one direction. The full version of proof is available online 2 .

Properties of One-direction Tseitin Translation. Let A[B] be a source negation normal form (NNF) formula containing a subformula B. Note that an arbitrary formula can be converted to the NNF formula by using De Morgan's law. Let p be an auxiliary variable that does not appear in A[B]. Let A be an NNF formula such that A is obtained by all the appearances of B with the positive literal p. Then, a one-directional Tseitin translation of A[B] for B can be performed by A ∧ (¬p ∨ B). We clarify its properties with the following two propositions, which is based on the proof for equi-satisfiability in the literature [Tamura et al., 2010].

Proposition 4. Suppose that α is a partial model of A∧(¬p∨ B), and α′ is a partial model obtained by excluding p from α. Then α′ (A[B]) = 1 holds.

2 https://doi.org/10.5281/zenodo.15462044

Proposition 5. Let α be a model such that α(A[B]) = 1. Let α′ be a model that extends α to the domain that includes p such that α′ (p) = α(B). Then α′ (A ∧ (¬p ∨ B)) = 1 holds.

Properties of PModelEnum. From here, we will also represent an assignment α as a set of literals for simplicity. For instance, suppose that ϕ is a formula over variables p, q, r. Then, an assignment p = 1, q = 1, r = 0 will be represented as {p, q, ¬r}. A partial assignment where p = 1, q = 1 will be represented as {p, q}. In addition, α(ϕ) = 1 means that any assignment satisfies ∧l∈α l also satisfies ϕ. The following proposition holds between M(ϕ) and PModelEnum(ϕ, ∞).

Proposition 6. Given an arbitrary formula ϕ, for any α ∈ M(ϕ), there exists a partial model α′ ∈ PModelEnum(ϕ,∞) such that α′ ⊆ α and α′ (ϕ) = 1.

Proof. We prove it by contradiction. Assume that for some α ∈ M(ϕ), there does not exist a partial model α′ ∈ PModelEnum(ϕ,∞) such that α′ ⊆ α and α′ (ϕ) = 1. For such α, by Proposition 5, there exists a model β (which is α extended with auxiliary variables) that is a model of 1TT(ϕ). Then, let β′ be a prime implicant of 1TT(ϕ). Obviously, β′ ⊆ β and β′ (1TT(ϕ)) = 1 hold. For such β′ , by Proposition 4, there exists a model α′ that satisfies α′ (ϕ) = 1. Clearly, α′ ⊆ α holds, which contradicts the assumption. Therefore, for any α ∈ M(ϕ), there exists a partial model α′ ∈ PModelEnum(ϕ,∞) such that α′ ⊆ α and α′ (ϕ) = 1.

Then, we get the following proposition that gives the logical equivalence between DNF (ϕ) and DNF P (ϕ).

Proposition 7. Given an arbitrary formula ϕ, any model of DNF(ϕ) is also a model of DNF P (ϕ) and vice versa.

Proof. Let Φ1 be DNF(ϕ) and Φ2 be DNF P (ϕ). (⇒) According to Proposition 6, for any α ∈ M(ϕ) there exists α′ ∈ PModelEnum(ϕ,∞) such that α′ ⊆ α. Thus, if Φ1 has a satisfied disjunct, then Φ2 also has a satisfied disjunct. Therefore, any model of Φ1 is also a model of Φ2. (⇐) According to Proposition 6, α′ (ϕ) = 1. Therefore, there exists a complete model α ∈ M(ϕ), which is obtained by extending α′ to a complete model. Therefore, any model of Φ2 is also a model of Φ1.

4.5 Hybrid Translation

Direct translation could generate shorter clauses, but it needs many auxiliary variables, which explode the search space. Indirect translation does not need any auxiliary variables but could generate an exponential number of lengthy clauses. Therefore, we develop a hybrid method that combines the strengths of both encodings as follows:

where HybridTranslation is detailed in Algorithm 2 and co is a cutoff value. The two translation methods are evaluated for each update function. The switching criterion is the number of models enumerated. If it becomes greater than co, Direct translation is used; otherwise, Indirect translation is used.

Algorithm 2 HybridTranslation

Input: vi, fi, co

Output: CNF formula

• 1: (isOK+,M +) := PModelEnum(¬vi ∧ fi, co)
• 2: (isOK-,M -) := PModelEnum(vi ∧¬fi, co)
• 3: if (isOK+ & isOK-) then return ¬(M + ∪ M -)
• 4: else return (vi ∨¬li) ∧ (¬vi ∨ li) ∧ TT(li ↔ fi)

Since one of the two translation methods is used for each update function, the following proposition is obvious.

Proposition 8. Given BN (V, F), any complete model of the formula ΨH restricted to the original variables of Ψ is a complete model of Ψ and any complete model of Ψ can be extended to a complete model of ΨH.

Indirect translation generally performs better than Direct translation, but it requires model enumeration, which can be time-consuming. We explored various co values across some instances and found that 1000 strikes an effective balance, covering most BN rules within a reasonable translation time. Therefore, we used co = 1000 in our experiments.

5 Experiments

Environment. We conducted all experiments on a machine with Core i5 12400 (2.5GHz) of CPU and 64 GB of RAM. The time limit was 30 minutes per solver per benchmark.

Benchmark. We have used all the benchmarks employed in [Trinh et al., 2023], which consist of a wide variety of Boolean networks. Biodivine Boolean Models (BBM) [Pastva et al., 2023] include 230 networks sourced from all popular model repositories of Biological Systems such as GINsim, CellCollective and Biomodels, as well as individual publications. Pseudo-random (PRandom) consists of 400 random networks with 1000-5000 variables, preserving the structure of real-world networks, created in the literature [Trinh et al., 2023]. This benchmark provides larger BN instances than the real-world networks. The selected BNs (Selected) are 13 instances of real-world BNs that were collected in [Trinh et al., 2023].

Compared Solvers. As discussed in Section 3, since no attractor counter exists to the best of the authors' knowledge, we use the state-of-the-art enumeration tools for comparison: AEON [Benes et al., 2022], fASP-c and fASP-s [Trinh et al., 2023], PyBoolNet [Klarner et al., 2017], and SAF [Soh et al., 2023a]. The latter one uses BioLQM [Naldi, 2018] to translate each BN into AN. The runtime provided does include such translation time.

No Printing Options. To ensure the best option for counting, we configured all enumeration tools to not print attractors but only the number of attractors found, thereby avoiding I/O overhead. If we want to obtain all the solutions as a plain text file, the computational performance significantly slows down. For instance, AEON solved 224 instances without printing solutions within 1800 seconds, but it only solved 204 instances when printing solutions. This also indicates that the enumeration of real large BNs is impractical.

Proposed Method Setting. The three translation methods are implemented in Scala. We evaluated three state-of-the-art exact model counters from the 2023 competition, SharpSATTD [Korhonen and J¨arvisalo, 2023], d4 [Lagniez and Marquis, 2017], and GPMC [Ryosuke Suzuki and Sakai, 2017], on our SACP instances. SharpSAT-TD, the model-counting track winner, provided the best results. Additionally, we offer an enumeration solver using Hybrid translation and the AllSAT solver BDDMINISATALL [Toda and Soh, 2016]. We use D.C., I.C., and H.C. for Direct, Indirect, and Hybrid counting solvers, and H.E. for the Hybrid enumeration solver in the tables below. The runtime provided does include all translation times. Another proposed baseline-method is the one using the existing SAT-based method for ANs [Soh et al., 2023a]. We configured SAF with BioLQM and SharpSAT-TD enabling to count singleton attractors of BNs (named SAF.C).

Consistency and Resource Availability. We have checked all numbers of attractors are consistent between all solvers. All benchmarks, solvers, and log files are available online2 .

Experimental results. Table 1 summarized the numbers of instances solved by each solver that are classified by A) benchmark series, B) range of the number of attractors |A|, and C) range of the number of nodes n. The first and second columns denote each category and the number of instances it contains. The remaining columns denote the number of instances solved by existing methods and proposed methods. In A), the proposed translation-based method was able to count the number of singleton attractors in almost half of the instances (between 313 and 330). In contrast, about one third of them could be enumerated by AEON (244), H.E. (240), and fASP-s (235). In B), the proposed method solved the largest number of instances over all ranges. In particular, the proposed method is the only tool capable of counting all singleton attractors when BNs contain more than 1030 of them. In C), the proposed method solved the largest number of instances over all ranges. In particular, the proposed method is the only tool capable of computing all singleton attractors in BNs whose number of nodes is greater than 2000. Note that AEON and fASP-s solved, respectively, 6 and 9 instances in that range, but all of those are UNSAT, which means there are no singleton attractors. Additionally, real-world instances have fewer than 500 nodes, except for two in Selected. All P-Random instances have more than 1000 nodes.

Table 2 lists BN instances where the number of singleton attractors was previously unknown and presents the computed values using the proposed approach. The first column denotes the network name and its source literature. The second and third columns denote the number of variables and the number of attractors identified in BN, respectively. This table demonstrates the ability of the proposed method to uncover previously unknown attractors in BN, and also highlights a challenging benchmark instance.

Impact of the SAT Translation on CNF Formulas. Table 3 shows the comparisons between the three proposed translation methods: Direct, Indirect and Hybrid. The first column denotes the initial of the benchmark set. The 2nd to 4th columns show the number of instances that cannot be translated within the time limit of 1800 seconds. The 5th column shows the average ratio of the number of variables between Direct and Indirect translation. Similarly, the 6th column shows the average ratio of the number of variables between Direct and Hybrid translations. That is, if the figures are less than 1 then the number is smaller than the number of Direct translation. The following columns similarly denote the ratio of the number of clauses and the number of literals. We can see that Indirect translation cannot translate three instances within the time limit while Direct and Hybrid translation can translate all instances. We can also see that the number of variables and clauses of Indirect and Hybrid translation are smaller than the ones of Direct translation. However, the number of literals is larger than in Direct translation, except for BBM instances. The update functions in BBM instances often consist of a small number of variables but many operators such as (v1 ∧(v2 ∨v3)∧¬v4)∨(v2 ∧v3 ∧¬v4). In these cases, Direct translation generates more literals than Indirect or Hybrid translation.

Table 3: Number of timeouts by each translation and average ratios between Direct and Indirect/Hybrid translations

6 Conclusion

In this work, we proposed a SAT-based method for counting singleton attractors in large Boolean networks (BNs). It determines the number of singleton attractors for 9 real-world and 73 crafted BNs that had never been revealed before, a result that holds some impact in the field of biological network research. Furthermore, the translation methods for converting arbitrary propositional formulas into compact CNF are considered applicable to other problems in the field of AI. As demonstrated by the results of the proposed baseline solver, our findings are not merely the outcome of applying a model counter but are instead a result of the careful refinements and techniques introduced in our translation approach. Future work is refinements on the algorithm, including leveraging existing logic minimizers, such as Espresso 3 , to further reduce the DNFs obtained through Indirect translation.

3 https://github.com/Gigantua/Espresso

Acknowledgements

This work was initiated under the collaborative research supported by the "PHC Sakura" program (43009SC, JPJSBP120193213), implemented by MESRI and JSPS, which ended in 2021. This work was also supported by JSPS KAKENHI (JP23K11047, JP24H00686, JP25K03097), and by ROIS NII Open Collaborative Research 2023 (23FP04).