BDD-BASED ALGORITHM FOR SCC DECOMPOSITION OF EDGE-COLOURED GRAPHS

. Edge-coloured directed graphs provide an essential structure for modelling and analysis of complex systems arising in many scientiﬁc disciplines (e.g. feature-oriented systems, gene regulatory networks, etc.). One of the fundamental problems for edge-coloured graphs is the detection of strongly connected components, or SCCs. The size of edge-coloured graphs appearing in practice can be enormous both in the number of vertices and colours. The large number of vertices prevents us from analysing such graphs using explicit SCC detection algorithms, such as Tarjan’s, which motivates the use of a symbolic approach. However, the large number of colours also renders existing symbolic SCC detection algorithms impractical. This paper proposes a novel algorithm that symbolically computes all the monochromatic strongly connected components of an edge-coloured graph. In the worst case, the algorithm performs O ( p · n · log n ) symbolic steps, where p is the number of colours and n is the number of vertices. We evaluate the algorithm using an experimental implementation based on binary decision diagrams (BDDs). Speciﬁcally, we use our implementation to explore the SCCs of a large collection of coloured graphs (up to 2 48 ) obtained from Boolean networks – a modelling framework commonly appearing in systems biology.


Introduction
In many scientific disciplines, the processing of massive data sets represents one of the most important computational tasks. A variety of these data sets can be modelled in terms of very large multi-graphs, augmented by a specific collection of application-dependent edge attributes. These attributes are often abstractly referred to as colours, and the resulting formalism is called an edge-coloured graph [BJG97,BCLF79]. Geographic information systems, telecommunications traffic, or internet networks are prime examples of data that are best represented as such edge-coloured graphs.
For instance, in social networks, coloured edges can be used to link together groups of nodes related by some specific criteria (Sports, Health, Technology, Religion, etc.). In software engineering, one often speaks about feature-oriented systems [CHS + 10]. In this case, colours represent possible combinations of features, altering the system's behaviour.
Our interest in processing huge edge-coloured graphs is primarily motivated by applications taken from systems biology [BBK + 12, GGG + 17] and genetics [Dor94] where we have to deal not only with giant graphs as measured by the number of vertices and edges but also with large sets of colours. In this case, the graph colours represent valuations of numerous parameters that influence the dynamics of a biological system [BBK + 12, BPC + 10, RCB05].
Fundamental graph algorithms such as breadth-first search, spanning tree construction, shortest paths, decomposition into strongly connected components (SCCs), etc., are building blocks of many practical applications. For the edge-coloured graphs, the primary research focus so far has been on some of the "classical" coloured graph problems, like the determination of the chromatic index, finding sub-graphs with a specified colour property (the coloured version of the k-linked problem), alternating edge-coloured cycles and paths, rainbow cliques, monochromatic cliques and cycles, etc. [ADF + 08, AA07, AG97, BJG97, TW07,KL08].
To the best of our knowledge, we are not aware of any work on SCC decomposition specifically for edge-coloured graphs, even though this problem has many important applications. For example, in biological systems, strongly connected components represent the so called attractors of the system. In this case, a specific focus is given to terminal (or bottom) SCCs, but non-terminal (transient) SCCs can also be detrimental to the system's long-term behaviour [PIS + 21]. Overall, SCCs play an essential role in determining the system's biological properties, since they may correspond, for example, to the specific phenotypes expressed by a cell [CC16].
The valuation of parameters (e.g. the presence of certain genes or external stimuli) in such systems is then represented as edge colours in the state-transition graph. The knowledge of SCCs and how their structure depends on parameters is vital for understanding various biological phenomena [DAERR16,LWA + 16]. Other applications where investigation of attractors is crucial include predictions of the global climate change [SRR + 18] or predictions of spreading of infectious diseases such as COVID-19 [Mat20].
There is a serious computational problem related to the processing of massive edgecoloured graphs (or even the non-coloured ones) that significantly affects the tractability of SCC decomposition. The graphs often cannot be handled using standard (explicit) representations, since they are too large to be kept in the main memory. Various approaches have been considered to deal with such giant graphs: distributed-memory computation, symbolic data structures for graph representation, or storing the graphs in external memory. We review these approaches in more detail in the related work section.
In [BBB + 17, BBP + 19] we have initially attacked the SCC decomposition problem for massive edge-coloured graphs by developing a parallel, semi-symbolic algorithm for detection of bottom SCCs. The algorithm uses symbolic structures to represent sets of parameters, while the graph itself is represented explicitly. However, the results have shown that the parallel semi-symbolic algorithm is often not sufficient to tackle graphs representing real-world problems practically. These findings have motivated us to propose a new, entirely symbolic approach.
In this paper, we consider edge-coloured multi-digraphs, i.e., multi-digraphs such that each directed edge has a colour and no two parallel (i.e., joining the same pair of vertices) edges have the same colour. Here, we refer to such graphs simply as coloured graphs. For coloured graphs, we can define several notions of strongly connected components involving colours. We consider the simplest case, where the SCCs are monochromatic, that is all their edges have the same colour [Kir14]. This choice is motivated by the application in systems biology, as mentioned above.
Contribution. We propose a novel fully symbolic algorithm for detecting all monochromatic strongly connected components in a coloured graph. This algorithm is in practice significantly faster than what is achievable by naïvely executing a symbolic SCC decomposition algorithm for each colour separately. This is because in many applications, the edges are largely shared among individual colours [BBK + 12] and our algorithm is capable of exploiting this fact. The algorithm conceptually follows the lock-step reachability approach by Bloem et al. [BGS00] for purely monochromatic digraphs. The key new ingredients behind our algorithm are a careful orchestration of the forward and backward reachability for different colours, and a colour-aware selection of the pivot set.
Structure of the paper. In Section 1, we recall the notions of strongly connected components and edge-coloured digraphs, and we state the coloured SCC decomposition problem. In Section 2, we first briefly introduce the forward-backward decomposition algorithm and the lock-step algorithm for monochromatic graphs. After that, we present the coloured SCC decomposition algorithm together with the proof of correctness and complexity analysis.
In Section 3, we introduce Boolean networks, discuss their symbolic encoding, and show how they can be translated into coloured graphs suitable for SCC-decomposition. Subsequently, Section 4 discusses several practical improvements to the main algorithm (saturation, trimming, and parallelism) which help it scale to larger models, and thus be more practically viable. Finally, Section 4 evaluates the main algorithm (including the improved variants) using a collection of large, real-world Boolean networks. A conclusion is provided in the last section.
This article is an extended version of an article that appeared in the Proceedings of TACAS 2021 [BBPŠ21]. We extend the information provided in the TACAS Proceedings with more in-depth technical details of the algorithm and related proof, including explanation of its key steps. Moreover, we extend the implementation and evaluation sections to give the reader more information on how the performance of the algorithm can be improved, and how the algorithm performs using a variety of real-world case studies.
Related Work. The detection of SCCs in (monochromatic) digraphs is a well-known problem computable in linear time. Best serial (explicit) algorithms are Kosaraju-Sharir [Sha81] and Tarjan [Tar72], which are both inherently based on depth-first search. However, these algorithms do not scale for large graphs, e.g., those encountered in model-checking, when using explicit graph representation. Therefore, alternative approaches to such SCC decomposition have been proposed (e.g. I/O efficient, parallel, or symbolic algorithms).
The algorithm of Jiang [Jia93] gives an I/O-efficient alternative based on a combination of depth-first and breadth-first search.
Efficient parallel, distributed-memory algorithms avoid the inherently sequential DFS step [Rei85] in several different ways. The Forward-Backward algorithm [FHP00] employs a divide-and-conquer approach relying on picking a pivot state and splitting the graph in three independent (SCC-closed) parts. The approach of Orzan [Orz05] uses a different distribution scheme called a colouring transformation, employing a set of prioritised colours to split the graph into many parts at once. The OWCTY-Backward-Forward (OBF) approach is proposed in [BCVDP11]. It recursively splits the graph in a number of independent sub-graphs called OBF slices and applies to each slice the One-Way-Catch-Them-Young (OWCTY) technique. In [SRM14], the authors utilise variants of the Forward-Backward and Orzan's algorithms for optimal execution on shared-memory multi-core platforms. Finally, Bloemen et al. [BLvdP16] present an on-the-fly parallel algorithm utilising a swarm of DFS searches, showing promising speed-up for large graphs containing large SCCs. On another end, GPU-accelerated approaches to computing SCCs have been addressed for example in [BBBv11,HRO13,LZCY14,WKB14]. Computing SCCs of (monochromatic) digraphs symbolically is another way to handle giant graphs and has been thoroughly explored in literature. As in the case of efficient parallelisation, depth-first search is not feasible in the symbolic framework [GPP08]. In consequence, many DFS-based algorithms cannot be easily revised to work with symbolically represented graphs. An algorithm based on forward and backward reachability performing O(n 2 ) symbolic steps was presented by Xie and Beerel in [XB00]. Bloem et al. present an improved O(n · log n) algorithm in [BGS00]. Finally, an O(n) algorithm was presented by Gentilini et al. in [GPP03,GPP08]. This bound has been proven to be tight in [CDHL18]. In [CDHL18], the authors argue that the algorithm from [GPP03] is optimal even when considering more fine-grained complexity criteria, like the diameter of the graph and the diameters of the individual components. Ciardo et al. [ZC11] use the idea of saturation [CMS06] to speed up state exploration within the Xie-Beerel algorithm, and show a saturation-based technique for computing the transitive closure of the graph's edge relation.
Besides these generic algorithms, there have also been symbolic SCC decomposition methods to deal with large graphs generated specifically by Boolean networks [MPQY19,YMPQ19]. However, these primarily target detection of bottom SCCs. Methods in this area are also often incomplete, for example focusing on detection of single-state or small bottom SCCs [ZHA + 07]. As such, they generally perform better than an exhaustive symbolic SCC detection in their respective application domains, but are inherently limited in scope.

Problem Definition
As we have already stated in the introductory section, the SCC decomposition problem for edge-coloured graphs has remained mostly unexplored until now. We thus start this paper by introducing and formalising the notion of coloured SCC decomposition itself and state some of its basic properties.
Before giving exact definitions, it might be instructive to discuss the substance of the coloured SCC decomposition intuitively. There are several ways of capturing the notion of a "coloured connected component". One of them is that of a colour-connectivity first introduced by Saad [Saa92]. It is based on alternating paths in which successive edges differ in colour. However, there is no unique, universally acceptable notion of a coloured component.
In the biological applications we have in mind (i.e. Boolean networks), we want to identify a coloured component as a coloured collection of SCCs-a collection where for every colour there is a set of all relevant monochromatic SCCs. Such a setting leads us to represent SCCs in the form of a relation. To that end, we first introduce such a relation for monochromatic graphs (Section 1.1) and afterwards extend it to edge-coloured graphs (Section 1.2). The relation-based approach gives us also the advantage of allowing a feasible symbolic encoding of the problem.
1.1. Graphs and Strongly Connected Components. Let us first recall the standard definitions of a directed graph and its strongly connected components: We are going to use the word graph to mean directed graph in the following. We write u → v when (u, v) ∈ E and u → * v when (u, v) ∈ E * , the reflexive and transitive closure of E. We say that v is reachable from u if u → * v. The reachability relation allows us to decompose a graph into strongly connected components, defined as follows: If the graph G is clear from the context, we can simply write SCC(v). A set of vertices S ⊆ V is said to be SCC-closed if every SCC W is either fully contained inside S (W ⊆ S), or in its complement (W ⊆ V \ S). Notice that given a vertex v, the set of all vertices reachable from v, as well as the set of all vertices that can reach v, are both SCC-closed.
A pivotal problem in computer science is to find the SCC decomposition of G. As mentioned above, we represent the decomposition in the form of an equivalence relation R scc such that the individual SCCs are exactly the equivalence classes of R scc . The relation-based formulation of the SCC decomposition problem is the following: Note that SCC(u) can be obtained by fixing the first attribute of R scc , i.e. SCC(u) = {v | (u, v) ∈ R scc }. We refer to such operation as section and denote it in the following way: SCC(u) = R scc (u, _) (the concept is properly formalised later as part of Fig. 1). Here, u is the specific value of an attribute at which the section is taken, and _ is used in place of the attributes that remain unchanged. Such notation naturally extends to arbitrary relations.
1.2. Coloured SCC Decomposition Problem. We now lift the formal framework to the coloured setting. An edge-coloured graph can be seen as a succinct representation of several different graphs, all sharing the same set of vertices. To emphasise the difference from the standard graphs (i.e. Definition 1.1), we sometimes call the standard graphs monochromatic.
Definition 1.4. An edge-coloured directed multi-graph (coloured graph for short) is a tuple G = (V, C, E) where V is a set of vertices, C is a set of colours and E ⊆ V × C × V is a coloured edge relation.
We also write u c − → v whenever (u, c, v) ∈ E and use c − → * to denote the reflexive and transitive closure of there is a path from u to v using only c-coloured edges. By fixing a colour c ∈ C and keeping only the c-coloured edges (with the colour attribute removed), we obtain a monochromatic graph G(c) = (V, {(u, v) | (u, c, v) ∈ E}). We call this graph the monochromatisation of G with respect to c. Intuitively, one can view the elements of C as a type of graph parametrisation where the edge structure of the graph changes based on the specific c ∈ C.
The SCC decomposition relation R scc is extended to the coloured SCC decomposition relation R scc by relating every colour c ∈ C with all SCCs of the monochromatisation G(c). In consequence, the SCC decomposition problem is then lifted to the coloured SCC decomposition problem as follows: Problem 1.5 (Coloured SCC decomposition). Given a coloured graph G = (V, C, E), find the coloured SCC decomposition relation From this definition, we can immediately observe the following properties about the relationship of R scc with the terms which we have defined before: and R scc are symmetric with regards to V ). From this, it should be immediately apparent that R scc contains all components of the underlying monochromatisations.

Algorithm
Conceptually, our algorithm follows the lock-step reachability approach by Bloem [BGS00] for monochromatic graphs. The lock-step algorithm itself is based on the basic forwardbackward decomposition algorithm [XB00]. In this section, we first briefly introduce these two algorithms to explain better the key ideas behind our approach and, in particular, to explain the main difficulties encountered in employing the concepts of these algorithms to edge-coloured graphs. Although the algorithms were originally presented as producing a set of SCCs, we reformulate them slightly using the equivalent relation-based approach as explained in the previous section. After that, we present the coloured SCC decomposition algorithm. However, before we dive into the algorithmics, let us first briefly discuss the computation model we are using.
2.1. Symbolic Computation Model. As a complexity measure of our algorithm, we consider the number of symbolic steps, or more specifically, symbolic set and relation operations that the algorithm performs. As is customary, we assume that sets of vertices (V ) and colours (C) can be represented symbolically (for example, using reduced ordered binary decision diagrams [Bry86]) as well as any relations over these sets. In particular, we often talk about coloured vertex sets, by which we mean the subsets of V × C.
Aside from normal set operations (union, intersection, difference, product and element selection), we also require some basic relational operations, all of which we outline in Figure 1. These extra operations tend to appear in other applications as well (such as symbolic model checking [BCM + 92]), and are thus typically already available in mature symbolic computation packages.
Finally, there are several derived operators that are partially specific to our application to coloured graphs. However, these can be constructed using standard set and relation operations. The intuitive meaning of the derived operators is as follows: Colours returns all the colours that appear in the given coloured vertex set. Pre and Post compute the preand post-image of a (monochromatic or coloured) set of vertices, i.e. the set of successors or predecessors of all the vertices in the given set, respectively. Finally, Join takes a coloured vertex set A and computes the set Figure 1: Summary of symbolic operations that appear in the presented algorithms. The derived operations can be implemented using the standard and relational operations. However, typically they also have a slightly more efficient direct implementations.

Forward-Backward Algorithm.
To symbolically compute the SCCs of a graph G = (V, E), Xie and Beerel [XB00] observed that for The algorithm thus picks an arbitrary pivot v ∈ V , and divides the vertices of the graph into four disjoint sets: . This is illustrated graphically in Figure 2 (left). The set W is then immediately reported as an SCC of the graph, and added into the component relation: It is easy to see that every other SCC is fully contained within one of the three remaining sets (they are SCC-closed), and the algorithm thus recursively repeats this process independently in each set.
The correctness of the algorithm follows from the initial observation and the fact that every vertex eventually appears in W (either as a pivot or as a result of F ∩ B). In the worst case, the algorithm performs O(|V | 2 ) symbolic steps, since every vertex is picked as  Figure 2: Illustration of the difference between the forward-backward algorithm (left) and the lock-step algorithm (right). On the left, we fully compute both backward (B) and forward (F ) reachable sets from the pivot v, identifying W as F ∩ B. On the right, without loss of generality, assume F is fully computed first. It thus becomes converged (Con) and the computation of B (N on) is stopped before it is fully explored.
a pivot at most once and the computation of F and B requires at most O(|V |) Pre/Post operations.

Lock-step Algorithm.
To improve the efficiency of the forward-backward algorithm, the lock-step approach [BGS00] uses another important observation: To compute W , it is not necessary to fully compute both F and B; only the smaller (in terms of diameter) of the two sets needs to be entirely known. With this observation, the computation of F and B can be modified in the following way: Instead of computing F and B one after the other, the computation is interleaved in a step-by-step manner (dovetailing). When one of the sets is fully computed, the computation of the second set is stopped. Let us call the computed set converged and denote it by Con, and the unfinished set non-converged and denote it by Non. This situation is illustrated in Figure 2 (right). However, even when Con is fully known, we still need to finish the computation of states in Non that are inside Con to discover the whole component W . This is necessary if there are vertices w in W whose forward distance from v (i.e. the length of the path v → * w) is short while their backward distance (the length of the path w → * v) is long, or vice versa. Such vertices are thus only discovered in one of the two reachability procedures and still need to be discovered by the other one to identify the whole component. However, an important observation is that only the vertices already inside Con need to be considered in this phase.
After this, the SCC can be identified and reported just as in the forward-backward algorithm. Finally, the recursion now continues in sets Con \ W and V \ Con. This is due to Non being not fully computed; we cannot guarantee that no SCC overlaps outside of Non (Non is not necessarily SCC-closed).
The algorithm is still correct because every vertex is eventually either picked as a pivot or discovered in some W . Furthermore, due to the way Con and Non are computed guarantees that W is still a whole SCC. In terms of complexity, the algorithm performs O(|V | · log |V |) symbolic steps in the worst case. To see why this is true, we may observe that every vertex appears in W exactly once, and that the smaller of the two sets Con \ W and V \ Con, let us call it S, is always smaller than |V | 2 . The authors then argue that the price of every iteration can be attributed (up to a multiplicative constant) to the vertices in S ∪ W and that every vertex appears in S at most O(log |V |)-times.
2.4. Coloured Lock-step Algorithm. When developing an algorithm for coloured graphs, one needs to deal with multiple challenges which do not appear for monochromatic graphs and require careful consideration. In the following, we refer to the pseudocode in Algorithm 1. An important observation is that the structure of components in the graph can change arbitrarily with respect to the graph colours. In consequence, our algorithm cannot simply operate with sets of graph vertices as the normal algorithm would. To that end, we use the notion of coloured vertex sets as introduced in Section 2.1 where the symbolic operations we perform on these sets have been described.
Pivot selection. Initially, the algorithm starts with all vertices and colours, i.e. the full set V × C. However, as the components are discovered, the intermediate results V may contain different vertices appearing only for certain subsets of C. As a result, we often cannot pick a single pivot vertex that would be valid for all considered colours. Instead, we aim to pick a pivot set P ⊆ V × C such that for every colour that still appears in V, the set contains exactly one vertex. Alternatively, one can also view the pivot set as a (partial) function from C to V . This is done in the Pivots function. In the following discussion of the algorithm, we write c-coloured pivot to mean the vertex u such that (u, c) is found in the coloured set returned by Pivots in the current iteration (for all colours still present in V).
Please note that the presented Pivots routine is rather naive, as it has to explicitly iterate all the pivot vertices, whose number can be substantial in the worst case. However, as presented, it should be easy to implement for basically any type of coloured graphs, regardless of the underlying representation. In the implementation section, we show Pivots can be re-implemented in the domain of BDDs such that it is guaranteed to always require only O(log |V |) symbolic operations.
Coloured lock-step (phase one). The lock-step reachability procedure also cannot operate as in a standard graph. First of all, there can be colours where the forward reachability converges first, as well as colours where this happens for backward reachability. The algorithm thus has to account for both options simultaneously. Second, for each colour, the reachability can converge in a different number of steps. To deal with this problem, we introduce the F lock and B lock variables. These store the mutually disjoint sets of colours for which the forward and backward reachability procedures have already converged. The lock-step procedure then terminates when F lock and B lock contain all the colours that appear in V.
Throughout the algorithm, we keep track of several coloured-set variables. The first two, F and B, represent the forward and backward reachable sets, respectively. This means that for every colour Furthermore, if c ∈ F lock then F contains exactly all such pairs; similarly for B lock and B.
We say that a coloured vertex pair (v, c) has been forward expanded or backward expanded in the current iteration of the algorithm, if there has been a call to the Post or Pre symbolic operation with (v, c) being an element of the coloured set argument. To track which reachable coloured vertices are to be expanded later, also called the frontiers of the reachability sets,  Algorithm 1: Symbolic Coloured SCC Decomposition The frontier of F is the union F open ∪ F paused . The sets F open and F paused are disjoint: F open involves those colours for which the lock-step reachability procedure has not finished yet, i.e. the colours that are neither in F lock nor in B lock , while F paused represents the part of the frontier whose exploration is currently paused due to the fact that its colours are in B lock . Note that there may be no pair (v, c) of the forward frontier with c ∈ F lock as that means that the exploration of the c-coloured forward-reachable set is complete. A symmetric role is played by the sets B open and B paused .
In the first while loop (lines 10-21), we compute the reachability sets in the lock-step manner. Once a reachability set is completed for some colours (i.e., there are no vertices to expand with those colours), we add the colours to the corresponding F lock or B lock variable. Note that we ensure that F lock and B lock remain disjoint even if the two reachability procedures converged at the same time for certain colours-see line 14. We use F lock and B lock to split the newly computed frontier sets into the parts that are to be expanded in the next iteration (F open , B open ) and the parts currently left unexpanded (F paused , B paused ).
Note that during the computation of Post and Pre on lines 11 and 12, we intersect the resulting set with V. This step is not necessary for correctness, but as the algorithm divides V × C into smaller sets in each recursive call to Decomposition, it can happen that the set of states reachable from V is substantially larger than V itself. In such cases, this intersection effectively restricts the computation of Post and Pre to the sub-graph of G induced by V.
Component identification (phase two). After the first while loop terminates, we compute the set Con that is an analogue for the converged set of the original lock-step algorithm (line 22). As already suggested above and unlike the original algorithm, this set cannot be just F or B, but is instead a mixture of both, depending on the converged colours. To compute this set, we use the F lock and B lock variables.
Once Con is computed, F open and B open are restarted using the converged portion of F paused and B paused (lines 23 and 24). The second while loop (lines 25-30) can then complete the unfinished forward and backward reachability set, now restricted to the inside of the converged set. The intersection of F and B then forms a coloured set W with the property that for all c ∈ Colours(V), W(_, c) is a strongly connected component of G(c). We create the corresponding relation using the Join operation, add this relation to the resulting R scc , and recursively call the whole procedure with V \ Con and Con \ W as the base sets.
Comments on the coloured approach. Let us note that there is possibly another approach to processing coloured graphs. Instead of trying to work with all colours still appearing in the coloured vertex set at once, we could fork a new recursive procedure whenever the colour set splits due to the differences in the graph structure. For example, instead of picking multiple coloured vertices as pivots, one could pick a single vertex with a valid subset of colours and then address the remaining colours in a separate recursive call. Similarly, instead of a single recursive Decomposition call with Con \ W, we could consider two calls, one with the F portion of Con and the other with the B portion of Con (note that these are colour-disjoint since each colour can converge only in one of the two sets).
While such an approach could be to some extent beneficial in a massively parallel environment where each recursive call can be executed independently on a new CPU, the amount of forking in large systems will soon become unreasonable. More importantly, it defeats the purpose of the symbolic representation, which aims to minimise the number of symbolic operations.  Example. The execution of one iteration of the algorithm is illustrated in Figure 3. Here, we have an edge-coloured graph with six vertices and two colours (red and blue). The top-most picture represents the initial situation after we have chosen the pivots; in this case, {(b, blue), (b, red )}. The next four rows illustrate the first phase (the first while loop) of the algorithm. After the second iteration of the loop, the blue colour becomes locked in F lock , and thus (f, blue) is not expanded in the backward reachability procedure. This is illustrated by its dashed outline. After the third iteration of the loop, the red colour becomes locked in B lock and thus the first phase ends. In the second phase, both reachability procedures continue from the paused coloured vertices (dashed outlines); the result is seen in the fifth row. The intersection of the two reachable sets (i.e. the coloured set W) is then illustrated in the bottom-most picture. The algorithm would now continue with the coloured sets V \ Con = {(a, blue), (d, blue)} and Con \ W = {(c, red )}.

Correctness and Complexity of the Coloured Lock-step Algorithm.
Theorem 2.1. Let G = (V, C, E) be a coloured graph. The coloured lock-step algorithm terminates and computes the coloured SCC decomposition relation R scc .
Proof. We first show that the set W computed in line 31 indeed contains one SCC for every colour c ∈ Colours(V) and that the recursive calls of Decomposition preserve the property that V is SCC-closed with respect to all colours. Let us assume that V is SCC-closed and let us take an arbitrary c ∈ Colours(V). The function Pivots chooses a set that contains exactly one pair whose colour is c, let us call this pair (v, c). Let us further assume that c is assigned into F lock first (the case with B lock is completely symmetric).
Let us now choose an arbitrary vertex w such that v and w are in the same SCC of G(c), i.e. v → * w and w → * v. As the first while loop finishes, F contains all the pairs of the form (u, c) ∈ V where u is reachable from v in G(c). Thus, it also contains (w, c) due to the fact that V is SCC-closed. Now, either (w, c) ∈ B, or there exists a vertex x such that w → * x, x → * v in G(c) and x ∈ B paused . This means that (w, c) is added to B in the second while loop. In both cases, both (v, c) and (w, c) are then added to W. As the vertex choices were arbitrary, this proves that the SCC of v in G(c) is contained in W. Furthermore, if (y, c) ∈ W for an arbitrary y, then v → * y and y → * v in G(c), which means that y is in SCC(G(c), v). This proves that W contains exactly one SCC for every colour c ∈ Colours(V).
We now argue that Con is SCC-closed with respect to all colours. This immediately implies that both V \ Con and Con \ W are SCC-closed. Let us assume that there is a colour c ∈ Colours(V) and two vertices v, w in the same SCC of G(c) such that (v, c) ∈ Con, but (w, c) ∈ Con. Let us assume that c ∈ F lock (as above, the case of B lock is completely symmetrical). Then (v, c) ∈ F after the first while loop finishes. This also means that (w, c) ∈ F as the forward reachability procedure is completed for c and thus (w, c) ∈ Con, a contradiction.
What remains is to show that the algorithm terminates and that every SCC is eventually found. Termination is trivially proved by the fact that size of the set V always decreases in recursive calls: both W and Con are non-empty because they contain the initial pivot set as a subset. Clearly, a representant of every SCC of every monochromatisation G(c) is eventually chosen as a pivot. Together with the above reasoning, this implies that the algorithm is correct. Theorem 2.2. Let |V | be the number of vertices in the coloured graph and let |C| be the number of colours. The coloured lock-step algorithm performs at most O(|C| · |V | · log |V |) symbolic steps.
Proof. Let us first note that all the derived operations defined in Figure 1 use only a constant number of the basic symbolic operations. As we are considering asymptotic complexity here, we can view all the operations in Figure 1 as elementary symbolic steps. We first make the observation that each vertex may be chosen as a part of the pivot set at most |C| times. Clearly, once a vertex is included in the pivot set with a set of colours C , then, {v} × C is a subset of first Con, and later W (due to the monotonicity of the construction of F and B). Therefore, the elements of {v} × C do not appear in subsequent recursive calls. Since a single vertex-colour pair cannot be returned by Pivots twice, it means that the total cumulative complexity of all the calls to the Pivots routine is bounded by O(|C| · |V |). We can therefore exclude them from the rest of the complexity analysis.
We now consider the complexity of a single call to Decomposition without the subsequent recursive calls. Let us now select one of the colours for which the lock-step reachability procedure (lines 10-21) finished last, i.e., one of the colours that have been added to F lock or B lock in the final iteration of the loop. Let us call this colour c. Recall that σ 2 (c, X ) is the set of vertices with colour c in a coloured set X .
Let us denote by W the monochromatic SCC discovered for c, i.e. W := σ 2 (c, W), and let S be the smaller of σ 2 (c, V \ Con) and σ 2 (c, Con \ W). Clearly S contains at most |V |/2 vertices. Let k = |S ∪ W |. We now argue that the number of symbolic steps in a given call (without the recursive calls) is bounded by O(k). This is because in a lock-step algorithm, the call to Decomposition must explore the discovered SCC itself (i.e. W ), and the smaller of the forward or backward reachable sets from this SCC (i.e. S) -intuitively, its complexity should be thus bounded by the size of these two sets.
Assume w.l.o.g. that c ∈ F lock (a completely symmetric argument solves the case c ∈ B lock ). Then after the first while loop finishes, we have σ 2 (c, Con) = σ 2 (c, F). If S is σ 2 (c, Con \ W) then k is the size of σ 2 (c, F) (and thus also σ 2 (c, Con)), since σ 2 (c, F) consists of σ 2 (c, Con \ W) (the set S) and σ 2 (c, W) (the discovered SCC). Each iteration of the first while loop puts at least one vertex with colour c into F; otherwise c would not have finished in the last iteration. This means that the loop runs for at most k iterations. This also means that the size of σ 2 (x, X ) for all colours x and X ∈ {F, B} is also bounded by k after the first while loop finishes, which in turn means that the second while loop cannot make more than O(k) steps.
If S is σ 2 (c, V \ Con) instead, let us define B := σ 2 (c, B) right after the first while loop has finished. We know that B ⊆ S ∪ W : if a vertex v was in B \ S, then it would have to be in Con (i.e. (v, c) ∈ Con). Due to our initial assumption of c ∈ F lock (w.l.o.g), we then also have (v, c) ∈ F which dictates v ∈ W . Consequently, we see that any vertex v ∈ B must be either in S or in W , arriving at B ⊆ S ∪ W . Again, each iteration of the first while loop puts at least one vertex with colour c into B; otherwise c would have been in B lock before it appeared in F lock . Similarly to the previous case, this means that both while loops run for at most O(k) steps.
The rest of the argument uses amortised reasoning, in a way similar to the proof in [BGS00]. Note that each vertex is going to be an element of the set W as described above at most |C| times (once for each colour). Furthermore, each vertex is going to be an element of the set S as described above at most |C| · log |V | times: for each colour, the vertex can be an element of the smaller of the two sets at most log |V | times. As the cost of each single call can be charged to the vertices in S ∪ W as explained above, it is sufficient to charge each vertex the total cost of |C| + |C| · log |V |. Together, this means that the total number of symbolic steps is bounded by O(|C| · |V | · log |V |).
Note that the upper bound established by Theorem 2.2 is no better than the one we would get if we split the coloured graph into its monochromatic constituents and processed each separately using the original lock-step algorithm [BGS00]. We remark, however, that the practical complexity of the coloured approach can be much smaller. Indeed, the complexity analysis in the previous proof focused on a single colour, omitting the fact that SCCs for many other colours are found at the same time. In cases where the edges are largely shared among the colours, which is true in many applications, the coloured algorithm has the potential to significantly outperform the parameter-scan approach. The situation is similar to that of the coloured model checking; see the observations made in [BBK + 12].

Symbolic Computation with Boolean Networks
The algorithm as presented in the previous section is completely agnostic to the properties of the underlying system, as long one provides an implementation of all the necessary symbolic operations. However, to empirically test its performance, we need to pick such an implementation, which typically entails analysis of a specific class of systems.
In this paper, we consider Boolean networks [Kau69, RCB05, SKI + 20, Tho73], specifically asynchronous Boolean networks, which represent a popular discrete modelling framework in systems biology [BČŠ13, GBS + 15]. Due to incomplete biological knowledge, the dynamics of a Boolean network can by often only partially known. This uncertainty can be then captured using coloured directed graphs. In this section, we introduce Boolean networks and show how they can be translated into coloured graphs suitable for SCC-decomposition.
Asynchronous Boolean networks are especially challenging for symbolic analysis. It is a well-known fact, that using symbolic structures (e.g BDDs) to explore very large state spaces gives good results for synchronous systems, but shows its limits when trying to tackle asynchronicity (see e.g. [CT05]).
3.1. Boolean networks with inputs. A Boolean network (BN), as the name suggests, consists of n Boolean variables s 1 , . . . , s n which together describe the state of the network. The dynamics of the network can also depend on additional m Boolean inputs c 1 , . . . , c m (sometimes also called constants, or logical parameters), whose value is assumed to be fixed, but generally unknown. The valuations of these inputs correspond to the colours of our Kripke structure.
Each network variable s i is equipped with a Boolean update function b i : {0, 1} n × {0, 1} m → {0, 1} that updates the variable based on the state of the network, and the values of its inputs. We assume that the variables are updated asynchronously, meaning that during every state transition, exactly one variable is updated.  Note that in practice, we often work within a subset of biologically relevant colours, denoted as V alid (i.e. not every possible valuation of c 1 , . . . , c m may be biologically admissible). In the algorithms, this is implicitly reflected such that the set of all possible colours C corresponds to the set V alid instead of the set of all possible valuation (i.e. {0, 1} m ) if demanded by the application at hand.

3.2.
Partially specified Boolean networks. A Boolean network with inputs allows us to easily encode a wide range of biochemical systems in a machine friendly format. However, for systems with a high degree of uncertainty, it often fails to capture this uncertainty in a way understandable to a human reader.
To mitigate this issue, we consider partially specified Boolean networks that allow us to explicitly mark parts of the update functions as unknown. Specifically, let us assume that f , . . . are symbols standing in for some uninterpreted (fixed but arbitrary) Boolean functions (here, a i denotes their arity). A partially specified Boolean network then consists of n Boolean variables and p uninterpreted Boolean functions. In such a network, every update function b i is specified as a Boolean expression that can use the function symbols f 1 , . . . , f p .
This type of formalism is often easier to comprehend, as the uncertainty in dynamics is tied to the update functions instead of inputs (if desired, input can be still expressed using uninterpreted functions of arity zero). It is not immediately clear how such a network should be represented symbolically though.
One option is to translate a partially specified network into a BN with inputs. Any uninterpreted function f in the j-th row of its truth table. Formally, this translation can be achieved using a repeated application of the following expansion rule: f (α 1 , . . . , α a ) ≡ (α 1 ⇒ f 1 (α 2 , . . . , α a )) ∧ (¬α 1 ⇒ f 2 (α 2 , . . . , α a )) Here, f 1 and f 2 are fresh uninterpreted functions of arity a − 1, and α i are arbitrary Boolean expressions. Using this rule, we can always convert a partially specified network to a Boolean network with inputs. The number of inputs will be exponential with respect to the arity of the employed uninterpreted functions though (since each application of the rule replaces one uninterpreted function with two, and the depth of the recursive expansion is the arity a).
For example, consider the following partially specified Boolean network: It uses three uninterpreted Boolean functions f 2 , and f (0) 3 . After performing the aforementioned expansion, and simplifying the resulting expressions slightly for readability, we obtain the following network with logical inputs: Vol. 18:1 BDD-BASED ALGORITHM FOR SCC DECOMPOSITIONOF EDGE-COLOURED GRAPHS 38:17 Here, each c i j corresponds to one truth table row of f i , such that j describes the input vector corresponding to said row (i.e. c 1 [0,1] represents the value of f 1 (0, 1)).
3.3. Symbolic Representation of BNs. As a symbolic representation, a natural choice are Reduced Ordered Binary Decision Diagrams (ROBDD, or simply BDD) [Bry86], which can concisely encode Boolean functions or relations of Boolean vectors. Specifically, out implementation leverages the internal tools and libraries provided by the tool AEON [BBPŠ20]. Since a Boolean network consists of n Boolean variables and m Boolean inputs, any subset of V , C, or a relation X ⊆ V × C (a coloured set of vertices) can be seen as a Boolean formula over the network variables and inputs. That is, each network variable and logical input corresponds to one decision variable of the BDD. Here, a pair (s, c) belongs to such a relation iff it represents a satisfying assignment of this formula X. For relations of higher arity, fresh decision variables are created for each component of the relation. Standard set operations as described in Fig. 1 then correspond to logical operations on such formulae (∧ ≡ ∩, ∨ ≡ ∪, etc.).
Relation operations are similarly implementable using BDD primitives. In particular, existential quantification of a single decision variable (e.g. ∃s i .X or ∃c j .X) is a native operation on BDDs. Consequently, existential quantification on relations (as well as Colours) is simply a quantification over all decision variables encoding the specific relation component (i.e. all network variables for V , or all logical inputs for C). Finally, Swap only influences the way in which a BDD is interpreted -the actual structure of the BDD is unaffected.
To encode the network dynamics, notice that every update function b i can be directly represented as a separate BDD. From such BDDs, we can build one large BDD describing the whole coloured transition relation, which is traditionally used for the computation of Pre and Post. But the symbolic representation of such relation is often prohibitively complex for asynchronous systems. Instead, we compute Pre and Post using partial results for individual variables, which uses more symbolic operations but is less likely to cause a blow-up in the size of the BDD: Here, [s i → ¬s i ] is the standard substitution operation, which we use to flip the value of variable s i in the resulting formula if it does not agree with the output of b i . Note that this operation can be also implemented structurally directly on the BDD by exchanging the children of decision nodes conditioning on s i . Also note that sub-formulae that do not depend on X can be pre-computed once for the whole run of the algorithm, and the version of Pre and Post for monochromatic graphs can be implemented in exactly the same way.

Implementation
Finally, let us discuss a number of technical improvements which our algorithm employs in practice, and whose impact we consider in the evaluation section.
4.1. Pivot Selection. In Algorithm 1, we gave a naive implementation of the Pivots(X ) function. Here, we show how to implement it for BDDs in a much more concise way. Note that our approach uses the notation we established earlier for Boolean networks, but is generally applicable to any set or relation of bit-vectors represented using BDDs. First, notice that for a single network variable, we can define a similar operation, which we call Pick(i, X ): Here, we first restrict X only to the valuations which have s i = false, and then invert the value of s i (resulting in s i being always true in the set). Once we subtract these valuations from X , the resulting set then contains a valuation with s i = true only if it does not contain the same valuation with s i = false. Intuitively, for any valuation of the remaining BDD decision variables (i.e. s j and c j in our case) that is in X , we just picked a single unique value of s i (while preferring the value s i = false).
However, observe that we cannot simply apply Pick to every network variable alone to obtain the result of Pivots. Intuitively, the problem lies in the fact that Pick selects a witness for each variable in isolation, while Pivots considers all network variables as interconnected. We resolve this problem using a different equation, one which eliminates the picked variable in the recursive invocation: Pivots(X ) = F(X , s 1 , . . . , s n ) F(X , s 1 ) = Pick(1, X ) F(X , s 1 , . . . , s k ) = Pick(k, X ) ∩ F(∃s k .X , s 1 , . . . , s k−1 ) In this equation, the final case F(X , s 1 ) is clearly correct, since it simply defers to Pick(i, X ). However, to understand why the recursive case F(X , s 1 , . . . , s k ) is correct, observe the following: Assume that the set Y = F(∃s k .X , s 1 , . . . , s k−1 ) is computed correctly. That is, for any valuation of the remaining variables, Y contains a single unique incomplete witness valuation of variables s 1 , . . . , s k−1 . Now, since the BDD representing ∃s k .X does not depend on s k , each such unique witness must be included in Y twice: once with s k = true and once with s k = false. In other words, a single witness valuation of s 1 , . . . , s k−1 must be tied to two different valuations of the remaining variables, and these valuations are differentiated only by the variable s k . Now, one of these two valuations is necessarily included in the set Pick(k, X ). The other is either missing from X altogether, or is eliminated by Pick(k, X ). As such, computing Pick(k, X ) ∩ Y extends the witness from s 1 , . . . , s k−1 to s 1 , . . . , s k by eliminating one of the two aforementioned occurrences of the incomplete witness.
Observe that, as opposed to the original naive implementation of Pivots, this implementation only requires O(n) (i.e. O(log |V |)) symbolic operations in any case. 4.2. Saturation. In [CMS06], and later in greater detail within [ZC11], Ciardo et al. show that when the system is asynchronous, it may be much easier to compute reachable sets (and consequently SCCs) by applying only one transition (e.g. denoted t 1 ) at a time. Once applying t 1 cannot add new states to the reachable set, another transition (e.g denoted t 2 ) can be considered, respecting the order in which the affected variables appear in the symbolic data structure (Ciardo et al. employ multivalued decision diagrams, but the principle also applies to BDDs). If the application of other transitions causes that we can again add new states using t 1 , the process starts anew and t 1 is "saturated" again.
In the comparison presented in [ZC11], only the Xie-Beerel O(|V | 2 ) algorithm is used with saturation enabled, while the lock-step algorithm is used as given in [BGS00]. However, we argue that saturation can be also beneficial in the lock-step algorithm.
Asymptotic complexity. Unfortunately, combining lock-step with saturation disrupts the O(|V | · log |V |) asymptotic complexity of the algorithm. To see why this is the case, observe that classical symbolic reachability (i.e. a fixed-point algorithm iterating the Post procedure) requires O(|V |) steps to explore a graph. Meanwhile, a reachability procedure employing saturation needs O(|V ||T |) operations, where |T | is the number of distinct transitions. This is caused by the fact that saturation needs to check up to |T | transitions to discover a vertex. For example, consider an asynchronous graph employing transitions t 1 , . . . , t n such that t 1 and t n are alternated on a path of length O(|V |). Between considering t 1 and t n , saturation will attempt each of the |T | transitions, which are useless on this path, but still consume a symbolic operation.
Consequently, this |T | factor trickles down into the complexity of both the Xie-Beerel and lock-step algorithms if saturation is used, as both ultimately rely on some form of reachability to discover the graph vertices. The complexity of the coloured algorithms is then similarly affected.
Saturation and lock-step. The main idea of how saturation is applied in a coloured lockstep algorithm (for Boolean networks) is shown in Algorithm 2. The algorithm presents a helper function which performs one reachability step, similar to what is performed by the Post function. However, in this algorithm, only one transition is fired for each colour (we assume the iteration follows the order of variables as they appear in the symbolic representation, which benefits saturation). Additionally, a set R of colours that could not perform a step is computed. A similar procedure can be considered for backwards reachability, simply replacing VarPost with VarPre.
Note that there is a slight discrepancy between Algorithm 2 and the intuitive description of saturation that we gave earlier. In particular, we see that during a NextStep operation, a transition for each variable is triggered at most once, as opposed to the original description, where a transitions are fired repeatedly. This is caused by the simple nature of Boolean networks: In a BN, a single transition always modifies a single Boolean variable. Consequently, no new states can be discovered by firing a single transition multiple times in sequence. For other asynchronous systems, VarPost may need to be modify to apply the corresponding transition repeatedly. Additionally, note that we use the set R to ensure that VarPost (i.e. firing of a single transition) is executed only for colours for which we have not found a successor yet using some of the previously considered transitions. This is necessary to ensure that in each invocation of NextStep, each colour present in F is either advanced by one step (using exactly one transition), or is reported as converged within the returned set R.
Using this process, we can replace the Pre/Post procedures in the main lock-step algorithm (lines 11 and 12 of Algorithm 1). The R sets computed here are then used to update F lock and B lock (lines 13 and 14), as they exactly represent the converged colours that do not need further computation. A similar modification is necessary for the second while loop (lines 25-29), but here the sets of remaining colours R are not needed.
Algorithm 2: Main idea of the lock-step-saturation approach. The algorithm extends F with one additional reachability step, and returns a set of colours locked in this iteration (R). 4.3. Trimming and Parallelism. Most graphs typically contain a large number of trivial SCCs that introduce unnecessary overhead to the main algorithm. To avoid this overhead, we additionally perform a trimming step before each invocation of Decomposition. Trimming consists of repeatedly removing all vertices which have no outgoing or no incoming edges and is employed by most symbolic SCC algorithms on standard directed graphs as well.
The coloured analogue of trimming is straightforward, as it can be achieved using Pre and Post operations just as in the non-coloured case. For a coloured set of vertices V, operation Post(G, Pre(G, V) ∩ V) ∩ V returns only the vertices which have at least one predecessor in V. The successor variant simply exchanges the Post and Pre operations.
As such, applying this operation to each V until a fixed-point is reached before Decomposition is invoked eliminates the undesired trivial SCCs. Since the total number of steps performed collectively by all such fixed-point computations is bounded by |C||V | (the total number of removable vertex-colour pairs), this does not impact the overall asymptotic complexity of the algorithm.
In some cases, we have observed that the symbolic representation is able to handle the SCC computation but explodes during trimming. The algorithm then times-out during trimming, even though useful information about SCCs could be obtained if the trimming was skipped or postponed. To avoid this issue, we enforce an extra condition that a trimming procedure is terminated prematurely if the computed BDDs are more than twice the size (in terms of BDD decision nodes) of the initial set.
Additionally, the lock-step algorithm can be rather trivially parallelised. The recursive Decomposition calls operate on independent coloured vertex sets and can be therefore deferred to separate threads. Since the body of the Decomposition method is rather complex, this can be done easily with a queue guarded by a mutex which is shared between all threads (i.e. the synchronisation overhead is negligible due to the long running time of Decomposition). Finally, a simple termination detection procedure is needed to ensure that idle threads do not terminate prematurely while decomposition is still running.
Note that most BDD packages are not internally thread-safe, as they share decision node memory across different BDD objects. In our experiments, this aspect is handled by cloning the set V corresponding to each recursive invocation, plus the symbolic representation of the BN necessary to compute Post and Pre. As such, the memory used to represent BDDs manipulated by each thread is completely independent from other threads.

Experimental Evaluation
To test the algorithm, we compiled a benchmark set of Boolean networks from the CellCollective [HKM + 12] and GINsim [CNT12] model databases. Since the models in these databases contain fully specified networks, uninterpreted functions were introduced into existing models by pseudo-randomly erasing parts of the existing update functions.
While this process is to some extent artificial, we believe it to be a good approximation of the model development process, where at some point, the structure of the network is already established, but its dynamics are still not fully determined. Using this process, we obtained a collection of networks ranging between 2 20 and 2 50 in the size of the coloured graph (i.e. |V × C|). Note that for each graph, we consider only a subset of possible input valuations that is biologically relevant with respect to the established network structure. For example, the first model (i.e. [SOHMMA17]) admits 2 48 input valuations, but only 2 19 are biologically relevant due to constraints on function monotonicity.
A complete overview of the employed models is given in Table 1. For each model, we give the number of discovered non-trivial components as an interval, because each colour can correspond to a different number of components. We employ a 24h timeout for all experiments.
The experiments were performed on a 32-core AMD Threadripper workstation with 64GB of RAM memory. All tested models are available in our source code repository. 3 Note that the smaller models (< 2 30 ) should be easy to process even on a less powerful machine; however, the larger models can require substantial amount of memory.
For each model, we have tested the lock-step algorithm as presented in the main part of this paper (Lock-step in Table 2), an enhanced version with saturation enabled (Satur. in Table 2), and a parallel implementation which also includes saturation (Parallel in Table 2). In all algorithms, we employ the trimming optimisation.
From the results, we can see that parallelisation improves the performance of the algorithm significantly: in case of models with a large number of SCCs, we see an up-to 30x speed-up, comparing Parallel and Satur. in Table 2. On the other hand, when the number  of SCCs is small (such as [OLB + 08]), the speed-up is understandably minimal, since the number of independent recursive calls is also small. As expected, the total number of SCCs has a significant impact on the performance of the algorithm (e.g. [Iro09] and [CMR + 15]) overall, since the number of calls to Decomposition increases. Furthermore, we see that our "coloured saturation" indeed provides a performance benefit. However, this improvement is mostly incremental.
After further analysis, we discovered that the whole algorithm is often limited by the performance of the trimming procedure, rather than reachability procedures though. In particular, the use of saturation has significantly reduced the size of symbolic representation during computation of reachability, however the symbolic representation still performs rather poorly (at least for Boolean networks) during trimming. This limits the performance of the whole method, since all the considered graphs contain a large portion of trivial SCCs. Furthermore, in many cases the number of iterations needed to completely trim a set of states is substantial. This leads us to believe there is still space for improvement in terms of SCC detection in large Boolean networks, even without parameters.
Finally, we examined the benefit of processing all colours simultaneously versus a naive parameter scan approach, where each monochromatic case is handled separately. To do so, we considered various pseudo-random monochromatisations of the studied models and processed these using our algorithm. Here, we observe that for the four models with at least 20 variables, no computation for any of the monochromatic models finished in under one second (with T-cell differentiation typically requiring more than one minute due to the relatively large number of components).
Consequently, we can extrapolate that computing the full coloured SCC decomposition using such naive parameter scan would require more than 10 ours for each model (and 10+ days in the case of T-cell differentiation). This approach could be to some extent beneficial in a massively parallel environment (hundreds or thousands of CPUs), but the coloured approach clearly scales better in setups where resources are more limited.

Conclusions
This paper presents a fully symbolic algorithm for detecting all monochromatic strongly connected components in edge-coloured graphs. The work has been motivated by systems sciences, namely systems biology, where the need for efficient automated analysis of components in large graphs with a large sets of coloured edges is emerging. The algorithm combines several ideas inspired by existing state-of-the-art algorithms for SCC decomposition in a non-trivial way. We believe this is the first fully symbolic algorithm aiming to solve the problem efficiently.
The experimental evaluation has shown that the algorithm can handle large, real-world systems that would be otherwise too large to fit into the memory of a conventional workstation (> 2 32 ), and that the performance of the algorithm can be further improved using saturation and parallelisation. Finally, the algorithm has a strong potential to be significantly faster than iterating a standard algorithm for SCC decomposition executed on all monochromatic sub-graphs one-by-one.