ASPbased method for the enumeration of attractors in nondeterministic synchronous and asynchronous multivalued networks
 Emna Ben Abdallah^{1}Email authorView ORCID ID profile,
 Maxime Folschette^{2, 3},
 Olivier Roux^{1} and
 Morgan Magnin^{1, 4}
https://doi.org/10.1186/s1301501701112
© The Author(s) 2017
Received: 1 April 2016
Accepted: 26 July 2017
Published: 15 August 2017
Abstract
Background
This paper addresses the problem of finding attractors in biological regulatory networks. We focus here on nondeterministic synchronous and asynchronous multivalued networks, modeled using automata networks (AN). AN is a general and wellsuited formalism to study complex interactions between different components (genes, proteins,...). An attractor is a minimal trap domain, that is, a part of the statetransition graph that cannot be escaped. Such structures are terminal components of the dynamics and take the form of steady states (singleton) or complex compositions of cycles (nonsingleton). Studying the effect of a disease or a mutation on an organism requires finding the attractors in the model to understand the longterm behaviors.
Results
We present a computational logical method based on answer set programming (ASP) to identify all attractors. Performed without any network reduction, the method can be applied on any dynamical semantics. In this paper, we present the two most widespread nondeterministic semantics: the asynchronous and the synchronous updating modes. The logical approach goes through a complete enumeration of the states of the network in order to find the attractors without the necessity to construct the whole statetransition graph. We realize extensive computational experiments which show good performance and fit the expected theoretical results in the literature.
Conclusion
The originality of our approach lies on the exhaustive enumeration of all possible (sets of) states verifying the properties of an attractor thanks to the use of ASP. Our method is applied to nondeterministic semantics in two different schemes (asynchronous and synchronous). The merits of our methods are illustrated by applying them to biological examples of various sizes and comparing the results with some existing approaches. It turns out that our approach succeeds to exhaustively enumerate on a desktop computer, in a large model (100 components), all existing attractors up to a given size (20 states). This size is only limited by memory and computation time.
Keywords
Background
In the last decades, the emergence of a wide range of new technologies have made it possible to produce a massive amount of biological data (genomics, proteomics...). This leads to considerable developments in systems biology which takes profit from this data. In order to understand the nature of a cellular function or more broadly a living biological system (healthy or diseased), it is indeed essential to study not only the individual properties of cellular components, but also their interactions. The behavior and functionalities of the cells emerge from such networks of interactions.
Considering this paradigm, the longterm behavior of regulatory networks dynamics is of specific interest [1]. Indeed, at any moment, a system may fall into a trap domain, which is a part of its dynamics that cannot be escaped. While evolving, the system may eventually fall into a new and smaller trap domain, which reduces its possible future behaviors (making previous states no longer reachable). This phenomenon depends on biological disruptions or other complex phenomena. Such outline has been interpreted as distinct responses of the organism, such as differentiating into distinct cell types in multicellular organisms [2].
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 longterm behavior. For instance, the model of the cellular development of Drosophila melanogaster, was described using Boolean networks and their attractors [3, 4].
Various kinds of mathematical models have been proposed for the modeling of biological regulatory networks (BRNs). These models include neural networks, differential equations, Petri nets, Boolean networks (BN) as proposed by Kauffman [5], probabilistic Boolean networks, and other multivalued models such synchronous/asynchronous automata networks (AN). In this paper, we use the AN formalism [6, 7] to model BRNs. ANs especially encompass the framework of René Thomas [8].
Qualitative frameworks have received substantial attention, because of their capacity to capture the switching behavior of genetic or biological processes, and therefore, the study of their longterm behavior. This explains our choice of a qualitative representation for the identification of trap domains. In such a qualitative framework, a minimal trap domain can take two different forms: it can be either a steady state, which is one state from which the system does not evolve anymore, called also a fixed point; or an attractor, which is a minimal set of states that loops indefinitely and cannot be escaped.
The computational problem of finding all attractors in a BRN is difficult. Even the simpler problem of deciding whether the system has a fixed point, which can be seen as the smallest kind of attractor, is NPhard [9]. Based on this, many studies have proven that computing attractors in BRNs is also a NPhard problem [10, 11]. Although some methods exist with a lesser complexity, consisting for instance in randomly selecting an initial state and following a long enough trajectory, hoping to eventually finding an attractor, they are not exhaustive and may miss some (hard to reach) attractors.
Therefore, in the absence of more efficient exhaustive methods, it is still relevant to develop an approach to resolve the original NPhard problem of attractors identification. Such an approach consists in exhaustively examine all possible states of a network, along with all possible paths from each of these states. Obviously, this brute force method is very time and memory consuming: \(2^n\) initial states have to be considered for a Boolean model with n nodes; and multivalued networks raise this value even more. Furthermore, a sufficient number of computations have to be performed to ensure that all trajectories have been explored and all attractors are found. This high complexity justifies the use of a tool able to tackle such hard problems.
The simplest way to detect attractors is to enumerate all the possible states and to run simulation from each one until an attractor is reached [12]. This method ensures that all attractors are detected but it has an exponential time complexity, therefore its applicability is highly restricted by the network size.
Regarding BNs only, algorithms for detecting attractors have been extensively studied in the literature. Irons [13] proposes to analyze partial states in order to discard potential attractors more efficiently. This method improves the efficiency from exponential time to polynomial time for a subset of biological Boolean models that is highly dependent on the topology (indegree, outdegree, update functions) of the underlying network. Another method, called GenYsis [14], starts from one (randomly selected) initial state and detects attractors by computing the successor and predecessor states of this initial state. It works well for small BNs, but becomes inefficient for large BNs.
More generally, the efficiency and scalability of attractor detection techniques are further improved with the integration of two techniques. This first is based on binary decision diagrams (BDD), a compact data structure for representing Boolean functions. In a recent work [15], algorithms have been based on the reducedorder BDD (ROBDD) data structure, which further speeds up the computation time of attractor detection. These BDDbased solutions only work for BRNs of a hundred of nodes and also suffer from the infamous state explosion problem, as the size of the BDD depends both on the regulatory functions and the number of nodes in the BRN. The other technique consists in representing the attractor enumeration problem as a satisfiability (SAT) problem such as in [16]. The main idea is inspired by SATbased bounded modelchecking: the transition relation of the BRN is unfolded into a bounded number of steps in order to construct a propositional formula which encodes attractors and which is then solved by a SAT solver. In every step, a new variable is required to represent a state of a node in the BRN. It is clear that the efficiency of these algorithms largely depends on the number of unfolding steps and the number of nodes in the BRN.
In [17], the authors separated the rules that describe the network (the nodes and their interactions: activation or inhibition) from the rules that define its dynamics (for instance: a gene will be activated in the next state if all its activators are active or when at least one of its activators is active at the current state). This allows to obtain more flexible simulations, and the authors also chose to use the declarative paradigm answer set programming (ASP) [18] in order to have more liberty in the expression of evolution rules. They illustrated that specifying large networks with rather complicated behaviors becomes cumbersome and error prone in paradigms like SAT, whereas this is much less the case in a declarative approach such as theirs.
Our goal in this paper is to develop exhaustive methods to analyze a BRN modeled in AN. We address two kinds of issues: finding all possible steady states of a BRN and enumerating all attractors of a given size \(n \ge 2\). We focus on two widespread nondeterministic update schemes (synchronous and asynchronous) and use ASP to solve these aforementioned issues. Although this approach is not new (see above), the use of ASP can still be considered innovative in the field of dynamic properties analysis and our aim here is to assess its computational potential.
Nevertheless, the originality of our contribution is to consider AN models: this formalism does not restrict entities to have Boolean expression levels (active/inactive) as they can have multivalued ones. Complex interactions are modeled in an AN as automata transitions instead of generic influences. This expressiveness allows to represent a wide range of dynamical models with the AN framework, and the particular form of its local transitions can be well handled in ASP. Finally, this framework allows to represent nondeterministic synchronous models, contrary to previous works focusing on asynchronous or deterministic synchronous models.
We previously introduced some rough ideas of this approach in [19]. In the present paper, we have extended this work by focusing on AN models, that are more expressive than the previous process hitting framework [20]. We give a more detailed stateoftheart and a more indepth formalization of the problems tackled and show the merits of our approach on a case study and various benchmarks.
This paper is organized as follows. "Automata networks" presents the main definitions related to the AN and the particular constructs that we will seek: fixed points and attractors. "Answer set programming" briefly presents the ASP framework necessary to understand the encoding part. Section "Fixed points enumeration" details the part of our method that allows to present an AN model using ASP rules and find all the fixed points in such a model. Then, "Length n attractors enumeration" explains how to enumerate all attractors in a model still using ASP. In "Results" we give benchmarks of our methods on several models of different sizes (up to 100 components). Finally, “Conclusion and future direction” concludes and gives some perspectives to this work.
Preliminary definitions
Automata networks
Definition 1 introduces the formalism of automata networks (AN) [6] (see Fig. 1) which allows to model a finite number of discrete levels, called local states, into several automata. A local state is denoted \(a_i\), where a is the name of the automaton, corresponding usually to a biological component, and i is a level identifier within a. At any time, exactly one local state of each automaton is active, modeling the current level of activity or the internal state of the automaton. The set of all active local states is called the global state of the network.
The possible local evolutions inside an automaton are defined by local transitions. A local transition is a triple noted \(a_i \overset{\ell }{\rightarrow } a_j\) and is responsible, inside a given automaton a, for the change of the active local state (\(a_i\)) to another local state (\(a_j\)), conditioned by the presence of a set \(\ell \) of local states belonging to other automata and that must be active in the current global state. Such a local transition is playable if and only if \(a_i\) and all local states in the set \(\ell \) are active. Thus, it can be read as “all the local states in \(\ell \) can cooperate to change the active local state of a by making it switch from \(a_i\) to \(a_j\)”. It is required that \(a_i\) and \(a_j\) are two different local states in automaton a, and that \(\ell \) contains no local state of automaton a. We also note that \(\ell \) should contain at most one local state per automaton, otherwise the local transition is unplayable; \(\ell \) can also be empty.
Definition 1

\(\Sigma = \{a, b,\ldots \}\) is the finite set of automata identifiers;

For each \(a \in \Sigma \), \(\mathcal {S}_a = \{a_i,\ldots ,a_j\}\) is the finite set of local states of automaton a; \(\mathcal {S}= \prod _{a \in \Sigma }\mathcal {S}_a\) is the finite set of global states; \(\user2{LS} = \cup _{{a \in \Sigma }} {\mathcal{S}}_{a} \) denotes the set of all the local states.

For each \(a \in \Sigma \), \(\mathcal {T}_a = \{ a_i \overset{\ell }{\rightarrow } a_j \in \mathcal {S}_a \times \wp (\user2{LS} \setminus \mathcal {S}_a) \times \mathcal {S}_a \mid a_i \ne a_j \}\) is the set of local transitions on automaton a; \(\mathcal {T}= \bigcup _{a \in \Sigma } \mathcal {T}_a\) is the set of all local transitions in the model.
For a given local transition \(\tau = a_i \overset{\ell }{\rightarrow } a_j\), \(a_i\) is called the origin or \(\tau \), \(\ell \) the condition and \(a_j\) the destination, and they are respectively noted \(\mathsf {ori}(\tau )\), \(\mathsf {cond}(\tau )\) and \(\mathsf {dest}(\tau )\).
Example 1

\(\Sigma = \{a, b, c, d\}\),

\(\mathcal {S}_a = \{a_0, a_1\}\), \(\mathcal {S}_b = \{b_0, b_1, b_2\}\), \(\mathcal {S}_c = \{c_0, c_1\}\), \(\mathcal {S}_d = \{d_0, d_1, d_2\} \),

\(\mathcal {T}= \{ \begin{array}[t]{ll} a_0 \overset{\{c_1\}}{\longrightarrow } a_1, a_1 \overset{\{b_2\}}{\longrightarrow } a_0, &{} b_0 \overset{\{d_0\}}{\longrightarrow } b_1, b_0 \overset{\{a_1, c_1\}}{\longrightarrow } b_2, b_1 \overset{\{d_1\}}{\longrightarrow } b_2, b_2 \overset{\{c_0\}}{\longrightarrow } b_0, \\ c_0 \overset{\{a_1, b_0\}}{\longrightarrow } c_1, c_1 \overset{\{d_2\}}{\longrightarrow } c_0, &{} d_0 \overset{\{b_2\}}{\longrightarrow } d_1, d_0 \overset{\{a_0, b_1\}}{\longrightarrow } d_2, d_1 \overset{\{a_1\}}{\longrightarrow } d_0, d_2 \overset{\{c_0\}}{\longrightarrow } d_0 \}\text {.} \end{array}\)
The local transitions given in Definition 1 thus define concurrent interactions between automata. They are used to define the general dynamics of the network, that is, the possible global transitions between global states, according to a given update scheme. In the following, we will only focus on the (purely) asynchronous and (purely) synchronous update schemes, which are the most widespread in the literature. The choice of such an update scheme mainly depends on the considered biological phenomena modeled and the mathematical abstractions chosen by the modeler.
Update schemes and dynamics of automata networks
As explained in the previous section, a global state of an AN is a set of local states of automata, containing exactly one local state of each automaton. In the following, we give some notations related to global states, then we define the global dynamics of an AN.
The active local state of a given automaton \(a \in \Sigma \) in a global state \(\zeta \in \mathcal {S}\) is noted \({\zeta [a]}\). For any given local state \(a_i \in {\mathbf{LS}} \), we also note: \(a_i \in \zeta \) if and only if \({\zeta [a]} = a_i\), which means that the biological component a is in the discrete expression level labeled i within state \(\zeta \). For a given set of local states \(X \subseteq \mathbf {LS} \), we extend this notation to \(X \subseteq \zeta \) if and only if \(\forall a_i \in X, a_i \in \zeta \), meaning that all local states of X are active in \(\zeta \).
Furthermore, for any given local state \(a_i \in \mathbf {LS} \), \(\zeta \Cap a_i\) represents the global state that is identical to \(\zeta \), except for the local state of a which is substituted with \(a_i\): \({(\zeta \Cap a_i)[a]} = a_i \wedge \forall b \in \Sigma{\setminus}\{ a \}, {(\zeta \Cap a_i)[b]} = {\zeta [b]}\). We generalize this notation to a set of local states \(X \subseteq \mathbf {LS} \) containing at most one local state per automaton, that is, \(\forall a \in \Sigma , X \cap \mathcal {S}_a \le 1\) where \(S\) is the number of elements in set S; in this case, \(\zeta \Cap X\) is the global state \(\zeta \) where the local state of each automaton has been replaced by the local state of the same automaton in X, if there exists: \(\forall a \in \Sigma , (X \cap \mathcal {S}_a = \{ a_i \} \Rightarrow {(\zeta \Cap X)[a]} = a_i) \wedge (X \cap \mathcal {S}_a = \emptyset \Rightarrow {(\zeta \Cap X)[a]} = {\zeta [a]})\).
In Definition 2 we formalize the notion of playability of a local transition which was informally presented in the previous section. Playable local transitions are not necessarily used as such, but combined depending on the chosen update scheme, which is the subject of the rest of the section.
Definition 2
(Playable local transitions) Let \(\mathcal {AN}= (\Sigma ,\mathcal {S},\mathcal {T})\) be an automata network and \(\zeta \in \mathcal {S}\) a global state. The set of playable local transitions in \(\zeta \) is called \(P_\zeta \) and defined by: \(P_\zeta = \{ a_i \overset{\ell }{\rightarrow } a_j \in \mathcal {T}\mid \ell \subseteq \zeta \wedge a_i \in \zeta \}\).
The dynamics of the AN is a composition of global transitions between global states, that consist in applying a set of local transitions. Such sets are different depending on the chosen update scheme. In the following, we give the definition of the asynchronous and synchronous update schemes by characterizing the sets of local transitions that can be “played” as global transitions. The asynchronous update sets (Definition 3) are made of exactly one playable local transition; thus, a global asynchronous transition changes the local state of exactly one automaton. On the other hand, the synchronous update sets (Definition 4) consist of exactly one playable local transition for each automaton (except the automata where no local transition is playable); in other words, a global synchronous transition changes the local state of all automata that can evolve at a time. Empty update sets are not allowed for both update schemes. In the definitions below, we willingly mix the notions of “update scheme” and “update set”, which are equivalent here.
Definition 3
Definition 4
Once an update scheme has been chosen, it is possible to compute the corresponding dynamics of a given AN. Thus, in the following, when it is not ambiguous and when results apply to both of them, we will denote by \(U^{}\) a chosen update scheme among \(U^{\mathsf {asyn}}\) and \(U^{\mathsf {syn}}\). Definition 5 formalizes the notion of a global transition depending on a chosen update scheme \(U^{}\).
Definition 5
We note that in a deterministic dynamics, each state has only one successor. However, in case of nondeterministic dynamics, such as the asynchronous and synchronous update schemes of this paper, each state may have several possible successors.
Example 2
Figures 2 and 3 illustrate respectively the asynchronous and synchronous update schemes on the model of Fig. 1. Each global transition is depicted by an arrow between two global states. Only an interesting subset of the whole dynamics is depicted in both figures.
At this point, it is important to remind that the empty set never belongs to the update schemes defined above: \(\forall \zeta \in \mathcal {S}, \emptyset \notin U^{\mathsf {asyn}}(\zeta ) \wedge \emptyset \notin U^{\mathsf {syn}}(\zeta )\). The consequence on the dynamics is that a global state can never be its own successor. In other words, even when no local transition can be played in a given global state (i.e., \(P_\zeta = \emptyset \)), we do not add a “selftransition” on this state. Instead, this state has no successors and is called a fixed point, as defined later in this section.
Definition 6 explains what are inconflict local transitions, which are interesting in the scope of the synchronous update scheme. Two local transitions are inconflict if they belong to the same automaton and produce some nondeterminism inside this automaton. Such phenomenon arises when both local transitions have the same origin and compatible conditions, but their destinations are different; or, in other words, there exists a global state in which they are both playable. In such a case, they allow the automaton to evolve in two different possible local states from the same active local state, thus producing a nondeterministic behavior. This definition will be used in the discussion of the next section and in "Length n attractors enumeration".
Definition 6
Finally, Definition 7 introduces the notions of path and trace which are used to characterize a set of successive global states with respect to a global transition relation. Paths are useful for the characterization of attractors that are the topic of this work. The trace is the set of all global states traversed by a given path (thus disregarding the order in which they are visited). We note that a path is a sequence and a trace is a set.
Definition 7
(Path and trace) Let \(\mathcal {AN}= (\Sigma ,\mathcal {S},\mathcal {T})\) be an automata network, \(U^{}\) an update scheme and \(n \in \mathbb {N}\setminus \{ 0 \}\) a strictly positive integer. A sequence \(H = ( H_i )_{i \in \llbracket 0; n \rrbracket } \in \mathcal {S}^{n+1}\) of global states is a path of length n if and only if: \(\forall i \in \llbracket 0; n1 \rrbracket , H_i \rightarrow _{U^{}} H_{i+1}\). H is said to start from a given global state \(\zeta \in \mathcal {S}\) if and only if: \(H_0 = \zeta \). Finally, the trace related to such a path is the set of the global states that have been visited: \(\mathsf {trace}(H) = \{ H_j \in \mathcal {S}\mid j \in \llbracket 0; n \rrbracket \}\).
In the following, when we define a path H of length n, we use the notation \(H_i\) to denote the ith element in the sequence H, with \(i \in \llbracket 0; n \rrbracket \). We also use the notation \(H = n\) to denote the length of a path H, allowing to write: \(H_{H}\) to refer to its last element. We also recall that a path of length n models the succession of n global transitions, and thus features up to n + 1 states (some states may be visited more than once).
Example 3
When there is one or several repetitions in a given path of length n (i.e., if a state is visited more than once), its trace is then of size strictly lesser than n + 1. More precisely, one can compute the size of the trace corresponding to a given path by subtracting the number of repetitions in that path (Lemma 1). For this, we formalize in Definition 8 the notion of repetitions in a path, that is, the global states that are featured several times, designated by their indexes.
Definition 8
Lemma 1
Proof of Lemma 1
By definition of a set, and knowing that \(\mathsf {sr}(H)\) counts the number of states that exist elsewhere in H with a lesser index. \(\square \)
We note that if there is no repetition in a path of length n (\(\mathsf {sr}(H)=\emptyset \Rightarrow \mathsf {sr}(H)=0\)), then the number of visited states is exactly: \(\mathsf {trace}{(H)} = n+1\).
Example 4
We can check Lemma 1 on the path H given in Example 3. Indeed, \(\langle a_1, b_2, c_1, d_1 \rangle \) is featured 3 times at \(H_0\), \(H_2\) and \(H_6\). Then, according to the Definition 8, this state is repeated twice at \(H_2\) and \(H_6\) because the first visit of this state is not computed in \(\mathsf {sr}(H)\). In addition, the state \(\langle a_0, b_2, c_1, d_1 \rangle \) is featured twice in this path, at \(H_1\) and \(H_5\), therefore it is considered as repeated once at \(H_5\). Thus, \(\mathsf {sr}(H)=\{2,6,5\}\), \(\mathsf {sr}(H)=3\) and \(\mathsf {trace}(H) = 6 + 1  3 = 4\).
Determinism and nondeterminism of the update schemes
In the general case, in multivalued networks, both the asynchronous and synchronous update schemes are nondeterministic, which means that a global state can have several successors.
In the case of the asynchronous update scheme, the nondeterminism may come from inconflict local transitions, but it actually mainly comes from the fact that exactly one local transition is taken into account for each global transition (see Definition 3). Thus, for a given state \(\zeta \in \mathcal {S}\), as soon as \(P_\zeta  > 1\), several successors may exist. In the model of Fig. 1, for example, the global state \(\langle a_1, b_2, c_0, d_1 \rangle \) (in green on Fig. 2) has three successors: \(\langle a_1, b_2, c_0, d_1 \rangle \rightarrow _{U^{\mathsf {asyn}}} \langle a_0, b_2, c_0, d_1 \rangle \), \(\langle a_1, b_2, c_0, d_1 \rangle \rightarrow _{U^{\mathsf {asyn}}} \langle a_1, b_0, c_0, d_1 \rangle \) and \(\langle a_1, b_2, c_0, d_1 \rangle \rightarrow _{U^{\mathsf {asyn}}} \langle a_1, b_2, c_0, d_0 \rangle \).
In the case of the synchronous update scheme (see Definition 4), however, the nondeterminism on the global scale is only generated by inconflict local transitions (see Definition 6), that is, local transitions that create nondeterminism inside an automaton. For example, the model of Fig. 1 features two local transitions \(b_0 \overset{\{d_0\}}{\longrightarrow } b_1\) and \(b_0 \overset{\{a_1, c_1\}}{\longrightarrow } b_2\) that can produce the two following global transitions from the same state (depicted by red arrows on Fig. 3): \(\langle a_1, b_0, c_1, d_0 \rangle \rightarrow _{U^{\mathsf {syn}}} \langle a_1, b_1, c_1, d_0 \rangle \) and \(\langle a_1, b_0, c_1, d_0 \rangle \rightarrow _{U^{\mathsf {syn}}} \langle a_1, b_2, c_1, d_0 \rangle \). Note that for this particular case, these transitions also exist for the asynchronous scheme (also depicted by red arrows on Fig. 2).
Therefore, it is noteworthy that if every automaton contains only two local states (such a network is often called “Boolean”) then the synchronous update scheme becomes completely deterministic. Indeed, it is not possible to find inconflict local transitions anymore because for each possible origin of a local transition, there can be only one destination (due to the fact that the origin and destination of a local transition must be different). This observation can speed up the computations in this particular case.
Fixed points and attractors in automata networks
Studying the dynamics of biological networks was the focus of many works, explaining the diversity of existing frameworks dedicated to modeling and the different methods developed in order to identify some patterns, such as attractors [9, 11, 17, 21, 22]. In this paper we focus on several subproblems related to this: we seek to identify the steady states and the attractors of a given network. The steady states and the attractors are the two longterm structures in which any dynamics eventually falls into. Indeed, they consist in terminal (sets of) global states that cannot be escaped, and in which the dynamics always ends.
In the following, we consider a BRN modeled in AN \((\Sigma ,\mathcal {S},\mathcal {T})\), and we formally define these dynamical properties. We note that since the AN formalism encompasses Thomas modeling [8], all our results can be applied to the models described by this formalism, as well as any other framework that can be described in AN (such as Boolean networks, Biocham [23]...).
A fixed point is a global state which has no successor, as given in Definition 9. Such global states have a particular interest as they denote conditions in which the model stays indefinitely. The existence of several of these states denotes a multistability, and possible bifurcations in the dynamics [1].
Definition 9
It is notable that the set of fixed points of a model (that is, the set of states with no successor) is the same in both update schemes asynchronous and synchronous update [24, 25]: \(\forall \zeta \in \mathcal {S}, U^{\mathsf {asyn}}(\zeta ) = \emptyset \Longleftrightarrow U^{\mathsf {syn}}(\zeta ) = \emptyset .\)
Example 5
The statetransition graphs of Figs. 2 and 3 depict three fixed points colored in red: \(\langle a_1, b_1, c_1, d_0 \rangle \), \(\langle a_1, b_1, c_0, d_0 \rangle \) and \(\langle a_0, b_0, c_0, d_1 \rangle \). Visually, they can be easily recognized because they have no outgoing arrow (meaning that they have no successors). Although these figures do not represent the whole dynamics, but they allow to check that in both update schemes the fixed points are the same, at least on this subset of the overall behavior.
Another complementary dynamical pattern consists in the notion of nonunitary trap domain (Definition 10), which is a (nonsingleton) set of states that the dynamics cannot escape, and thus in which the system indefinitely remains. In this work, we focus more precisely on (nonsingleton) attractors (Definition 11), that are cyclic and minimal trap domains in terms of sets inclusion. In order to characterize such attractors, we use the notion of cycle (Definition 12), which is a looping path. Indeed, a cycle is a strongly connected component (Lemma 2), allowing us to give an alternative definition for an attractor (Lemma 3). Formally speaking, fixed points can be considered as attractors of size 1. However, in the scope of this paper and for the sake of clarity, we call “attractors” only nonunitary attractors, that is, only sets containing at least two states. This is justified by the very different approaches developed for fixed points and attractors in the next sections.
Definition 10
Definition 11
(Attractor) Let \(\mathcal {AN}= (\Sigma ,\mathcal {S},\mathcal {T})\) be an automata network and \(U^{}\) an update scheme. A set of global states \(\mathbf {A}\), with \(\mathbf {A} \ge 2\), is called an attractor (regarding scheme \(U^{}\)) if and only if it is a minimal trap domain in terms of inclusion.
Definition 12
Example 6
The path H of length 6 given in Example 3 is a cycle because \(H_0 = H_6\).
Lemma 2 states that the set of (traces of) cycles in a model is exactly the set of strongly connected components. Indeed, a cycle allows to “loop” between all states that it contains, and conversely, a cycle can be built from the states of any strongly connected component. This equivalence is used in the next lemma.
Lemma 2
(The traces of cycles are the SCCs) The traces of the cycles are exactly the strongly connected components (with respect to the global transition relation).
Proof of Lemma 2
(\(\Rightarrow \)) From any state of a cycle, it is possible to reach all the other states (by possibly cycling). Therefore, the trace of this cycle is a strongly connected component. (\(\Leftarrow \)) Let \(\mathbf {S} = \{ \zeta _{i} \}_{i \in \llbracket 0; n \rrbracket }\) be a strongly connected component in which the elements are arbitrarily labeled. Because it is a strongly connected component, for all \(i \in \llbracket 0; n \rrbracket \), there exists a path \(H^i\) made of elements of \(\mathbf {S}\) so that \(H^i_0 = \zeta _i\) and \(H^i_{H^i} = \zeta _{i+1}\) (or \(H^n_{H^n} = \zeta _0\) for \(i = n\)). We create a path \(\mathbf {C}\) by concatenation of all paths \(H^0, H^1, \ldots , H^n\) by merging the first and last element of each successive path, which is identical: \(\forall i \in \llbracket 0; n1 \rrbracket , H^i_{H^i} = \zeta _{i+1} = H^{i+1}_0\). \(\mathbf {C}\) is a cycle, because \(\mathbf {C}_0 = H^0_0 = \zeta _0 = H^n_{H^n} = \mathbf {C}_{\mathbf {C}}\). Furthermore, \(\forall i \in \llbracket 0; n \rrbracket , \zeta _i = H^i_0 \in \mathsf {trace}(\mathbf {C})\), thus \(\mathbf {S} \subseteq \mathsf {trace}(\mathbf {C})\). Finally, only states from \(\mathbf {S}\) have been used to build \(\mathbf {C}\), thus \(\mathsf {trace}(\mathbf {C}) \subseteq \mathbf {S}\). Therefore, \(\mathsf {trace}(\mathbf {C}) = \mathbf {S}\). \(\square \)
In Definition 11, attractors are characterized in the classical way, that is, as minimal trap domains. However, we use an alternative characterization of attractors in this paper, due to the specifics of ASP: Lemma 3 states that an attractor can alternatively be defined as a trap domain that is also a cycle, and conversely. In other words, the minimality requirement is replaced by a cyclical requirement.
Lemma 3
(The attractors are the trap cycles) The attractors are exactly the traces of cycles which are trap domains.
Proof of Lemma 3
(\(\Rightarrow \)) By definition, an attractor is a trap domain. It is also a strongly connected component, and thus, from Lemma 2, it is the trace of a cycle. (\(\Leftarrow \)) Let \(\mathbf {C}\) be both a cycle and a trap domain. From Lemma 2, \(\mathbf {C}\) is also a strongly connected component. Let us prove by contradiction that \(\mathbf {C}\) is a minimal trap domain, by assuming that it is not minimal. This means that there exists a smaller trap domain \(\mathbf {D} \subsetneq \mathbf {C}\). Let us consider \(x \in \mathbf {D}\) and \(y \in \mathbf {C} \setminus \mathbf {D}\). Because \(\mathbf {D}\) is a trap domain, it exists no path between x and y; this is in contradiction with \(\mathbf {C}\) being a strongly connected component (as both x and y belong to \(\mathbf {C}\)). Therefore, \(\mathbf {C}\) is a minimal trap domain, and thus an attractor. \(\square \)
As explained before, Lemma 3 will beused in "Length n attractors enumeration". Indeed, directly searching for minimal trap domains would be too cumbersome; instead, we enumerate cycles of length n in the dynamics of the model and filter out those that are not trap domains. The remaining results are the attractors formed of cycles of length n. The previous lemma ensures the soundness and completeness of this search for a given value of n.
Lemma 4
(Characterization of nonattractors) Let \(\mathbf {A} \subset \mathcal {S}\) be a set of states. If \(\exists \zeta _1 \in \mathbf {A}\) and \(\exists \zeta _2 \in \mathcal {S}\setminus \mathbf {A}\) such that \(\zeta _1 \rightarrow _{U^{}} \zeta _2\) then \(\mathbf {A}\) is not an attractor.
Proof of Lemma 4
By definition, \(\mathbf {A}\) is not a trap domain (Definition 10) and thus it is not an attractor (Definition 11). \(\square \)
Example 7

\(\{ \langle a_0, b_1, c_0, d_0 \rangle , \langle a_0, b_1, c_0, d_2 \rangle \}\) is depicted in blue and appears in both figures. It is a cyclic attractor, because it contains exactly one cycle.

\(\{ \langle a_0, b_2, c_1, d_0 \rangle , \langle a_0, b_2, c_1, d_1 \rangle , \langle a_1, b_2, c_1, d_1 \rangle , \langle a_1, b_2, c_1, d_0 \rangle \}\) is only present for the asynchronous update scheme and is depicted in yellow on Fig. 2. It is a complex attractor, that is, a composition of several cycles.

\(\{ \langle a_1, b_2, c_1, d_1 \rangle , \langle a_0, b_2, c_1, d_0 \rangle \}\) is, on the contrary, only present for the synchronous update scheme and is depicted in gray on Fig. 3. It is also a cyclic attractor.
The aim of the rest of this paper is to tackle the enumeration of fixed points ("Fixed points enumeration") and attractors ("Length n attractors enumeration") in an AN. For this, we use ASP ("Answer set programming") which is a declarative paradigm dedicated to the resolution of complex problems.
Answer set programming
For the present work, we used Clingo ^{1} [27] which is a combination of a grounder and a solver. In the rest of this paper, we use ASP to tackle the problems of enumerating all fixed points and attractors of a given AN model.
Fixed points enumeration
The first aspect of our work is the enumeration of a special type of trap domains: fixed points (also called stable states or steady states) which are composed of only one global state (see Definition 9). They can be studied separately from attractors because their enumeration follows a different pattern which is more specific to this problem. A previous version of this work using another framework (process hitting) is presented in [19]. Although the main idea is preserved, the work we present here is different because we are interested in the more expressive AN framework in which the transitions have a different form.
Translating automata networks into answer set programs
Before any analysis of an AN, we first need to express it with ASP rules. We developed a dedicated converter named AN2ASP ^{2} and we detail its principle in the following.
First, the predicate automatonLevel(A,I) is used to define each automaton A along with its local state I. The local transitions are then represented with two predicates: condition which defines each element of the condition along with the origin, and target which defines the target of the local transition. Each local transition is labeled by an identifier that is the same in its condition and target predicates. Example 8 shows how an AN model is defined with these predicates.
Example 8
Besides, all the local transitions of the network are defined in lines 7–21; for instance, all the predicates in line 7 declare the transition \(\tau _1 = a_0 \overset{\{c_1\}}{\longrightarrow } a_1\), which is labeled 1. We declare as many predicates condition as necessary in order to fully define a local transition \(\tau \) that has potentially several elements in its condition \(\mathsf {cond}(\tau )\). For instance, transition \(b_0 \overset{\{a_1, c_1\}}{\longrightarrow } b_2\) is defined in line 11 with label 4 and requires three of these predicates for \(b_0\), \(a_1\) and \(c_1\). Finally, in lines 4–5, predicate automaton gathers all existing automata names in the model, and predicate localTrans gathers all transition labels. The underscore symbol (_) in the parameters of a predicate is a placeholder for any value.
Since the names of the biological components may start with a capital letter, it is preferable to use the double quotes (“”) around the automata names in the parameters of all predicates to ensure that the automata names are understood as constants by the ASP grounder and not as variables.
Fixed points search
Line 29 constitutes a cardinality rule (as defined in "Answer set programming") whose consequence is the enumeration of all global states of the model in distinct answer sets. Each global state is defined by considering exactly one local state for each existing automaton from the shown ones defined in shownAutomatonLevel. Each global state is described using predicates fix(A,I), named in anticipation of the final fixed point results, where I is the active local state of automaton A.
The last step consists in filtering out any global state \(\zeta \), that is not a fixed point, among all generated states. In this case, it consists in eliminating all candidate answer sets in which at least one local transition can be played, that is, where \(P_\zeta \ne \emptyset \). Such a filtering part is ideally realized with the use of one or several constraints. As explained in "Answer set programming", a constraint removes all answer sets that satisfy its righthand part. Regarding our problem, an answer set representing a given global state must be filtered out if there exists at least one playable local transition in this state (line 33). A transition T is considered as unplayable in a state, that is, \(\texttt {T} \notin P_\zeta \), if at least one of its conditions is not satisfied. For this, predicate unPlayable(T) defined in line 31, flags a local transition as unplayable when one of its condition contains a local state that is different from the local state of the same automaton. This is used in the final constraint (line 33) which states that if there exists a local transition which is playable in the considered global state (i.e., \(\exists \texttt {T} \in \mathcal {T}, \texttt {T} \in P_\zeta \)) then this global state should be eliminated from the result answer sets (because it is not a fixed point). In the end, the fixed points of a considered model are exactly the global states represented in each remaining answer sets, described by the set of the atoms fix(A,I) which define each automaton local state.
Example 9
(Fixed point enumeration) The AN model of Fig. 1 contains 4 automata: a and c have 2 local states while b and d have 3; therefore, the whole model has \(2*2*3*3 = 36\) states (whether they can be reached or not from a given initial state). We can check that this model contains exactly 3 fixed points: \(\langle a_1, b_1, c_0, d_0 \rangle \), \(\langle a_1, b_1, c_1, d_0 \rangle \) and \(\langle a_0, b_0, c_0, d_1 \rangle \). All of them are represented in both Figs. 2 and 3. In this model, no other state verifies this property. We recall that the fixed points are identical for the synchronous and asynchronous update schemes [24].
Length n attractors enumeration
In the previous section we gave a method to enumerate all fixed points of a given model. In a sense, a fixed point can be considered as an attractor: it cannot be escaped and its size (\(n=1\)) makes it trivially minimal. However, attractors in the general case are made of several states. In the rest of this paper, we exclude onestate attractors (tackled in the last section "Fixed points enumeration"). We focus on attractors composed with severalstates (following Definition 11) and we describe how to obtain some or all the attractors of a given length in a model. Obtaining all attractors of any length can be theoretically tackled by gradually increasing the considered length.
 1.
Enumerate all paths of length n,
 2.
Remove all paths that are not cycles,
 3.
Remove all cycles that are not trap domains (i.e., keep only attractors).
Cycles enumeration
The approach presented here first enumerates all the paths of length n in the AN model (Definition 7).
In an ASP program, it is possible to instantiate constants whose values are defined by the user at each execution: this is the role of the lowercase n in step(0..n) (line 26), that represents the number of considered steps. For example, knowing the initial global state, step(0..5) will compute all paths of length 5 (thus containing 6 successive global states).
In our approach, we tackle separately the computation of the dynamics and the resolution of our problem (namely, attractors enumeration). We show in the following how to compute the evolution of the model through the asynchronous and the synchronous update schemes, as presented in "Update schemes and dynamics of automata networks". The piece of program that computes the attractors, given afterwards, is common to whatever update schemes.
All possible evolutions of the network (that is, the resulting paths after playing a set of global transitions) can be enumerated with a cardinality rule (explained in "Answer set programming") such as the ones in line 33 for the asynchronous update scheme, and line 39 for the synchronous update scheme. Such rules reproduce all possible paths in the dynamics of the model by representing each possible successor of a considered state as an answer set. This enumeration encompasses the nondeterministic behavior (in both update schemes).
In a nutshell, one should choose one of both pieces of program presented above, that is, either lines 39–41 for the asynchronous update scheme, or lines 39–41 for the synchronous one. The overall result of both of these pieces of programs is a collection of answer sets, where each answer is a possible path of length n (that is, computed in n steps) and starting from any initial state (at step 0).
Between two consecutive steps S and S+1, we witness that the active level of a given automaton B has changed with the predicate change in line 43, which stores the chosen local transition.
As given in Lemma 2, all remaining paths are strongly connected components. We finally need to verify if they are trap domains (Lemma 3) in order to discriminate attractors.
Attractors enumeration
Due to the nondeterministic behavior in the dynamics, each state in the statetransition graph of a given AN may have several successors. Therefore a cyclic path is not necessarily an attractor. The only certain exception is the case of the deterministic synchronous update scheme (such as in Boolean models, as explained in Section "Determinism and nondeterminism of the update schemes"). In this case, the computation may be stopped here because a cycle is necessarily an attractor. This result is used in [28–30].
In the rest of this section, we will tackle the more general and challenging case of nondeterminism. Indeed, in the general case, some local transitions may allow the dynamics to escape the cycle; in such case, the cycle would not even be a trap domain (see Lemma 4). For instance, in the partial statetransition graph of Fig. 2, we can spot many cycles of various lengths but not all of them are attractors. In particular, the initial global state is part of a cycle of length 2 which is not an attractor, and which trace is: \(\{ \langle a_1, b_2, c_0, d_1 \rangle , \langle a_1, b_2, c_0, d_0 \rangle \}\).
That is why another check is required to filter out all the remaining cycles that can be escaped (and are therefore not attractors). Once again, this filtering is performed with constraints, which are the most suitable solution. In order to define such constraints, we need to describe the behavior that we do not wish to observe: escaping the computed cycle. For this, it is necessary to differentiate between the effectively played local transitions (played) and the ”also playable” local transitions which were not played (alsoPlayable in line 61). Then, we verify at each step S, comprised between 0 and n, if these also playable local transitions make the system evolve or not to a new global state that is not a part of the cycle trace.
Example 10
Complex attractors
Up to this point, we managed to propose an ASP program that enumerates all the attractors in a given AN. Each attractor is the trace of a path of length n. In many cases, except for some complex attractors, this length n (which corresponds to the number of played global transitions in the path) is also equal to the number of visited states (i.e., the size of the trace). This is a trivial case of a minimal path covering a given attractor, that is, no path of lesser length can cover it. Indeed, as in the examples of attractors in Figs. 2 and 3, enumerating the paths of length 2 is enough to obtain all the attractors having two global states, and the same goes for the attractors of length 4. But without the constraint that we develop below (given in lines 70–93), when the program is asked to display the attractors covered by a path of length n, it will also return various paths of size lower than n by considering nonminimal paths, that is, containing unwanted repetitions inside the cycle, or even repetitions of the entire cycle. In the example of Fig. 3, for instance, with \(\texttt {n}=6\), the program returns the two attractors, although they both are of size 2. Indeed, each of them can be covered by a cycle of length 6: it consists of a cycle of size 2 repeated three times.
Therefore, the objective of this section is to exclude most cases where a cycle is nonminimal, such as the obvious one where it is entirely repeated, because such a case is useless with respect to the computation of attractors. Moreover, we would prefer that our method yields no answer set when no attractor traversed by a cycle of length n is found (even if nonminimal attractors on cycles of lesser length are found). We don’t formally claim here that our method eliminates all of these cases, but we aim at tackling most of these cases in order to sanitize the answer set as much as possible. For instance, an attractor \(\zeta _0 \rightarrow \zeta _1 \rightarrow \zeta _0\) of length \(\texttt {n}=2\) could be listed among the attractors of length \(\texttt {n}=4\) if it is repeated twice as the following path: \(\zeta _0 \rightarrow \zeta _1 \rightarrow \zeta _0 \rightarrow \zeta _1 \rightarrow \zeta _0\). Although all general assumptions regarding attractors are verified (it consists in a cycle and all the global transitions produce global states that are still cycle), we aim at willingly excluding it from the answers because it is not minimal in terms of length.
Our objective in the following is to propose a program that returns, as far as possible, all attractors of the model that actually correspond to a path of length n which is minimal. We propose the following rules to verify this property; each of them concludes with the atom isNotMinimal(n), which means that the considered cycle is not minimal. In the end, isNotMinimal(n) is used in the constraint of line 93 which eliminates all these unwanted cases together.
We first verify if there exists a path of length X < n without repetitions from the state of step 0 to step X, where X is the trace size of the cycle, that is, the number of different states in the path. Then we also verify if there is a transition from the state of step X to the state of step 0. If both properties are true, then there exists a path of size X < n that covers all the states of the attractor, and thus n is not the minimal path length of that attractor (lines 81–84).
Nevertheless, there is still a case of duplicate results that cannot be tackled by the previous constraint: the existence of several minimal cycles for the same attractor. Indeed, for one given attractor, it is possible to find several minimal covering cycles by changing the initial state, or the traversal (in the case of complex attractors). For instance, the hypothetical attractor \(\{ \zeta _0 ; \zeta _1 \}\) is captured by the two cycles: \(\zeta _0 \rightarrow \zeta _1 \rightarrow \zeta _0\) and \(\zeta _1 \rightarrow \zeta _0 \rightarrow \zeta _1\). This leads to repetitions which are not removed from the answers of our method.
The final result presented by each answer set is described by the collection of atoms active(ALs,S), where S denotes the label of one of the steps in the cycle, and ALs corresponds to one of the active local states.
The problem of finding attractors in a discrete network is NPhard, therefore the implementation that we gave in this section also faces such a complexity. However, ASP solvers (namely, Clingo in our case) are specialized in tackling such complex problems. Next section will be dedicated to the results of several computational experiments that we performed on biological networks. We show that our ASP implementation can return results in only a few seconds attractors of small size even on models with 100 components, which is considered large.
Results
In this section, we exhibit several experiments conducted on biological networks. We first detail the results of our programs on the AN model of Fig. 1. Then, we sum up the results of benchmarks performed on other models up to 100 components. In general, the time performances are good and the overall results confirm the applicability ASP for the verification of formal properties or the enumeration of special constructs on biological systems.
All experiments were run on a desktop PC with a Pentium VII 3 GHz processor and 16 GB memory.
Case study
We first conducted detailed experiments on the 4components model of Fig. 1. As detailed in "Automata networks", this network contains 4 automata and 12 local transitions. Its statetransition graph comprises 36 different global states and some of them are detailed in the partial statetransition graphs in Fig. 2 (for the asynchronous update scheme) and Fig. 3 (for the synchronous update scheme).

Asynchronous update scheme:

\(\texttt {n}=0\): \(\langle a_1, b_1, c_1, d_0 \rangle \),\(\langle a_1, b_1, c_0, d_0 \rangle \) and \(\langle a_0, b_0, c_0, d_1 \rangle \);

\(\texttt {n}=2\): \(\{ \langle a_0, b_1, c_0, d_0 \rangle , \langle a_0, b_1, c_0, d_2 \rangle \}\);

\(\texttt {n}=4\): \(\{ \langle a_1, b_2, c_1, d_1 \rangle , \langle a_0, b_2, c_1, d_1 \rangle , \langle a_0, b_2, c_1, d_0 \rangle , \langle a_1, b_2, c_1, d_0 \rangle \}\).


Synchronous update scheme:

\(\texttt {n}=0\): \(\langle a_1, b_1, c_1, d_0 \rangle \),\(\langle a_1, b_1, c_0, d_0 \rangle \) and \(\langle a_0, b_0, c_0, d_1 \rangle \);

\(\texttt {n}=2\): \(\{\langle a_0, b_1, c_0, d_0 \rangle , \langle a_0, b_1, c_0, d_2 \rangle \}\) and \(\{ \langle a_1, b_2, c_1, d_1 \rangle , \langle a_0, b_2, c_1, d_0 \rangle \}\).

Moreover, executing the programs with \(\texttt {n}\ne 2\) and \(\texttt {n}\ne 4\) returns no results, which means that the solver correctly terminates having found no answer set. This is expected because there is no attractor of length different than 2 and 4 for this model, and we excluded repeated cycles from the results (therefore, the attractors already found for lengths 2 and 4 are not found for \(\texttt {n}=6\) or \(\texttt {n}=8\), for instance). For this small network, all the results are computed in less than 0.05 second.
Benchmarks

Lambda phage: a regulatory network featuring some viral genes crucial in the decision between lysis and lysogenization in temperate bacteriophage lambda [31];

Trpreg: a qualitative model of regulated metabolic pathways of the tryptophan biosynthesis in E. coli [32];

Fissionyeast: a cell cycle model of Schizosaccharomyces pombe [33];

Mamm: a mammalian cell cycle model [34];

Tcrsig: a signaling and regulatory network of the TCR signaling pathway in the mammalian differentiation [35];

FGF: a drosophila FGF signaling pathway [36];

Thelper: a model of the Thelper cells differentiation and plasticity, which accounts for novel cellular subtypes [37].

The existing GINsim tool allows to export its models into the SBML qual formalism;

The existing LogicalModel library [39, 40] can convert SBML qual models into AN models;

Finally, our script AN2ASP converts AN models into ASP programs, following the principles detailed in "Translating automata networks into answer set programs".
Brief description of the models used in our benchmarks
Models  Model description  

\(\Sigma \)  \(\max _{a\in \Sigma }\{\mathcal {S}_a\}\)  \(\mathcal {T}\)  \(\mathcal {S}\)  
Example  4  3  12  36  
[31]  Lambda phage  4  4  46  48 
[32]  Trpreg  4  3  14  36 
[33]  Fissionyeast  9  3  43  3 \(\times \) 2\(^\text {9}\) \(=\) 1536 
[34]  Mamm.  10  2  34  2\(^{10}\) \(=\) 1024 
[35]  Tcrsig  40  2  85  2\(^\text {40}\) \(\simeq \) 10\(^{12}\) 
[36]  FGF  59  3  102  2\(^\text {31}\) \(\simeq \) 1.2 \(\times \) 10\(^\text {10}\) 
[37]  Thelper  101  3  316  2\(^\text {102}\) \(\simeq \) 5.7 \(\times \) 10\(^\text {31}\) 
Results of our fixed points enumeration implementation
Models  Fixed points enumeration for both update schemes  

\(\Delta ^{\mathsf {all}} t\) (ms)  #\(^{\mathsf {all}} \mathbf {F}\)  
Example  2  3 
Lambda phage  4  1 
Trpreg  6  2 
Fissionyeast  5  1 
Mamm.  3  1 
Tcrsig  5  8 
FGF  25  1536 
Thelper  170,642  5,875,504 
We note that all the results for the fixed points search have been compared and confirmed using GINsim [38] and Pint [39]. Regarding the attractor enumeration, we compared our results with Boolean network system (BNS) [16] for the synchronous update scheme on the Fissionyeast, Mamm., and Tcrsig models; and with GINsim [38] for the asynchronous update scheme on the Lambda phage, Trpreg, Fissionyeast and Mamm. models. In all cases, we found the same results. It is interesting to note that our method allows to return a response regarding attractors of small size even on big models. In contrast, other tools may take a very long time or even fail to answer. For instance, that happens with GINsim for the Tcrsig, FGF and Thelper models. Indeed, they are based on the computation of the complete transition graph even for the study of small attractors.
Our results could not be compared with, for example, the existing ASPG method [17]. Indeed, with this tool, the user has to choose an update rule on which the dynamic evolution will be based on. For instance, one rule consists in activating a gene when at least one of its activators is active while no inhibitor is; another one activates a gene when it has more expressed activators than inhibitors. Because the chosen activation rule is applied for all the components of the model, while the evolution rules in our AN semantics are specific to each component, the results of both tools cannot be strictly compared.
We recall that among the results output, some attractors may be listed several times in the answers, despite any filtering, as explained at the end of "Complex attractors". Indeed, the solver returns different answer sets for different paths that cover the same trace but differ in terms of initial global state. Therefore, in the results of Table 3, we focused on the conclusion and computation times of the search of any first found attractor of length n.
In case the user may need the exhaustive list of all attractors, our method can also list all the answers, including these repetitions. For instance, our method yields 4 answers for the Trpreg model and a cycle length of \(\texttt {n} = 4\) with the asynchronous update scheme, and the computation takes 47 ms; this typically represents an attractor of size 4 where each answer set represents a cycle starting from a different initial state. Regarding the Thelper model (the largest studied model with 101 automata), the search for all attractors of size \(\texttt {n} = 2\) with the synchronous update scheme takes about 275 s (\(\sim \)5 min) and returns 2,058,272 answers, while it takes only 57 s to return all the attractors of size n=12, (6144 answers). However, as explained before, these results mean that this model features strictly less than, for instance, 6144 attractors covered by a cycle of length 12, because each one is repeated several times.
Results of our attractors enumeration implementation
Models  n  Attractors enumeration  

Asynchronous scheme  Synchronous scheme  
\(\Delta t\) (ms)  \(\exists ?\mathbf {A}\)  \(\Delta t\) (ms)  \(\exists ?\mathbf {A}\)  
Example  2  7  Yes  7  Yes 
4  16  Yes  14  No  
8  98  No  75  No  
Lambda phage  2  14  Yes  14  Yes 
10  1352  No  842  No  
20  15,656  No  14,452  No  
Trpreg  2  8  No  7  No 
4  14  Yes  15  No  
20  3908  No  3808  No  
Fission yeast  2  16  No  16  Yes 
10  1011  No  807  No  
20  17,302  No  16,313  No  
Mamm.  2  12  No  12  No 
7  177  No  147  Yes  
10  720  No  605  No  
20  58,133  No  9253  No  
Tcrsig  2  26  No  25  No 
6  353  No  288  Yes  
10  2420  No  1841  No  
20  85,599  No  27,078  No  
FGF  2  38  No  36  No 
10  2080  No  1953  No  
20  30,861  No  29,838  No  
Thelper  2  180  No  125  Yes 
3  391  No  301  Yes  
4  782  No  1064  No  
6  4271  No  2372  Yes  
7  7909  No  3522  Yes  
9  26,443  No  7042  Yes  
10  44,924  No  12,208  Yes  
12  107,358  No  28,520  Yes  
20  4,230,836 \(\sim \) 1h17  No  187,105 \(\sim \) 3min  No 
Conclusion and future direction
In this paper, we presented a new logical approach to efficiently compute the list of all fixed points and attractors in biological regulatory networks. We formalized our approach using the AN framework, which is bisimilar to many logical networks [41]. All results given here can thus be applied to the widespread Thomas’ modeling [42] in the asynchronous scheme and to the Kauffman modeling in the synchronous scheme [43]. In addition, this framework can encompass any update rules, such as the ones represented in [44, 45].
We designed a dedicated method for computing steady states and other programs for nonunitary attractors of a given length and a chosen update scheme (synchronous or asynchronous). The originality of our work consists in the exhaustive enumeration of all attractors thanks to the use of ASP, a powerful declarative programming paradigm. The computational framework is based on the AN formalism presuming nondeterministic dynamics. Thanks to the encoding we introduced, and the powerful heuristics developed in modern solvers, we are able to tackle the enumeration of fixed points, cycles and attractors of large models. The major benefit of a such method is to get an exhaustive enumeration of all potential states while still being tractable for models with a hundred of interacting components. As the identification of attractors can give an insight to the longterm behavior of biological systems, tackling this issue is a challenge to which we cared to contribute to. Besides, we hope our work helps open new ways and tools to explore this field.
We plan to extend this work by considering adaptations and optimizations of the approach to address larger models. First, the “projection” feature of Clingo which displays only one answer set when several answer sets contain common predicates, is currently studied in order to filter out repeated attractors, that currently appear multiple times because they are covered by several possible cycles. Another trail consists in returning approximations of the results, that is, sometimes “missing” some answers, but with the benefit of a highly improved performance. Once again, applying various filters to the generated results may avoid redundancy and guide the solving process. Conversely, it may be possible to reduce the incremental aspect of the analysis process, for instance by searching for cycles of size lower than (and not only equal to) a given value, so that the user could directly start with higher values.
Of course, other extensions allowing to tackle other close problems would be of interest. For instance, the attractor inverse problem consists in building or enumerating networks possessing a given set of attractor properties, in order to answer to network inference matters. We would also like to extend these ASPbased methods to study other interesting properties of dynamical patterns such as the enumeration of basins of attraction, gardens of Eden or bifurcations [46].
All programs and benchmarks are available as additional files and at http://www.irccyn.ecnantes.fr/~benabdal/attractors.zip.
All programs and benchmarks are available as additional files and at http://www.irccyn.ecnantes.fr/~benabdal/attractors.zip.
Declarations
Authors’ contributions
EBA designed and implemented the ASP programs. EBA and MF did the formalization. EBA performed the experiments. EBA and MF wrote the paper. MM and OR supervised the work. All authors read and approved the final manuscript.
Acknowledgements
The authors would like to thank the reviewers for their careful reading of the paper and their many insightful comments and suggestions. The authors are grateful to Laurent Trilling for his fruitful comments regarding the use of ASP for dynamic verification.
Competing interests
The authors declare that they have no competing interests.
Availability of data and materials
The datasets analyzed during the current study are available in the GinSim repository, http://ginsim.org/models_repository. All analyzed data during this study are included in this published article, its supplementary information files and at http://www.irccyn.ecnantes.fr/~benabdal/attractors.zip.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
Funding
Our work is funded by the ANR HyClock project in France.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Wuensche A. Genomic regulation modeled as a network with basins of attraction. Pac Symp Biocomput. 1998;3:89–102.Google Scholar
 Huang S, Eichler G, BarYam Y, Ingber DE. Cell fates as highdimensional attractor states of a complex gene regulatory network. Phys Rev Lett. 2005;94(12):128701.View ArticlePubMedGoogle Scholar
 González A, Chaouiya C, Thieffry D. Logical modelling of the role of the hh pathway in the patterning of the drosophila wing disc. Bioinformatics. 2008;24(16):234–40.View ArticleGoogle Scholar
 Albert R, Othmer HG. The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila melanogaster. J Theor Biol. 2003;223(1):1–18.View ArticlePubMedGoogle Scholar
 Stuart A. Kaufmann. The origins of order: selforganization and selection in evolution. Oxford: Oxford University Press; 1993. p. 354.Google Scholar
 Folschette M, Paulevé L, Magnin M, Roux O. Sufficient conditions for reachability in automata networks with priorities. Theor Comput Sci. 2015;608:66–83.View ArticleGoogle Scholar
 Paulevé L. Goaloriented reduction of automata networks. In: International Conference on computational methods in systems biology. Lecture notes in bioinformatics, vol. 9859. Springer; 2016. p. 252–72.Google Scholar
 Thomas R. Regulatory networks seen as asynchronous automata: a logical description. J Theor Biol. 1991;153(1):1–23.View ArticleGoogle Scholar
 Zhang SQ, Hayashida M, Akutsu T, Ching WK, Ng MK. Algorithms for finding small attractors in Boolean networks. EURASIP J Bioinform Syst Biol. 2007;2007(1):1–13.View ArticleGoogle Scholar
 Klemm K, Bornholdt S. Stable and unstable attractors in Boolean networks. Phys Rev E. 2005;72(5):055101.View ArticleGoogle Scholar
 Akutsu T, Kosub S, Melkman AA, Tamura T. Finding a periodic attractor of a Boolean network. IEEE/ACM Trans Comput Biol Bioinform. 2012;9(5):1410–21.View ArticlePubMedGoogle Scholar
 Somogyi R, Greller LD. The dynamics of molecular networks: applications to therapeutic discovery. Drug Discov Today. 2001;6(24):1267–77.View ArticlePubMedGoogle Scholar
 Irons DJ. Improving the efficiency of attractor cycle identification in Boolean networks. Phys D: Nonlinear Phenom. 2006;217(1):7–21.View ArticleGoogle Scholar
 Garg A, Mendoza L, Xenarios I, DeMicheli G. Modeling of multiple valued gene regulatory networks. In: 2007 29th Annual International Conference of the IEEE engineering in medicine and biology society. IEEE; 2007. p. 1398–404.Google Scholar
 Zhao Z, Liu CW, Wang CY, Qian W. Bddbased synthesis of reconfigurable singleelectron transistor arrays. In: Proceedings of the 2014 IEEE/ACM International Conference on computeraided design. IEEE Press; 2014. p. 47–54.Google Scholar
 Dubrova E, Teslenko M. A SATbased algorithm for finding attractors in synchronous Boolean networks. IEEE/ACM Trans Comput Biol Bioinform. 2011;8(5):1393–9.View ArticlePubMedGoogle Scholar
 Mushthofa M, Torres G, Van de Peer Y, Marchal K, De Cock M. ASPG: an ASPbased method for finding attractors in genetic regulatory networks. Bioinformatics. 2041;481.Google Scholar
 Baral C. Knowledge representation, reasoning and declarative problem solving. Cambridge: Cambridge University Press; 2003.View ArticleGoogle Scholar
 Ben Abdallah E, Folschette M, Roux O, Magnin M. Exhaustive analysis of dynamical properties of biological regulatory networks with answer set programming. In: 2015 IEEE International Conference on bioinformatics and biomedicine (BIBM). IEEE; 2015. p. 281–85.Google Scholar
 Paulevé L, Chancellor C, Folschette M, Magnin M, Roux O. Analyzing large network dynamics with process hitting. Log Model Biol Syst. 2014:125–66.Google Scholar
 Skodawessely T, Klemm K. Finding attractors in asynchronous Boolean dynamics. Adv Complex Syst. 2011;14(03):439–49.View ArticleGoogle Scholar
 Berntenis N, Ebeling M. Detection of attractors of large Boolean networks via exhaustive enumeration of appropriate subspaces of the state space. BMC Bioinform. 2013;14(1):1.View ArticleGoogle Scholar
 Calzone L, Fages F, Soliman S. Biocham: an environment for modeling biological systems and formalizing experimental knowledge. Bioinformatics. 2006;22(14):1805–7.View ArticlePubMedGoogle Scholar
 Klarner H, Bockmayr A, Siebert H. Computing maximal and minimal trap spaces of Boolean networks. Nat Comput. 2015;14(4):535–44.View ArticleGoogle Scholar
 de Espanés PM, Osses A, Rapaport I. Fixedpoints in random Boolean networks: the impact of parallelism in the barabásialbert scalefree topology case. Biosystems. 2016;150:167–76.View ArticleGoogle Scholar
 Gelfond M, Lifschitz V. The stable model semantics for logic programming. In: ICLP/SLP; 1988. p. 1070–080.Google Scholar
 Gebser M, Kaminski R, Kaufmann B, Ostrowski M, Schaub T, Wanko P. Theory solving made easy with Clingo 5. Wadern: Schloss DagstuhlLeibnizZentrum fuer Informatik; 2016.Google Scholar
 Dubrova E, Teslenko M. A SATbased algorithm for computing attractors in synchronous Boolean networks; 2009. arXiv preprint arXiv:0901.4448.
 Qu H, Yuan Q, Pang J, Mizera A. Improving bddbased attractor detection for synchronous Boolean networks. In: Proceedings of the 7th AsiaPacific Symposium on Internetware. ACM; 2015.Google Scholar
 Hayashida M, Tamura T, Akutsu T, Zhang SQ, Ching WK. Algorithms and complexity analyses for control of singleton attractors in Boolean networks. EURASIP J Bioinform Syst Biol. 2008;2008(1):1.View ArticleGoogle Scholar
 Thieffry D, Thomas R. Dynamical behaviour of biological regulatory networks—ii. Immunity control in bacteriophage lambda. Bull Math Biol. 1995;57(2):277–97.PubMedGoogle Scholar
 Simao E, Remy E, Thieffry D, Chaouiya C. Qualitative modelling of regulated metabolic pathways: application to the tryptophan biosynthesis in E. coli. Bioinformatics. 2005;21(suppl 2):190–6.View ArticleGoogle Scholar
 Davidich MI, Bornholdt S. Boolean network model predicts cell cycle sequence of fission yeast. PloS ONE. 2008;3(2):1672.View ArticleGoogle Scholar
 Fauré A, Naldi A, Chaouiya C, Thieffry D. Dynamical analysis of a generic Boolean model for the control of the mammalian cell cycle. Bioinformatics. 2006;22(14):124–31.View ArticleGoogle Scholar
 Klamt S, SaezRodriguez J, Lindquist JA, Simeoni L, Gilles ED. A methodology for the structural and functional analysis of signaling and regulatory networks. BMC Bioinform. 2006;7(1):1.View ArticleGoogle Scholar
 Mbodj A, Junion G, Brun C, Furlong EE, Thieffry D. Logical modelling of drosophila signalling pathways. Molecular BioSyst. 2013;9(9):2248–58.View ArticleGoogle Scholar
 AbouJaoudé W, Monteiro PT, Naldi A, Grandclaudon M, Soumelis V, Chaouiya C, Thieffry D. Model checking to assess thelper cell plasticity. Front Bioeng Biotechnol. 2014;2.Google Scholar
 Chaouiya C, Naldi A, Thieffry D. Logical modelling of gene regulatory networks with GINsim. Bact Mol Netw: Methods Protoc. 2012:463–79.Google Scholar
 Paulevé L. Pint, a static analyzer for dynamics of automata networks. In: 14th International Conference on computational methods in systems biology (CMSB 2016); 2016.Google Scholar
 Naldi A, Monteiro PT, Müssel C, Kestler HA, Thieffry D, Xenarios I, SaezRodriguez J, Helikar T, Chaouiya C, et al. Cooperative development of logical modelling standards and tools with colomoto. Bioinformatics. 2015;013.Google Scholar
 Chatain T, Haar S, Jezequel L, Paulevé L, Schwoon S. Characterization of reachable attractors using petri net unfoldings. In: International Conference on computational methods in systems biology. Springer. p. 129–42.Google Scholar
 Thomas R. Boolean formalization of genetic control circuits. J Theor Biol. 1973;42(3):563–85.View ArticlePubMedGoogle Scholar
 Kauffman SA. Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol. 1969;22(3):437–67.View ArticlePubMedGoogle Scholar
 Gershenson C. Updating schemes in random Boolean networks: do they really matter. In: Artificial life IX Proceedings of the Ninth International Conference on the simulation and synthesis of living systems. MIT Press; 2004. p. 238–43.Google Scholar
 Noual M, Sené S. Synchronism versus asynchronism in monotonic Boolean automata networks. Nat Comput. 2017. doi:10.1007/s1104701696088.Google Scholar
 FippoFittime L, Roux O, Guziolowski C, Paulevé L. Identification of bifurcations in biological regulatory networks using answerset programming. In: Constraintbased methods for bioinformatics Workshop; 2016.Google Scholar