 Methodology article
 Open Access
 Published:
Identifying dynamical modules from genetic regulatory systems: applications to the segment polarity network
BMC Bioinformaticsvolume 8, Article number: 413 (2007)
Abstract
Background
It is widely accepted that genetic regulatory systems are 'modular', in that the whole system is made up of smaller 'subsystems' corresponding to specific biological functions. Most attempts to identify modules in genetic regulatory systems have relied on the topology of the underlying network. However, it is the temporal activity (dynamics) of genes and proteins that corresponds to biological functions, and hence it is dynamics that we focus on here for identifying subsystems.
Results
Using Boolean network models as an exemplar, we present a new technique to identify subsystems, based on their dynamical properties. The main part of the method depends only on the stable dynamics (attractors) of the system, thus requiring no prior knowledge of the underlying network. However, knowledge of the logical relationships between the network components can be used to describe how each subsystem is regulated. To demonstrate its applicability to genetic regulatory systems, we apply the method to a model of the Drosophila segment polarity network, providing a detailed breakdown of the system.
Conclusion
We have designed a technique for decomposing any set of discretestate, discretetime attractors into subsystems. Having a suitable mathematical model also allows us to describe how each subsystem is regulated and how robust each subsystem is against perturbations. However, since the subsystems are found directly from the attractors, a mathematical model or underlying network topology is not necessarily required to identify them, potentially allowing the method to be applied directly to experimental expression data.
Background
Genetic regulatory systems are assumed to be 'modular', with specific combinations of genes and proteins responsible for different biological functions. Although there is no authoritative definition of a 'module', one common description is a group of genes, proteins and/or molecules that combine to carry out a relatively distinct function (distinguishable from the functions associated with other modules) [1]. Rather than being a protein complex or group of coexpressed genes, such a module can be viewed as the temporal activity (dynamics) of a group of genes/proteins that controls a specific function in different environmental conditions, cell types and/or tissues. For example, the temporal activity of genes/proteins controlling progression through the cell cycle (in many different environmental conditions and cell types) can be viewed as a module. Moreover, since genetic regulatory systems are hierarchical [2], it is likely that these modules also interact in a hierarchical manner.
It is this concept of a 'dynamical module' or 'subsystem' that we consider in this paper. An important challenge in biology is how to identify such subsystems and discover how they combine hierarchically to determine cell types, tissues and organisms. A dynamic approach is the most suitable starting point for decomposing a genetic regulatory system into modules, since it is the activity profiles that control biological functions. The topology of the underlying interaction network is just a description of the interactions associated with the activity profiles.
Most previous attempts at modular decomposition have relied on the topology of the underlying network, deduced from proteinprotein interaction and/or transcription factor binding data. From a generic network point of view, algorithms have been created which partition any network into topological modules [3]. From a biological viewpoint, transcription factor networks have been used to find 'network motifs' [4] and 'structural modules' [5], based solely on topological properties of the network. Although informative, topology based approaches do not necessarily allow the functions/dynamics associated with groups of genes/proteins to be inferred [6]; rather, activity/expression levels over a period of time are required. When expression data have been used, the resulting modules tend to be groups of genes whose expression levels change in unison. For example, 'transcription factor modules' [7, 8] are groups of genes with common transcription factor binding sites and expression patterns. However, for a biological function, it may be the case that a series of interactions are involved, with different genes being affected at different points in time (e.g. genes involved in cell cycle regulation). Similar issues arise when applying clustering algorithms to expression data, which also groups together genes whose expression levels change in unison. In this paper, we present a method for identifying subsystems ('dynamical modules'), given a set of discrete state, discrete time attractors. The subsystems are found directly from the attractors (Fig. 1), without the need for topological information, making the method applicable both to models and (potentially) to experimental expression data. To demonstrate the method, we consider a class of mathematical models called Boolean network models (defined in the Methods section). These models are used as a starting point because of their relative simplicity and the fact that different dynamics can be easily compared. We also consider how each subsystem is regulated and how subsystems can help investigate robustness in these systems. In order to test our method, we then apply it to a model of the Drosophila segment polarity network.
Results: Identifying subsystems and regulatory interactions between subsystems
Given a set of discretestate, discretetime attractors as an input, we have produced a method that optimally breaks the attractors up and identifies subsystems, each one

1.
Conserved across some set of attractors,

2.
Distinguishable from the stable dynamics for the rest of the system.
In order to demonstrate the key features of the resulting subsystems, we use the simple example network in Fig. 2. Here, A_{1}, ..., A_{4} are four attractors for the system, but it is evident that the subdynamics associated with some subnetworks (S_{1}, ..., S_{6}) can distinguish different sets of attractors and highlight parts of the system that are dynamically/functionally distinct. For example, the subsystem S_{1} is conserved across three attractors (A_{1}, A_{2} and A_{3}) and provides a way of distinguishing those attractors from the remaining attractors (A_{4}). Moreover, S_{1} can be viewed as dynamically distinct from all remaining subdynamics, since the dynamics of each remaining node (1, 4, 5 and 6) has a distinguishable profile when viewed alongside S_{1}. In particular,
(a) Nodes 1, 5 and 6 each exhibit different behaviours/dynamics alongside S_{1} (compare A_{3} with A_{1} and A_{2}).
(b) The activity of node 4 oscillates out of phase with the activity of nodes 5 and 6, when viewed alongside S_{1} in attractors A_{1} and A_{2}.
The subdynamics S_{1}, ..., S_{6} in Fig. 2 are all subsystems (according to our method) because each one can be viewed as dynamically distinct from all the other subdynamics in attractors A_{1} – A_{4}. The attractors themselves than can be viewed as combinations of subsystems.
The main method does not explicitly identify functional/regulatory relationships between subsystems. However, in the case of logical models such as Boolean network models, the logical functions can be used to identify suitable regulatory relationships.
For the remainder of this section we describe these new techniques in terms of Boolean network models (which we formally define in the Methods section). Algorithms for these new techniques are described in the Methods section of this paper, whilst formal proofs for all the methods can be found in Additional files 1 and 2. It should be noted that all of the definitions and algorithms can be extended to any set of discrete state, discrete time attractors, with or without a detailed mathematical model.
Partial state sequences
Much of this paper involves looking at the subdynamics in Boolean network models. Therefore, for a subset of nodes N ⊆ V, we consider partial states.
Definition 1. A partial state, x^{N}∈ {0, 1}^{N}is a set of Boolean states, one for each node n_{ i }∈ N ⊆ V. i.e. x^{N}= {s_{ i }: n_{ i }∈ N}.
Definition 2. A partial state sequence, $P=\{{x}_{0}^{N},{x}_{1}^{N},\mathrm{...},{x}_{q1}^{N}\}$, is an ordered set of partial states, for a node set N (⊆ V).
These partial state sequences can be used to represent subdynamics in Boolean network models, and therefore, can be used as a starting point when defining subsystems. In particular, we are looking for partial state sequences that occur (or cycle within) attractors.
Definition 3. A partial state sequence $P=\{{x}_{0}^{N},{x}_{1}^{N},\mathrm{...},{x}_{q1}^{N}\}$occurs in an attractor A = {z_{0}, ..., z_{p1}} if there exists integers b_{0}, ..., b_{p1}∈ {0, ..., q  1} for which the following are true

1.
For k = 0, ..., p  1, ${x}_{{b}_{k}}^{N}$ = {s_{ i }∈ z_{ k }: n_{ i }∈ N}.

2.
For each k ∈ {0,..., p  1} and j = k  1 (mod p), either
(a) b_{ k }= b_{ j }or (b) b_{ k }= b_{ j }+ 1 (mod q)

3.
Properties 1 and 2 are not true for any smaller partial state sequence ${P}^{\prime}=\{{y}_{0}^{N},{y}_{1}^{N},\mathrm{...},{y}_{{q}^{\prime}1}^{N}\}$ and integers c_{0}, ..., c_{p1}∈ {0, ..., q'  1} (q' <q).
As an example of a partial state sequence that occurs in an attractor, consider the partial state sequence
and the attractor A_{1} = {z_{0}, z_{1}, z_{2}, z_{3}} in Fig. 2. As the system enters the attractor A_{1}, we continually cycle through the component states over time (i.e. z_{0}, z_{1}, z_{2}, z_{3}, z_{0}, z_{1}, z_{2}, z_{3}, z_{0}, ...). Therefore, since

$${x}_{0}^{N}$$
is contained in z_{0} and z_{1} (i.e. ${x}_{0}^{N}$ = {s_{ i }∈ z_{0} : n_{ i }∈ N} = {s_{ i }∈ z_{1} : n_{ i }∈ N})

$${x}_{1}^{N}$$
is contained in z_{2} and z_{3} (i.e. ${x}_{1}^{N}$ = {s_{ i }∈ z_{2} : n_{ i }∈ N} = {s_{ i }∈ z_{3} : n_{ i }∈ N})
we also continually cycle through the partial states of P over time (i.e. ${x}_{0}^{N},{x}_{1}^{N},{x}_{0}^{N},{x}_{1}^{N},{x}_{0}^{N}$, ...). We note that the time taken to change from one partial state to the next is not considered (2 time steps in this example). Ignoring such time lags could be especially important in genetic regulatory systems, where the same process may take different lengths of time under different conditions (in different attractors). For example, in cell cycle regulation, some mutations can alter the time taken for a round of cell division, without necessarily altering the relationships that exist between other groups of genes/proteins [9]. Properties 1 and 2 allow us to capture the fact that P continually cycles within the the attractor A_{1} by noting that a sequence of integers {b_{0}, b_{1}, b_{2}, b_{3}} = {0, 0, 1, 1} exists that maps the partial states in P onto the attractor states in A_{1}. Property 2 ensures that the partial states in P occur in the correct order, whilst allowing different lengths of time between each change. Properties 2 and 3 ensure that P is the smallest set of partial states that occurs in the attractor A. This leaves a partial state sequence that just describes the 'order' in which the node states change in A (for nodes in N), ignoring any time lags associated with individual partial states and the number of times P cycles within A. Returning to the example in Fig. 2, S_{3} cycles twice in A_{1} and A_{2}. However, just two partial states (${x}_{0}^{N}$ = {s_{4} = 1} and ${x}_{1}^{N}$ = {s_{4} = 0}) are sufficient to capture the fact that S_{3} occurs in A_{1} and A_{2} (letting {b_{0}, b_{1}, b_{2}, b_{3}} = {0, 1, 0, 1} and {1, 0, 1, 0} respectively).
Given a node set N and a set of attractors, the Methods section of this paper gives detailed algorithms of how to identify the unique partial state sequences (within an order of rotation) that occur in each attractor. Essentially this is done by going through each attractor state and writing down the partial states involving the node set N. We then remove partial states when
(a) Multiple adjacent partial states are identical (in which case only one copy is kept)
(b) A sequence of states cycles many times within an attractor (in which case only one copy is kept).
leaving a partial state sequence that just describes the 'order' in which the node states change in that attractor.
Method for identifying subsystems
The new method of identifying subsystems is a two stage process that breaks up the system's attractors, to leave partial state sequences that are optimally distinguishable from one another. A summary of these two stages, along with relevant definitions, is given below using the simple example in Fig. 2 and Fig. 3. Detailed descriptions of the algorithms are given in the Methods section for those wanting to implement this method. A flow diagram of the procedures involved can be seen in Fig. 1
Stage 1: Partition sequences
The first stage of the method involves identifying every partial state sequence that satisfies Definition 5. Each such partition sequence (for a node set N) is both (a) conserved across a set of attractors and (b) has a fundamentally different dynamical profile from that associated with any larger node set M ⊇ N. We first introduce intersection sequences that are used in the definition of partition sequences and set up a hierarchy amongst them.
Definition 4. A partial state sequence P for a node set N is an intersection sequence, if there exists a subset of attractors $\u2102(\subseteq \mathbb{A})$ for which the following hold

1.
P occurs in every attractor A ∈ $\u2102$

2.
P does not occur in any attractor A ∉ $\u2102$

3.
Given a larger node set M ⊃ N, there is no partial state sequence P' (for the node set M) that satisfies condition 1.
If the above properties hold, we say P intersects at $\u2102$.
In order to identify every intersection sequence, node sets (N) are analysed to identify which partial state sequences occur in which attractors. Any such partial state sequence (P) and corresponding set of attractors ($\u2102$) will then satisfy properties 1 and 2 of Definition 4. These pairs {P, $\u2102$} are stored for each node set N and then compared during the procedure to remove those that fail property 3. In the Methods section of this paper we give a detailed algorithm for this as well as providing some notes on improving efficiency. In practice, many node sets N and/or attractors can be ignored because we know that the pairs {P, $\u2102$} will fail property 3.
The set of all intersection sequences provide a hierarchical breakdown of the system's dynamics (see Fig. 3A, B and Table 1). At the top of the hierarchy are partial state sequences that contain only a few nodes but are conserved over (relatively) many attractors. Then as you go down the hierarchy, extra nodes are added but the resulting partial state sequences are conserved over fewer attractors. Property 3 above ensures that each partial state sequence is optimal in the sense that no more extra nodes can be considered without the new extended partial state sequence occurring in strictly fewer attractors.
Although these sequences provide a neat hierarchical breakdown of the system's dynamics, some may be superfluous whilst important 'core' subdynamics may be missed. Therefore, we introduce 3 extra constraints when defining partition sequences.
Definition 5. A partial state sequence $P=\{{x}_{0}^{N},{x}_{1}^{N},\mathrm{...},{x}_{q1}^{N}\}$ is a partition sequence if it satisfies any of the following properties, for some set of attractor cycles $\u2102(\subseteq \mathbb{A})$
A : P is Core to $\u2102$
The following 3 properties hold for P
1.P occurs in an intersection sequence P', which intersects at $\u2102$ (P can equal P').

2.
If an intersection sequence Q (for a node set M) intersects at $\mathbb{D}$ (where $\mathbb{D}\cap \u2102\ne \varnothing $), then there exists an intersection sequence Q' (for a node set M' ⊇ M ∪ N) that occurs in every attractor A ∈ $\mathbb{D}\cap \u2102$

3.
1 and 2 are not true for any larger partial state sequence P" (for a node set N" ⊃ N)
B : P is Exclusive to $\u2102$
P is the only intersection sequence that intersects at $\u2102$.
C : P is Independently Oscillating
P intersects at $\u2102$ and cycles out of phase with another intersection sequence Q. i.e. ∃ Q that involves the node set M and intersects at $\mathbb{D}$, for which

1.
$\u2102\cap \mathbb{D}$ ≥ 2

2.
N ∪ M = V (the set of all nodes)
The three parts of this definition are used to describe three ways in which a subdynamic (partial state sequence) can be viewed as fundamental to the make up of the attractors and/or distinguishable from the remaining subdynamics. For the example in Fig. 2 and Fig. 3, Table 2 shows the partial state sequences that satisfy properties A, B and C above.
Two partition sequences P_{1} and P_{2} are not intersection sequences, but are core to the dynamics associated with different sets of attractors. For example, P_{2} is the core dynamic spanning attractors A_{1} and A_{2}, rather than P_{5} and P_{6}, which both occur in A_{1} and A_{2}. This is because P_{5} and P_{6} involve partially overlapping node sets N_{5} and N_{6} (i.e. N_{5} ⊈ N_{6}, N_{6} ⊈ N_{5}), and so neither corresponds to a subnetwork whose dynamics are core to attractors A_{1} and A_{2}. On the other hand P_{2} acts as a central point from which these intersections sequences branch off.
Although in this example every intersection sequence is also a partition sequences, this does not necessarily have to be the case. One example where this is not the case is given in Section S3.2 of Additional file 3. In order to identify every partition sequence, we use the full set of intersection sequences to identify every partial state sequence satisfying property A, B or C of Definition 5. These 3 tests are done independently and the partial state sequences found in each test are grouped together to give the full set of partition sequences. In the Methods section of this paper we give detailed algorithms for these procedures.
Stage 2: Subsystems
The second stage of the method compares the partition sequences, to identify every partial state sequence that satisfies Definition 6. The partition sequences allow the attractors to be broken up in a hierarchical manner. However, these sequences are not subsystems (for a start, all of the attractors themselves are partition sequences). Of more interest are the key subdynamics that are unique to a particular partition sequence and that distinguish one level in the hierarchy from the next. Therefore, we define subsystems to be those components that are unique to a partition sequence. i.e.
Definition 6. A partial state sequence S (for a node set N) is a subsystem if it is unique to a partition sequence. i.e. there exists a partition sequence P (for a node set M) for which

1.
S occurs in P.

2.
If another partition sequence P' (for a node set M' ⊂ M) occurs in P, then M' ∩ N = ∅.

3.
1 and 2 are not true for any partial state sequence S', for a larger node set N' ⊃ N.
In the Methods section of this paper we give detailed algorithm for identifying every such subsystem. This essentially involves running through every partition sequence P (for a node set M) and identifying the nodes (N) that do not belong to any other partition sequence P', that occurs in P. The dynamics associated with this node set N are then a subsystem.
As an example, consider P_{5} in Fig. 3. Looking at the remaining partition sequences we see that only P_{1} and P_{2} occur within it. Therefore, the nodes and dynamics in P_{5} that are unique to P_{5} are the node set N = {n_{5}, n_{6}} and the subsystem S_{5} (from Fig. 2D). Obviously, S_{5} satisfies Property 1 of Definition 6 because it occurs in the partition sequence P_{5}. Property 2 is satisfied since P_{1} and P_{2} are the only partition sequences occurring in P_{5} and the node sets do not intersect with N = {n_{5}, n_{6}} ({n_{2}, n_{3}} ∩ {n_{5}, n_{6}} = ∅ and {n_{1}, n_{2}, n_{3}} ∩ {n_{5}, n_{6}} = ∅). Property 3 is satisfied since any node set N' ⊆ M, N' ⊃ N would fail Property 2. For the example Boolean network model, all the subsystems can be seen in Fig. 2, whilst Fig. 3C and Table 3 demonstrate which subsystem corresponds to which partition sequences.
Regulation of subsystems
We say a collection of subsystems $\mathbb{S}$_{ x }= {S_{1}, ..., S_{ f }} triggers an individual subsystem S_{ y }in an attractor A, if the cooccurrence of subsystems S_{1}, ..., S_{ f }ensures the occurrence of S_{ y }in A.
For Boolean network models (or other discretestate models), such interactions can be identified by considering the Boolean functions (or Logical functions). In the case of the simple example in Fig. 2, it is possible to look at the Boolean functions and say that
(a) S_{2} can trigger the occurrence of S_{1} in attractors A_{1} and A_{2}.
(b) S_{1} can trigger the occurrence of itself (S_{1}) in attractors A_{1}, A_{2} and A_{3}.
This is since
(a) The prolonged activation of node 1 (s_{1} = 1) is sufficient to activate nodes 2 and 3 (s_{2} = 1, s_{3} = 1).
(b) If nodes 2 and 3 are both on at time t, then this sufficient to maintain their activation for all time steps t' ≥ t.
As can be seen above, different collections may act in different sets of attractors to fully explain the occurrence of a subsystem S_{ y }.
Definition 7. A set of subsystem collections ${\mathbb{S}}_{1},\mathrm{...},{\mathbb{S}}_{g}$regulates an (individual) subsystem S_{ y }if the following are true

1.
For i = 1, ..., g, ∃ an attractor A for which $\mathbb{S}$_{ i }triggers S_{ y }in A.

2.
If S_{ y }occurs in an attractor A, ∃ i ∈ {1, .., g} for which $\mathbb{S}$_{ i }triggers S_{ y }in A.
We call the set {${\mathbb{S}}_{1},\mathrm{...},{\mathbb{S}}_{g}$} the regulation set of S_{ y }.
In the Methods section, we describes how a suitable regulation set can be found for each subsystem S_{ y }, by looking at how each partial state within S_{ y }is triggered/regulated in each attractor. This approach uses the Boolean functions from a Boolean network model to identify which partial states can trigger the occurrence of those in S_{ y }.
For the simple example in Fig. 2, Table 4 shows a regulation set for each of the subsystems. As can be seen in this example it is often the case that subsystems play a role in their own regulation, although this does not necessarily have to be the case.
Even without the Boolean functions, it is possible to identify relationships between subsystems. On a simple observational level, a subsystem S_{ x }may be hierarchically linked to another subsystem S_{ y }, because S_{ x }only occurs in an attractor in conjunction with the 'higher order' S_{ y }. i.e.
Definition 8. Consider two subsystems S_{ x }and S_{ y }. S_{ x }is hierarchically linked to S_{ y }if the following are true

S_{ x }occurs in an attractor A ⇒ S_{ y }occurs in an attractor A
Furthermore, such a link can be viewed as direct if it is impossible to find a subsystem S_{ z }for which the following is true

1.
S_{ x }is hierarchically linked to S_{ z }

2.
S_{ z }is hierarchically linked to S_{ y }

3.
There exists attractors A_{1} and A_{2} for which
(a) S_{ y }occurs in A_{1} and A_{2}
(b) S_{ z }occurs in A_{1} but not A_{2}
(c) S_{ x }occurs in neither A_{1} nor A_{2}
This terminology can easily be extended to collections of subsystems.
Robustness of attractors and subsystems
Suppose S is a subsystem that involves a set of nodes N and occurs in a set of attractors $\u2102$. For any network state z_{ k }= {s_{1}, ..., s_{ v }} in any attractor A ∈ $\u2102$, a node n_{ f }can have its state s_{ f }perturbed to (s_{ f }+ 1) (mod 2), to give a new network state ${z}_{k}^{(f)}$ (say). By looking at how often such a network state ${z}_{k}^{(f)}$ converges to an attractor A' ∈ $\u2102$ containing the same subsystem, it is possible to measure how robust the subsystem S is to perturbations in the system.
Let r(S, L) be the probability that perturbing the state of a node n_{ f }∈ L in an attractor A ∈ $\u2102$ causes the system to converge to an attractor A' ∈ $\u2102$. i.e.
where
Then, using this measure, it is possible to measure Global Robustness = r(S, V), Internal Robustness = r(S, N), and External Robustness = r(S, M) (where V is the set of all nodes and M = V\N). This allows one to distinguish the effects of perturbing nodes within and outside of the subsystem. The measure r(S, L) can also be adapted to give the robustness of an individual attractor A (by letting S = A, $\u2102$ = {A} and L = V).
Results: Application to the segment polarity network
During the development of the fruit fly Drosophila melanogaster, the embryo becomes segmented, with the segment polarity genes responsible for specifying the number, spacing and polarity (direction) of these segments. In order to demonstrate the applicability of our technique to genetic regulatory systems, we analyse a previously published Boolean network model of the Drosophila segment polarity network from references [10] and [11]. This model corresponds to a 4cell ring, where intercellular interactions are allowed between cells 1 and 4. This approach was kept here since the segments start off 4 cells wide and the wild type patterns repeat every 4 cells (before cell proliferation). The Boolean functions for this model are shown in Table S4.1 of Additional file 4, whilst the main interactions are summarised in Fig. 4. This 4 cell model has 10 fixed point attractors. The wild type attractor A_{1} is shown in Fig. 5, and is characterised by WG and EN/HH expression in 1 cell wide stripes either side of the parasegment boundary (in cell 4 and 1 resp). Two other attractors, A_{2} and A_{3}, correspond to experimentally observed phenotypes and are also shown in Fig. 5. All 10 attractors are shown in Additional file 4 (Fig. S4.1).
Applying our method, we identified 19 subsystems that satisfied Definition 6, which are shown in Fig. 6 and Table 5. S_{ A }, which occurs in all 10 attractors, corresponds to the cellular response to SLP (expressed in cells 3 and 4 only). Of the remaining subsystems, S_{B 1}, S_{B 2}, S_{C 1}and S_{C 2}appear to be the most interesting since they capture a large proportion of the global dynamics. They correspond to different states and different cells associated with the same positive feedback loop, involving Wingless (WG), Engrailed (EN), Hedgehog (HH) and Cubitus interruptus (CIA). The feedback loop can either be ON (S_{B 1}, S_{C 1}) or OFF (S_{B 2}, S_{C 2}) and either involve intercellular interactions between cells 4 and 1 (S_{B 1}, S_{B 2}) or cells 2 and 3 (S_{C 1}, S_{C 2}). As mentioned earlier, the wild type attractor is characterised by WG and EN/HH expression in 1 cell wide stripes either side of the parasegment boundary (in cell 4 and 1 respectively). This is captured by S_{B 1}.
A previous study has considered the modular design of this system [12], but there are significant differences in the component 'modules'. In [12], modules were arbitrarily chosen to be important intracellular pathways involving (a) SLP, (b) EN/HH and (c) WG/CI. However, here we have found that the most important subsystems are intercellular combinations of these pathways. Firstly, (a) and (b) combine to form S_{ A }. Secondly, (b) and (c) combine across cell boundaries to give 2cell positive feedback loops S_{B 1}, S_{B 2}, S_{C 1}and S_{C 2}. Moreover, the method in the paper captures the states of these pathways.
Interactions exist between these 19 subsystems, which explain how each one is regulated and which of the 10 attractors they occur in. For each of the 19 subsystems, Table 6 shows the sets of subsystems that regulate them, by ensuring their occurrence in an attractor (see previous section). Of the 19, the five subsystems highlighted in Fig. 6 play the largest role in regulating other subsystems with S_{ A }, S_{B 1}, S_{B 2}, S_{C 1}and S_{C 2}involved in the regulation of 12, 9, 9, 9 and 9 subsystems (respectively). Using these five main subsystems as a starting point, it is possible to generate a hierarchical breakdown of the 10 attractors in the model (Fig. 7A). The 10 attractors all contain S_{ A }, then split up into 4 main groups depending on the occurrence of S_{B 1}, S_{B 2}, S_{C 1}and S_{C 2}. In two cases, where $\mathbb{S}$_{1} = {S_{B 1}, S_{C 1}} or $\mathbb{S}$_{2} = {S_{B 2}, S_{C 2}} occur, this is sufficient to distinguish the attractors A_{2} and A_{3} (respectively) from the others. In the remaining two cases, where $\mathbb{S}$_{3} = {S_{B 1}, S_{C 2}} and $\mathbb{S}$_{4} = {S_{B 2}, S_{C 1}} occur, two additional levels of specification determine the individual attractors. Additional file 4 provides further details of the interactions involved in specifying each individual attractor (and ensuring ${\mathbb{S}}_{1},\mathrm{...},{\mathbb{S}}_{4}$ occur in an attractor).
Since only 3 of the 10 attractors from the model correspond to observed phenotypes, it is important to ensure that the main subsystems aren't an unrealistic artefact of the model. Reapplying our method of identifying subsystems to the reduced set of observed attractors {A_{1}, A_{2}, A_{3}}, we find that the 5 main subsystems in Fig. 6 are preserved. The only change is that the new subsystems incorporated some of the old smaller subsystems. In particular we get the following 5 main subsystems

$${S}_{A}^{\ast}$$
= S_{ A }∪ S_{D 2}∪ S_{E 2}∪ S_{H 2}

$${S}_{B1}^{\ast}$$
= S_{B 1}∪ S_{F 2}∪ S_{J 2}

$${S}_{B2}^{\ast}$$
= S_{B 2}∪ S_{F 1}∪ S_{J 1}

$${S}_{C1}^{\ast}$$
= S_{C 1}∪ S_{G 2}

$${S}_{C2}^{\ast}$$
= S_{C 2}∪ S_{G 1}
Moreover, as can be seen in Fig. 7B, the 5 main subsystems still provide a breakdown of the attractors. Focusing on subsystems allows us to look at the segment polarity network from a different angle. The first thing that stands out is that the subsystems exhibit a symmetry across the cell 1/2 boundary, whereby every subsystem in cell 4/1 has a symmetric counterpart in cell 3/2 (compare the pairs S_{ A }S_{ A }, S_{B 1}S_{C 1}, S_{B 2}S_{C 2}, S_{D 1}S_{E 1}, S_{D 2}S_{E 2}, S_{F 1}S_{G 1}, S_{F 2}S_{G 2}, S_{H 1}S_{I 1}, S_{H 2}S_{I 2}, S_{J 1}S_{J 1}and S_{J 2}S_{J 2}). Moreover, as can be seen in Table 6, these symmetric counterparts are analogously regulated. The protein SLP is crucial to setting up this symmetry, since it sets off a chain of interactions that nullifies hh, HH and en in cells 3 and 4 (see S_{ A }in Fig. 6). This then blocks signalling between cells 3 and 4, since (1) HH and WG are the only proteins involved in cell – cell signalling (in this model) and (2) en is the only gene (in this model) that can receive signals from WG. Therefore, it appears that one of the effects of SLP is to partition the cells in the embryo into 'isolated' 4 cell wide blocks (at this developmental stage) with the cell 1/2 boundary at the centre of this block. SLP also imposes restrictions on interactions within the neighbouring cells 1 and 2. This appears to be sufficient to partition the 4 cells into two, 2cell blocks (cell 4/1 and cell 2/3) that have relatively independent subdynamics. These two, 2cell blocks are forced to choose dynamics from S_{B 1}, S_{B 2}, S_{C 1}and S_{C 2}, which in turn specify the attractor chosen. Important interactions do occur across the 1/2 cell boundary and subdynamics are conserved across this boundary. However, the driving force behind the dynamics and specification of the network appear to come from the two, 2cell wide blocks. Whether this symmetry is inherent in the system or an artefact of the model remains to be verified, since only three attractors have been observed experimentally, and so more data may be required. In order to gain added insight into both the subsystems and the attractors, we calculate the robustness of each attractor and subsystem (see Tables 7 and 8). The robustness score is between 0 and 1 and measures how often the attractor/subsystem can survive after a perturbation to any node state. S_{B 1}and S_{C 1}are maximally robust in that no single node perturbations can destroy them. This also implies that S_{B 1}and S_{C 1}draw in local subdynamics that only marginally differ from them. Therefore, since S_{B 1}and S_{C 1}both occur in A_{2}, this partially explains why A_{2} is so dominant in the state space, with over 98 % of network states converging to it. On the other hand S_{B 2}and S_{C 2}are vulnerable to perturbations, especially to nodes within the subsystems themselves. The wild type attractor is A_{1}, which contains S_{B 1}and S_{C 2}. Only a small proportion of the state space (0.01 %) converges to this attractor but it is still relatively robust (0.67). It appears that this robustness is primarily due to S_{B 1}, the same subsystem responsible for the characteristic WG and EN/HH expression either side of the parasegment boundary (in cell 4 and 1 resp).
Discussion and Conclusions
Described in this paper is a framework for identifying subsystems ('dynamical modules') from a Boolean network model. This method of identifying subsystems is applicable for systems with either fixed point attractors, cyclic attractors or both. The methods are designed to be applicable to genetic regulatory systems. Therefore, we have applied the method to an existing model of the Drosophila segment polarity network. We identify novel subsystems acting across cell boundaries, and demonstrate how they regulate each other to give the attractors. Analysis of the subsystems also allows us to predict which ones underlie robustness in the wild type attractor.
Our methods are ideally suited to analysing multistable systems, where different stable states/conditions are captured within different attractors. For a node set with two or more different dynamics in multiple attractors, the different dynamics will form part of different intersection sequences, partition sequences and subsystems. Therefore, some identified subsystems will capture the key components associated with particular stable states (for some set of nodes). Meanwhile, other subsystems will capture shared/conserved dynamics. For example, in the Drosophila segment polarity network model, subsystems S_{B 1}, S_{B 2}, S_{C 1}and S_{C 2}(Fig. 6) capture different states of 2 positive feedback loops, which underlie multistationarity in this particular system. On the other hand, the subsystem S_{ A }captures a conserved substate that can cooccur alongside all of the remaining subsystems.
Once the subsystems have been identified, they can provide novel insight into aspects of the model. Firstly, Boolean functions from the model can be used to determine how each subsystem is regulated. This knowledge can then be used to provide a more detailed hierarchical breakdown of the original attractors, allowing the key differences and similarities between different attractors to be highlighted. A second advantage of looking at subsystems is that they can provide a new way of looking at robustness and adaptability. An external perturbation may switch the system from one attractor to another. However, despite this, a subsystem may be robust and remain unchanged. Looking at the effect of external perturbations on subsystems can highlight which parts of the system are most robust, or how certain subsystems can adapt/change without affecting others.
Although we have focussed on Boolean network models in this paper, the method can be adapted to other discretestate, discretetime models/systems. Moreover, since subsystems (in our method) are found directly from the attractors, a mathematical model or underlying network topology is not even necessarily required to identify them. Hierarchical links (see Results) between these subsystem can also be found without a mathematical model. These links allow dependencies between the subsystems to be discovered, providing insight into potential functional links. However, to gain a full understanding of how the subsystems interact, and the nature of these links, details of interactions between nodes (genes/proteins) are currently required. Importantly, our method can be applied directly to data on the attractors of a system when these data are incomplete. This can result either from only a subset of the full set of attractors being known, or from information being available for only a subset of nodes. In these cases, the method will still break up the attractors and identify relevant subdynamics. However, the more attractors there are, the more detailed and informative the eventual subsystems will be.
For many systems, data may be incomplete and/or noisy, leading to models that are also incomplete and/or unreliable. Even in the case of the extensively studied Drosophila segment polarity network, components may be missing and not fully understood. The most important factor affecting the quality of the results from our method (i.e. whether we find useful subsystems) is the reliability of the data/attractors. Even if there are some inaccuracies in a model or network topology, we can still make some reliable assertions, by focussing on the reliable data/attractors. In the case of the Drosophila segment polarity network model studied in this paper, 3 of the 10 attractors correspond to observed phenotypes. As we have already shown, applying the method directly to these 3 attractors gives us 5 'primary' subsystems (${S}_{A}^{\ast},{S}_{B1}^{\ast},{S}_{B2}^{\ast},{S}_{C1}^{\ast},{S}_{C2}^{\ast}$, see Fig. 7B). Then, once extra attractors from any model are taken into account, these 'primary' subsystems are split up so that each new subsystem only comes from a single 'primary' subsystem (e.g. in the existing model, the subsystem ${S}_{A}^{\ast}$ splits into S_{ A }, S_{D 2}, S_{E 2}, S_{H 2}). Therefore, for any future model of this system (with extra components or interactions), the subsystems obtained from our method should also relate back to these 5 primary subsystems ${S}_{A}^{\ast},{S}_{B1}^{\ast},{S}_{B2}^{\ast},{S}_{C1}^{\ast},{S}_{C2}^{\ast}$. Therefore, modules extracted from the currently available data/models are still informative, and those that only rely on observed data can be viewed as reliable (even if they are larger/less detailed than the 'true' subsystems). New data and improved models will just increase the precision of the results/subsystems. Another way in which we can assess subsystems (from a current model) is by looking at how robust they are. As well as looking at how robust each subsystem is to perturbations in node states, we can assess how robust each subsystem is to changes in the model itself (such as changes in Boolean functions or the addition of new components). We envisage that this method can be extended in a number of ways. Firstly, the method can be applied to signal transduction pathways by having 'n' attractors representing the state of the pathway(s) under 'n' different environmental/cellular conditions. These attractors can then be analysed with the new method to identify subsystems within pathway(s). Additionally, there is no need to restrict ourselves to studying cellular systems and we envisage that the method will be applicable in other fields.
The main method we have introduced could be applied directly to large scale datasets, to extract novel functional information. For an individual experiment, data on the state of a set of nodes (genes/proteins) are taken from a particular tissue, developmental stage and environmental condition (e.g. normal, high/low temperature, light/dark). Moreover, expression data are typically taken at discrete time points and are often converted to binary (e.g. 1 = 'expressed', 0 = 'not expressed'). Therefore, the results of each individual experiment would correspond to an individual discretestate discretetime attractor. Although we have not applied our method to a large datasets in this paper, we believe that following methodology could be used as a starting point

1.
Carry out multiple experiments in different environmental conditions/tissues/developmental stages, to give a range of discretestate attractors,

2.
Apply the method described in this paper, directly to this set of attractors (without using any model).
Then, identified subsystems would then be subsets of proteins, whose collective dynamics are conserved across data sets. In principle, this method could be applied to either time series data or fixed point data. However, since it is more difficult to obtain accurate data on time courses rather than steady states, we believe that the approach would be most likely to yield valuable information for fixed point/steady state data. We note here, that data from mutants are not necessary for such an approach. In fact, mutant data may not be as suitable, since each mutant would correspond to a different system. This type of analysis would be complementary to clustering, which primarily looks for groups of genes whose expression levels change in unison.
Reliability of data will be a major issue if trying to apply this method directly to experimental data. Expression data are very noisy and the less reliable the data, the less certain we will be about any identified subsystems. However, if can still be the case that novel groupings of genes/proteins could be found. Moreover, in the future, we envisage that these methods can be adapted so that we estimate the 'probability' of certain partial state sequences occurring in attractors, and then use this data to identify likely subsystems.
One limitation of this method is that it may have difficulties with systems with lots of cyclic attractors, each one containing a similar (but not identical) subdynamic. This is because we will get lots of similar intersection sequences (Definition 4). Since this definition provides the hierarchical backbone of the method, the method could struggle to identify informative partition sequences (Definition 5) or subsystems (Definition 6).
The method works best when (1) every attractor is a fixed point or (2) different cyclic subdynamics interact hierarchically or independently. This implies that the dynamics are partitioned in a strict hierarchical manner and there is no ambiguity when selecting subsystems. In case (1), an additional advantage is that all the fixed point attractors still occur when any asynchronous/stochastic updating scheme is used in a Boolean network model (once the model reaches a fixed point attractor, no update can cause the system to leave that attractor). Therefore, since subsystems are found directly from the attractors, the subsystems would also be the same. The situation of having all fixed point attractors is often relevant when considering cell type specification in developmental systems, where a cell must settle on one of a number of fixed states.
Methods
Boolean network models
For the purposes of this paper, a Boolean network model is a discrete time, deterministic, synchronous process acting on a directed network of v nodes V = {n_{1}, ..., n_{ v }}. At each discrete time step t ≥ 0, each node n_{ i }∈ V has a Boolean state s_{ i }(t) ∈ {0, 1} and these collectively form a network state x = x(t) = (s_{1}(t), ..., s_{ v }(t)). The model progresses, from one time step to the next, by synchronously updating these Boolean states according to a set of Boolean functions f = (f_{1}, ..., f_{ v }), as follows
 x(t + 1) = f(x(t)) = (f_{1}(x(t)), ..., f_{ v }(x(t))).
As time progresses, x(t) eventually gets trapped in an attractor A = {z_{0}, ..., z_{p1}}. i.e. there is a time point t' for which

For all t ≥ t', x(t) = z_{ k }(where k = t  t' (mod p)).
Here, each z_{ i }∈ A is called an attractor state
For a given model, there are typically multiple attractors and these correspond to the stable dynamics of the system. For the purposes of this paper $\mathbb{A}$ = {A_{1}, ..., A_{ r }} is the set of all attractors. A method to identify every attractor in a given Boolean network model is given in [13]. An example of a Boolean network model, along with a sample of its attractors, can be seen in Fig. 2.
It is possible to have models where nodes are updated in an asynchronous fashion. The method of identifying subsystems (below) can still be applied to such models as long as you focus on a specific set of attractors associated with such an updating scheme. The case where node updates are chosen stochastically is not explicitly considered here.
Algorithms
Below, we describe the main algorithms used in our new method. As can be seen in Fig. 1, there are multiple procedures/steps in our method, and so we describe each of these steps separately; namely

Identify every intersection sequence (satisfying Definition 4),

Identify every partition sequence (satisfying Definition 5),

Identify every subsystem (satisfying Definition 6).
Finally, if a suitable mathematical model is available, we show how regulation sets can be found for each subsystem (satisfying Definition 7). Formal proofs for all of the procedures can be found in Additional file 1 (Main Method) and Additional file 2 (Regulation of Subsystems).
However, before describing the main procedure, we need to be able to identify the partial state sequences that occur in each attractor, given a node set N. Therefore, we first describe some procedures relating to partial state sequences.
Algorithms: Partial state sequences
Given a node set N and an attractor A = {z_{0}, z_{1}, ..., z_{p1}}, the following procedure identifies a partial state sequence $P=\{{x}_{0}^{N},\mathrm{...},{x}_{q1}^{N}\}$ that occurs in A (i.e. the 3 properties of Definition 3 are satisfied)
Procedure 1.
Initially let k = 0, b_{0} = 0 and ${x}_{0}^{N}$ = {s_{ i }∈ z_{0} : n_{ i }∈ N}. The enter the following loop
Step 1: If k = p  1, let q* = b_{p1}+ 1 and go to step 6
Step 2: Let j = k and increment k by 1 (let k = k + 1)
Step 3: If ${x}_{{b}_{j}}^{N}$ = {s_{ i }∈ z_{ k }: n_{ i }∈ N}, then let b_{ k }= b_{ j }and go to step 1 (otherwise go to step 4)
Step 4: Let b_{ k }= b_{ j }+ 1
Step 5: Let ${x}_{{b}_{k}}^{N}$ = {s_{ i }∈ z_{ k }: n_{ i }∈ N} and go to step 1
Step 6: If ${x}_{{b}_{p1}}^{N}={x}_{{b}_{0}}^{N}$and q* > 1, reduce q* by 1 (let q* = q*  1)
Step 7: Let q be the smallest integer for which both
(a) qq* (this can be q = q*)
(b) ${x}_{f}^{N}={x}_{g}^{N}$, whenever f ≤ b_{p1}, g ≤ b_{p1}and f (mod q) = g (mod q)
Step 8: For k = 0, ..., p  1, let b_{ k }= b_{ k }(mod q)
end of procedure
A formal proof for this procedure can be be found in Additional file 1 (see Theorem S1.4 in section S1.1.1). However, here, we give a brief justification.
At the end of this procedure $P=\{{x}_{0}^{N},\mathrm{...},{x}_{q1}^{N}\}$ occurs in A and the 3 properties of Definition 3 are satisfied. Steps 2–5 ensure that properties 1 and 2 are satisfied and the partial states in P cycle within the attractor A in the correct order (as the attractor progresses over time). Steps 3, 6, 7 and 8 ensures property 3, so that P is the smallest possible set of partial states that cycles within A. In particular
(a) Steps 3 and 6 ensures no two adjacent partial states in P are identical
(b) Step 7 ensures that if a sequence of states cycles many times within an attractor, only one copy is kept.
This leaves a partial state sequence that just describes the 'order' in which the node states change in A (for nodes in N).
This procedure can be easily modified to look at a partial state sequence ${P}_{x}=\{{x}_{0}^{N},{x}_{1}^{N},\mathrm{...},{x}_{q1}^{N}\}$ (for a node set N) occurring in another partial state sequence ${P}_{y}=\{{y}_{0}^{M},\mathrm{...},{y}_{r1}^{M}\}$ (where M ⊇ N).
Given a node set N and a set of attractors $\u2102$, we need to be able to find a set of partial state sequences P_{1}, ..., P_{ k }that are all distinguishable from one another and optimally partition $\u2102$ into smaller sets ${\u2102}_{1},\mathrm{...},{\u2102}_{k}$.
One such way is to apply the following procedure (Procedure 2). This will identify partial state sequences P_{1}, ..., P_{ k }and sets of attractors ${\u2102}_{1},\mathrm{...},{\u2102}_{k}$ that satisfy properties A F below
A: For i = 1, ..., k, P_{ i }involves the node set N (i.e. ${P}_{i}=\{{x}_{{i}_{0}}^{N},\mathrm{...},{x}_{{i}_{q1}}^{N}\}$)
B: For i = 1, ..., k, P_{ i }occurs in every attractor A ∈ $\u2102$_{ i }
C: For i = 1, ..., k, P_{ i }does not occur in any attractor A ∉ $\u2102$_{ i }
D: For any i, j (1 ≤ i <j ≤ k), ${\u2102}_{i}\cap {\u2102}_{j}=\varnothing $
E: ${\u2102}_{1}\cup \mathrm{...}\cup {\u2102}_{k}=\u2102$
F: Given the node set N, there are no other partial state sequences P' ∉ {P_{1}, ..., P_{ k }} that occur in any attractor A ∈ $\u2102$ (unless P' and some P_{ i }contain the same partial states in the same order, within a rotation)
Procedure 2.
Begin with the node set N and set of attractors $\u2102$ and then carry out the following steps
Step 1: For every attractor A_{ j }∈ $\u2102$, apply Procedure 1 to N and A_{ j }, to get a partial state sequence Q_{ j }that occurs in A_{ j }.
Step 2: Put the Q_{ j }'s into groups i = 1, ..., k, whereby two partial state sequences ${{Q}^{\prime}}_{x}$, ${{Q}^{\prime}}_{y}$ go in the same group ⇔ ${{Q}^{\prime}}_{x}$ is equivalent to ${{Q}^{\prime}}_{y}$ (i.e. ${{Q}^{\prime}}_{x}$ contains the same partial states in the same order, within a rotation). Here, k is the minimum number of groups required to hold every Q_{ j }.
Step 3: For each group, i, let
i) P_{ i }= any Q_{ j }in the group i
ii) $\u2102$_{ i }= {A_{ j }: Q_{ j }is part of the group i}
end of procedure
A formal proof for this procedure can be be found in Additional file 1 (see Theorem S1.10 in section S1.1.3). However, here, we give a brief justification of why properties A – F are satisfied at the end. If a partial state sequence $P=\{{x}_{0}^{N},{x}_{1}^{N},\mathrm{...},{x}_{q1}^{N}\}$ occurs in an attractor A then so will q  1 other equivalent partial state sequences that contain the same partial states in the same order (within a rotation). e.g. ${P}^{\prime}=\{{x}_{1}^{N},\mathrm{...},{x}_{q1}^{N},{x}_{0}^{N}\}$ Moreover, if P occurs in an attractor A then no other partial state sequences for the same node set N can (other than the equivalent ones). Then, because

Each partial state sequence Q_{ j }(from Step 1) occurs in the attractor A_{ j }and involves a node set N,

Partial state sequences found in Step 1 are only grouped together (in Step 2) if they are equivalent,

Attractors are only grouped together in Step 3, if equivalent partial state sequences occur in them,
properties A, B, C and F must be satisfied. Properties D and E are satisfied because each attractor is put into exactly one set in Step 3.
Algorithms: Intersection sequences (Stage 1)
Identifying every intersection sequence is equivalent to finding every partial state sequence that satisfies the 3 properties of Definition 4 (for some set of attractors $\u2102$, say).
The method for identifying every intersection sequence can be visualised by considering the tree in Fig. 8. Searching through a tree analogous to this one (for a network with nodes V = {n_{1}, ..., n_{ v }}), means that every node set N can be visited at some point. Then, using Procedure 2, we can identify the partial state sequences that occur in each attractor (for each node set N). Then, after the tree has been fully examined, we can pick out the partial state sequences (P) and sets of attractors $\u2102$ that satisfy the 3 properties of Definition 4. In reality, many branches of the tree can be ignored, leading to improvements in efficiency (discussed below).
We first give the procedure for identifying every intersection sequence, then give an example and finally discuss ways to make the process more efficient.
Procedure 3
First, consider the tree in Fig. 8 and note that for every node set N, there exists a path from left to right (starting at '') that corresponds to it. Therefore, searching through a tree analogous to the one in Fig. 8 (for a network with nodes V = {n_{1}, ..., n_{ v }}), every node set N can be visited at some point.
The procedure searches through the tree (as in Fig. 8) and carries out the following steps for each node set N. At the end of the procedure the set S contains every intersection sequence
Step 0 (Initialisation): Let S = ∅ and let N = ∅ ()
Step 1: Move onto the next node set N in the tree (as in Fig. 8).
Step 2: For the node set N, apply Procedure 2 to identify partial state sequences P_{1}, ..., P_{ k }and sets of attractors ${\u2102}_{1},\mathrm{...},{\u2102}_{k}$ satisfying
A: For i = 1, ..., k, P_{ i }involves the node set N (i.e. ${P}_{i}=\{{x}_{{i}_{0}}^{N},\mathrm{...},{x}_{{i}_{q1}}^{N}\}$)
B: For i = 1, ..., k, P_{ i }occurs in every attractor A ∈ $\u2102$_{ i }
C: For i = 1, ..., k, P_{ i }does not occur in any attractor A ∉ $\u2102$_{ i }
D: For any i, j (1 ≤ i <j ≤ k), ${\u2102}_{i}\cap {\u2102}_{j}=\varnothing $
E: ${\u2102}_{1}\cup \mathrm{...}\cup {\u2102}_{k}=\mathbb{A}$ (the set of all attractors)
F: Given the node set N, there are no other partial state sequences P' ∉ {P_{1}, ..., P_{ k }} that occur in any attractor A ∈ $\mathbb{A}$
Step 3: For i = 1, ..., k, add the pair {P_{ i }, $\u2102$_{ i }} to the set S
Step 4: For i = 1, ..., k, check S to see if there is any pair {$Q=\{{y}_{0}^{M},\mathrm{...},{y}_{r1}^{M}\},\mathbb{D}$} for which either of the following are true
(a) M ⊂ N and $\mathbb{D}={\u2102}_{i}$
(b) M ⊃ N and $\mathbb{D}={\u2102}_{i}$
If (a) is true, remove {Q, $\mathbb{D}$} from S. If (b) is true, remove {P_{ i }, $\u2102$_{ i }} from S
Step 5: If the tree has been completely searched, end procedure. Otherwise, return to step 1.
end of procedure
A formal proof for this procedure can be be found in Additional file 1 (see Theorem S1.14 in section S1.2.1). However, here, we give a brief justification.
At the end of the procedure, S gives a complete set of intersection sequences (satisfying the 3 properties of Definition 4). Step 2 ensures every partial state sequence that satisfies properties 1 and 2 are identified for each node set N. Step 4 then ensures that only those satisfying property 3 remain in S.
As an example of how the procedure finds intersection sequences, consider the model in Fig. 2 and the intersection sequence P_{3} (from Fig. 3). When the node set N = {n_{2}, n_{3}, n_{4}} is analysed in Step 2 of the algorithm (using Procedure 2) we find two distinct partial state sequences that satisfy properties 1 and 2 of Definition 4, namely

P_{3} that intersects at $\u2102$_{3} = {A_{1}, A_{2}, A_{3}}

P_{ x }= {s_{2} = 0, s_{3} = 0, s_{4} = 0} that intersects at $\u2102$_{ x }= {A_{4}}
The remainder of the algorithm involves deciphering whether of these partial state sequences satisfy the final property (Property 3 of Definition 4). This can be done by looking at the way the attractors are split for different node sets. For the node set N = {n_{2}, n_{3}, n_{4}} (above) we split the attractors as follows  {A_{1}, A_{2}, A_{3}} + {A_{4}}. For larger nodes sets M ⊃ N
M_{1} = {n_{1}, n_{2}, n_{3}, n_{4}} splits the attractors as follows  {A_{1}, A_{2}} + {A_{3}} + {A_{4}}
M_{5} = {n_{2}, n_{3}, n_{4}, n_{5}} splits the attractors as follows  {A_{1}} + {A_{2}} + {A_{3}} + {A_{4}}
M_{6} = {n_{2}, n_{3}, n_{4}, n_{6}} splits the attractors as follows  {A_{1}} + {A_{2}} + {A_{3}} + {A_{4}}
etc
Therefore, when P_{3} and $\u2102$_{3} = {A_{1}, A_{2}, A_{3}} are considered in Step 4, it can never be removed from the set S. This is because when a larger node set M ⊃ N is considered, we will never have a set of attractors $\u2102={\u2102}_{3}$ (or $\mathbb{D}\supseteq {\u2102}_{3}$). Therefore, P_{3} satisfies the property 3 and is an intersection sequence. However, P_{ x }is not an intersection sequence since the corresponding set of attractors ({A_{4}}) is still associated with a larger node set M ⊃ N. Therefore it will be removed from the set S in Step 4 (either when N or M is analysed in the procedure).
We now explain how the procedure can be made more efficient
Improving efficiency(1)
Consider node sets N ⊂ M and two partial state sequences $P=\{{x}_{0}^{N},\mathrm{...},{x}_{q1}^{N}\}$ and ${P}^{\prime}=\{{y}_{0}^{M},\mathrm{...},{y}_{r1}^{M}\}$ that both occur in some attractor A. Then, P' must contain P and hence only occur in an attractor whenever P does (proved in Lemma S1.12 of Additional file 1).
Therefore, if N ⊂ M ⊂ V, and P only occurs in a single attractor A*, neither P nor P' can be intersection sequences. This is since the attractor A* is itself a partial state sequence (for a larger node set V) that occurs in exactly the same set of attractors (i.e. just A*), and so property 3 of Definition 4 fails.
Therefore, if Step 2 of Procedure 3 identifies a partial state sequence ${P}_{i}=\{{x}_{0}^{N},\mathrm{...},{x}_{q1}^{N}\}$ that only occurs in a single attractor $\u2102$_{ i }= {A*}, we know that reanalysing this attractor {A*} for any node set M ⊃ N is pointless.
We show one way in which this knowledge can improve efficiency in Procedure 3. Suppose, attractors ${{A}^{\prime}}_{1},\mathrm{...},{{A}^{\prime}}_{f}$ were returned as single attractors in Step 2, when analysing an earlier node set P ⊆ N (including P = N). Then, when a path is extended to the right in the tree (Fig. 8) from a node set N to a node set M ⊃ N, Procedure 2 in Step 2 need only be applied to the smaller set of attractors
Moreover, if every attractor is returned as a single attractor in Step 2, when analysing an earlier node set P ⊆ N (including P = N), there is no need to extend the path to look at node sets M ⊃ N. In Fig. 8, this is equivalent to ignoring all longer paths that include extra nodes to the right. For example, if N = {n_{1}, n_{3}}, there would be no need to look at longer paths (from left to right) that give node sets M = {n_{1}, n_{3}, n_{4}}, M = {n_{1}, n_{3}, n_{5}} or M = {n_{1}, n_{3}, n_{4}, n_{5}}.
However, because all of the attractors are intersection sequence, the full node set V = {n_{1}, ..., n_{ v }} should still be fully analysed in Steps 2 – 4 (possibly at the very end of the procedure).
Improving efficiency(2)
As can be seen in Fig. 8, some nodes appear less than others, with the least frequent nodes visited earlier in the tree. Therefore, it is likely to be advantageous to reindex nodes in the tree during the search. At any stage during the search, nodes along paths to the right (from a node set N) can be reindexed without impairing our ability to search the tree. For example, once N = {n_{1}, n_{3}} has been reached, reindexing nodes {n_{4}, n_{5}} to {n_{5}, n_{4}} still allows us to reach the same node sets M = {n_{1}, n_{3}, n_{4}}, M = {n_{1}, n_{3}, n_{5}} and M = {n_{1}, n_{3}, n_{4}, n_{5}}, as before. However, they would be visited in a different order (M = {n_{1}, n_{3}, n_{5}} then M = {n_{1}, n_{3}, n_{4}} then M = {n_{1}, n_{3}, n_{4}, n_{5}}).
Once a node set N has been analysed, reindexing so that the next node n_{ j }to be visited maximises c (below) will speed up the search
 For the sets of attractors ${\u2102}_{1},\mathrm{...},{\u2102}_{k}$ identified in Step 2 (for the new node set M = N ∪ {n_{ j }}), $\u2102$_{ i }= {A_{ i }} is a single attractor for c (≤ k) different values of i
Although this involves carrying out Step 2 multiple times (to compare different n_{ j }'s), selecting an n_{ j }that gives lots of single attractor $\u2102$_{ i }'s will mean less analysis later on (as discussed above). The quicker we can reach a stage where every set $\u2102$_{ i }in Step 2 is a single attractor, the more of the tree can be ignored during the rest of the search.
Algorithms: Partition sequences (Stage 1)
In order to identify every partition sequences, we start with the full set of intersection sequences and all the corresponding sets of attractors (obtained from Procedure 3). i.e. every pair {P, $\u2102$} such that P intersects at $\u2102$. Then, these intersection sequences are used to identify all the partial state sequences that satisfy any of properties A, B or C of Definition 5. Using the complete list of intersection sequences, properties A, B and C are considered independently and partition sequences stored after each test.
Part A: Core components
From Procedure 3, we get a set S that contains the complete set of intersection sequences, along with the set of attractors each one intersects at (if {P', $\u2102$} S, then P' intersects at $\u2102$).
Using this set S as an input, the following procedure identifies every partial state sequence that is core to some set of attractors (i.e every partial state sequence P satisfying Definition 5A)
Procedure 4. A
Initially, let the set T = ∅ (empty set)
Then, for every intersection sequence ${P}^{\text{'}}=\{{y}_{0}^{\mathrm{N\text{'}}},{y}_{1}^{\mathrm{N\text{'}}},\mathrm{...},{y}_{r1}^{\mathrm{N\text{'}}}\}\in \{{P}^{\text{'}},\u2102\}\in S$, carry out the following steps
Step 1: From the complete set of intersection sequences (S), identify every Q_{ i }(for the node set M_{ i }) for which
(a) Q_{ i }intersects at $\mathbb{D}$_{ i }, where ${\mathbb{D}}_{i}\cap \u2102\ne \varnothing $
(b) There is no intersection sequence Q* (for a larger node set M* ⊃ M_{ i }) that intersects at ${\mathbb{D}}^{\ast}\supseteq {\mathbb{D}}_{i}\cap \u2102$
Step 2: Let k be the number of partial state sequences from Step 1
Step 3: Let N = M_{1} ∩ ... ∩ M_{ k }(N ⊆ N' since P' is itself identified in Step 1)
Step 4: If N ≠ ∅ in Step 3, find a partial state sequence $P=\{{x}_{0}^{N},{x}_{1}^{N},\mathrm{...},{x}_{q1}^{N}\}$ that occurs in P' (see Procedure 1).
Step 5: If N ≠ ∅ in Step 3, add the pair {P, $\u2102$} to the set T
end of procedure
At the end of the procedure, T contains every partial state sequence P that is core to some set of attractors $\u2102$ (i.e every partial state sequence P satisfying Definition 5A). Essentially, we take each intersection sequence P' in turn and then rerun though the set of all intersection sequences to find those (Q_{ i }, for a node set M_{ i }) that satisfy the following
(a) Q_{ i }cooccurs with P' in at least one attractor (Q_{ i }can be P').
(b) Q_{ i }isn't contained in a larger intersection sequence (Q*, say) that cooccurs with P' in the same attractors
Then, the node set N = M_{1} ∩ ... ∩ M_{ k }is the core set of nodes underlying P' and the attractors it occurs in ($\u2102$). The partial state sequence P that satisfies property A is then just partial state sequence (for the node set N) that occurs in P'
A formal proof for this procedure can be be found in Additional file 1 (see Theorem S1.17 in section S1.2.2). However, here, we give a brief justification.
The 3 properties of Definition 5A are satisfied, for all the identified partial state sequences P. Property 2 is satisfied because of the following. If an intersection sequence Q (for a node set M) intersects at a set of attractors $\mathbb{D}$ (where $\mathbb{D}\cap \u2102\ne \varnothing $), then there is an intersection sequence Q' (for a node set M') identified in step 1 for which

M ⊆ M' ⊆ N (because of Steps 1 and 3)

Q' occurs in every attractor A ∈ $\mathbb{D}\cap \u2102$ (because of Step 1)
and so

M ∪ N ⊇ M'

Q' occurs in every attractor A ∈ $\mathbb{D}\cap \u2102$
Property 1 is satisfied because of Step 4. Property 3 is satisfied because of the way N is chosen in Step 3. We note that every partial state sequence satisfying Definition 5A is identified because (a) Property 1 implies that such a sequence must occur in an intersection sequence and (b) every intersection sequence is analysed independently.
Part B: Exclusive (Procedure 4B)
This is simply done by searching through all intersection sequences (in S) and identifying those that satisfy Definition 5B.
Part C: Independently Oscillating (Procedure 4C)
This is simply done by searching through every pair of intersection sequences (in S) to see which pairs satisfy Definition 5C. Where such a pair is found, both of them are partition sequences.
Algorithms: Subsystems (Stage 2)
In order to identify every subsystem, we start with the full set of partition sequences from Stage 1. i.e. every partial state sequence satisfying properties either A, B and C of Definition 5, The following procedure identifies every subsystem (satisfying Definition 6)
Procedure 5.
Initially, let the set U = ∅ (empty set)
Then, for every partition sequence $P=\{{y}_{0}^{M},{y}_{1}^{M},\mathrm{...},{y}_{r1}^{M}\}$ (identified in A, B or C of the previous section), carry out the following steps
Step 1: From the complete set of partition sequences, identify every partition sequence P_{ i }(for a node set M_{ i }) for which
(a) M_{ i }⊂ M
(b) P_{ i }and P both cooccur in some attractor A
(This also implies that P_{ i }occurs in P)
Step 2: Let k be the number of partition sequences from Step 1
Step 3: Let N = M\(M_{1} ∪ ... ∪ M_{ k })
Step 4: If N ≠ ∅, identify $S=\{{x}_{0}^{N},{x}_{1}^{N},\mathrm{...},{x}_{q1}^{N}\}$ that occurs in P (see Procedure 1)
Step 5: If N ≠ ∅, add S (identified in Step 4) to the set U
end of procedure
For each partition sequence P (for a node set M, say), we look to see which other partition sequences ${{P}^{\prime}}_{1},\mathrm{...},{{P}^{\prime}}_{k}$ (for node sets ${{M}^{\prime}}_{1},\mathrm{...},{{M}^{\prime}}_{k}$) occurs within it. Then the node set N, which is unique to P, can be calculated as N = M\(M_{1} ∪ ... ∪ M_{ k }). If N ≠ ∅, the partial state sequence S (for the node set N) that occurs in P, is stored as a subsystem.
A formal proof for this procedure can be be found in Additional file 1 (see Theorem S1.19 in section S1.3). However, here, we give a brief justification.
At the end of the procedure, U contains every subsystem P satisfying the 3 properties of Definition 6. Essentially, property 1 is satisfied because of Step 4. Property 2 is satisfied because of Step 1 and the choice of N in step 3. Property 3 is satisfied because N is the largest set for which M_{ i }∩ N = ∅ for all i = 1, ..., k. We note that every partial state sequence satisfying Definition 6 is identified because (a) Property 1 implies that such a sequence must occur in an partition sequence and (b) every partition sequence is analysed independently.
Algorithms: Regulation of subsystems
In order to show how regulation sets can be identified for each subsystem, we first focus on how partial states are regulated in Boolean network models (although all the methods are applicable to other discretestate discretetime logical models).
It may be that a partial state x^{N}controls some Boolean functions {f_{ i }: n_{ i }∈ M} and ensure the occurrence of y^{M}in the following time step. i.e x^{N}is a predecessor of y^{M}(or x^{N}triggers the occurrence of y^{M}). For example, in Fig. 2, by following the Boolean functions it is possible to say that

The occurrence of ${x}_{3}^{N}$ = {s_{2} = 1, s_{3} = 1} can be triggered by the occurrence of ${x}_{2}^{N}$ = {s_{1} = 1, s_{2} = 1}.

The occurrence of ${x}_{2}^{N}$ = {s_{1} = 1, s_{2} = 1} can be triggered by the occurrence of ${x}_{1}^{N}$ = {s_{1} = 1}.
or ${x}_{1}^{N}$ = {s_{1} = 1} triggers the occurrence of ${x}_{3}^{N}$ = {s_{2} = 1, s_{3} = 1} after two time steps.
This notion of predecessors can be extended to look for predecessors within attractors and subsystems. Going backwards around an attractor it is possible to identify which partial states are responsible for triggering the occurrence of other partial states in the same attractor state at a later point in time. Returning to the example in Fig. 2, and looking at attractor A_{1}, it is possible to go backwards around the attractor via two different routes to say
Route 1:
The occurrence of {s_{2} = 1, s_{3} = 1} in z_{0} can be triggered by the occurrence of {s_{1} = 1, s_{2} = 1} in z_{3}
The occurrence of {s_{1} = 1, s_{2} = 1} in z_{1} can be triggered by the occurrence of {s_{1} = 1} in z_{2}
The occurrence of {s_{1} = 1} in z_{2} can be triggered by the occurrence of {s_{1} = 1} in z_{1}
The occurrence of {s_{1} = 1} in z_{3} can be triggered by the occurrence of {s_{1} = 1} in z_{0}
Route 2:
The occurrence of {s_{2} = 1, s_{3} = 1} in z_{0} can be triggered by the occurrence of {s_{2} = 1, s_{3} = 1} in z_{3}
The occurrence of {s_{2} = 1, s_{3} = 1} in z_{1} can be triggered by the occurrence of {s_{2} = 1, s_{3} = 1} in z_{2}
The occurrence of {s_{2} = 1, s_{3} = 1} in z_{2} can be triggered by the occurrence of {s_{2} = 1, s_{3} = 1} in z_{1}
The occurrence of {s_{2} = 1, s_{3} = 1} in z_{3} can be triggered by the occurrence of {s_{2} = 1, s_{3} = 1} in z_{0}
and so we have

1.
{s_{1} = 1} triggers the occurrence of {s_{2} = 1, s_{3} = 1} in the attractor state z_{0}

2.
{s_{2} = 1, s_{3} = 1} triggers the occurrence of {s_{2} = 1, s_{3} = 1} in the attractor state z_{0}
(Note: Similar conclusions could also be made for z_{1}, z_{1} and z_{3})
So, more generally, given a partial state y^{M}that occurs in an attractor state z ∈ A, we want to be able to find a suitable collection of predecessors ${x}_{1}^{{N}_{1}},\mathrm{...},{x}_{k}^{{N}_{k}}$ (that trigger the occurrence of y^{M}in z).
Procedure 6. Given a partial state y^{M}that occurs in an attractor state z ∈ A, this procedure identifies partial states ${x}_{1}^{{N}_{1}},\mathrm{...},{x}_{k}^{{N}_{k}}$ that trigger the occurrence of y^{M}in z.
Starting from the partial state y^{M}in z, this procedure involves going backwards around the attractor A (as in the above example) and identifying new predecessors at each step (using a method adapted from [13]). The procedure ends when no more new predecessors ${x}_{j}^{{N}_{j}}$ (that trigger the occurrence of y^{M}in z) can be found.
Since much of this method involves adapting a previously published algorithm (from [13]), we do not give full details here. More details of this procedure can be found in section S2.1 of Additional file 2 (in particular, Procedures S2.8 and S2.11)
end of procedure
We now want to apply these procedures to subsystems. In particular, we want to find collections of subsystems ${\mathbb{S}}_{x}=\{{S}_{{x}_{1}},\mathrm{...},{S}_{{x}_{y}}\}$ whose cooccurrence in an attractor triggers chain of interactions that results in the occurrence of S_{ y }.
Returning to the above example (from Fig. 2), we carry out Procedure 6 on S_{1} = {s_{2} = 1, s_{3} = 1} in all of the attractor states in A_{1}, A_{2} and A_{3} (which contain S_{1}), which gives the following results.
A_{1} and A_{2}: S_{1} = {s_{2} = 1, s_{3} = 1} can trigger the occurrence of S_{1} = {s_{2} = 1, s_{3} = 1} in the attractor states z_{0}, z_{1}, z_{2} and z_{3}
A_{1} and A_{2}: S_{2} = {s_{1} = 1} can trigger the occurrence of S_{1} = {s_{2} = 1, s_{3} = 1} in the attractor states z_{0}, z_{1}, z_{2} and z_{3}
A_{3}: S_{1} = {s_{2} = 1, s_{3} = 1} can trigger the occurrence of S_{1} = {s_{2} = 1, s_{3} = 1} in the attractor states z_{0} and z_{1}
Therefore, by looking at all of the attractor states (individually) in each attractor, we can say that
(a) S_{2} can trigger the occurrence of S_{1} in attractors A_{1} and A_{2}.
(b) S_{1} can trigger the occurrence of itself (S_{1}) in attractors A_{1}, A_{2} and A_{3}.
This is then sufficient to give a regulation set for S_{1}; namely $\mathbb{S}$_{1} = {S_{1}} and $\mathbb{S}$_{2} = {S_{2}}. In this very simple example, the same subsystems/partial states are involved in every attractor state within an attractor. However, it may be case that different subsystems are involved in triggering different partial states (from a subsystem S_{ y }) in different attractor states, and so it is necessary to consider each attractor state individually.
Procedure 7 (below) demonstrates a method of identifying a regulation set for a subsystem ${S}_{y}=\{{y}_{0}^{M},\mathrm{...},{y}_{p1}^{M}\}$. This is done by looking at every attractor state in every attractor (containing S_{ y }) to see how what triggers the occurrence of each relevant partial state ${y}_{i}^{M}$.
However, since partial states within subsystems could be subject to different time lags in different attractors (see Definition 3), we first note that we need to consider the precise dynamics (the instances) of subsystems in attractors. i.e.
Definition 9. Consider a collection of subsystems $\mathbb{S}$ = {S_{1}, ..., S_{ f }} where every ${S}_{i}=\{{x}_{{i}_{0}}^{{N}_{i}},\mathrm{...},{x}_{{i}_{{q}_{i}}}^{{N}_{i}}\}\in \mathbb{S}$ involves a node set N_{ i }and occurs in the attractor A = {z_{0}, ..., z_{p1}}
The instance of $\mathbb{S}$ in A is the partial state sequence ${z}_{0}^{M},\mathrm{...},{z}_{p1}^{M}$, where

1.
M = N_{1} ∪ ... ∪ N_{ f }

2.
For k = 0, ..., p  1, ${z}_{k}^{M}$ = {s_{ x }∈ z_{ k }: n_{ x }∈ M}.
We can then use this terminology to describe whether a collection of subsystems $\mathbb{S}$ = {S_{1}, ..., S_{ f }} triggers an individual subsystem S_{ y }in an attractor A. i.e If the cooccurrence of subsystems S_{1}, ..., S_{ f }ensures the occurrence of S_{ y }in A.
Definition 10. Suppose we have

1.
An attractor A = {z_{0}, ..., z_{p1}}

2.
A collection of subsystems ${\mathbb{S}}_{x}=\{{S}_{{x}_{1}},\mathrm{...},{S}_{{x}_{y}}\}$ where
(a) ${S}_{{x}_{1}},\mathrm{...},{S}_{{x}_{f}}$ all occur in A
(b) ${x}_{0}^{N},\mathrm{...},{x}_{p1}^{N}$ is the instance of $\mathbb{S}$_{ x }in A

3.
An individual subsystem S_{ y }where
(a) S_{ y }occurs in A
(b) ${y}_{0}^{M},\mathrm{...},{y}_{p1}^{M}$ is the instance of S_{ y }in A
Then $\mathbb{S}$_{ x }triggers the occurrence of S_{ y }in A if the following holds for every i ∈ {0, ..., p  1}

$${x}_{i}^{N}$$
triggers the occurrence of ${y}_{i}^{M}$ in the attractor state z_{ i }∈ A.
We now give the procedure for identifying a regulation set of a subsystem S_{ y }(Procedure 7). At the start of the procedure, we take a subsystem S_{ y }(involving a node set M_{ y }) and a set of attractors $\u2102$_{ y }, where
(a) S_{ y }occurs in every attractor A ∈ $\u2102$_{ y }
(b) S_{ y }does not occur in any attractor A ∉ $\u2102$_{ y }
Moreover, we assume we know every subsystem $S=\{{x}_{0}^{N},\mathrm{...},{x}_{q1}^{N}\}$, along with the node set involved (N) and a list of attractors it occurs in. Each subsystem and node set N is a byproduct of the method of identifying subsystems (described previously in this Methods section). A list of attractors can be found by applying Procedure 2 to the node set N (and the full set of attractors). This information is used in Step 2 of the procedure (below).
Procedure 7.
Initially, let the set R = ∅ (empty set)
For every A_{ i }= {z_{0}, ..., z_{p1}} ∈ $\u2102$_{ y }carry out the following steps.
Step 1: Let the set R_{ i }= ∅. Let the sets U_{0}, ..., U_{p1}= ∅
Step 2: Identify every subsystem T_{1}, ...., T_{ h }that occurs in A_{ i }. Moreover, let M_{1}, ..., M_{ h }be the node sets involved in T_{1}, ...., T_{ h }(resp)
Step 3: Identify the instance of S_{ y }in A_{ i }. i.e. ${y}_{0}^{{M}_{y}},\mathrm{...},{y}_{p1}^{{M}_{y}}$.
(The procedure for this is obvious from Definition 9, given the node set M_{ y }and attractor A_{ i })
Step 4: For j = 0, ..., p  1, carry out Procedure 6 to identify predecessors of ${y}_{j}^{{M}_{y}}$ in z_{ j }∈ A_{ i }(i.e. the partial states that trigger the occurrence of ${y}_{j}^{{M}_{y}}$ in z_{ j }). The resulting predecessors are added to the set U_{ j }
Step 5: For every possible combination of partial states ${x}_{0}^{{N}_{0}},\mathrm{...},{x}_{p1}^{{N}_{p1}}$ satisfying

$${x}_{j}^{{N}_{j}}$$
∈ U_{ j }(for j = 0, ..., p  1)
carry out the following
(a) Let N = N_{0} ∪ ... ∪ N_{p1}
(b) Let $\mathbb{S}$ = {T_{ a }: M_{ a }∩ N ≠ ∅}
(c) Add $\mathbb{S}$ to the set R_{ i }
Step 6: Remove all subsystem collections $\mathbb{S}$ from R_{ i }that contain other subsystem collections $\mathbb{S}$' ∈ R_{ i }. (i.e. $\mathbb{S}\supset {\mathbb{S}}^{\prime}$)
Step 7: Add the subsystem collections in R_{ i }to the set R
end of procedure
At the end of this procedure, every subsystem collection $\mathbb{S}$ ∈ R_{ i }triggers the occurrence of S_{ y }in the attractor A_{ i }(i.e. Definition 10 is satisfied). For every attractor state z_{ j }∈ A_{ i }, we identify the partial state ${y}_{j}^{{M}_{y}}$ ∈ S_{ y }that occurs in z_{ j }(Step 3). Then, starting from the partial state ${y}_{j}^{{M}_{y}}$ in the attractor state z_{ j }, we go backwards around the attractor, identifying suitable predecessors at each time step (this backwards process can go multiple times around the attractor). This will then give us a list of partial states that can trigger the occurrence of ${y}_{j}^{{M}_{y}}$ in z_{ j }(added to U_{ j }Step 4). Having done this for every attractor state z_{ j }∈ A_{ i }, Step 5 then pulls the results together to find collections of subsystems $\mathbb{S}$ = {S_{1}, ..., S_{ f }} that triggers the occurrence of S_{ y }in A (i.e. those collections that satisfy Definition 10). In particular, finding collections for which the following is true.
For every attractor state z_{ j }∈ A_{ i }, there exists a partial state ${x}_{j}^{{N}_{j}}$ that satisfies
(a) ${x}_{j}^{{N}_{j}}$ triggers the occurrence of ${y}_{j}^{{M}_{y}}$ in z_{ j }
(b) ${x}_{j}^{{N}_{j}}$ only involves nodes and states from the subsystems {S_{1}, ..., S_{ f }}
Step 6 just ensures that only the most informative collections are kept and there is no redundancy. Finally, since these collections of subsystems are added to the set R in Step 7, for every attractor (which contains S_{ y }), R = {${\mathbb{S}}_{1},\mathrm{...},{\mathbb{S}}_{g}$} is a regulation set and satisfies the properties of Definition 7.
A formal proof for this procedure can be be found in Additional file 2 (see Theorem S2.20 in section S2.2.2). However, this procedure shows just one way of a regulation set, for a subsystem S_{ y }. There may be more than one possible regulation set for a subsystem and some may be more descriptive than others. In Section S2.2.1 and S2.2.2 of Additional file 2, we give some some extra constraints (and corresponding procedures) that can be applied when looking for regulation sets.
References
 1.
Hartwell L, Hopfield J, Leibler S, Murray A: From molecular to modular cell biology. Nature (supp) 1999, 402: C47C52. 10.1038/35011540
 2.
Yu H, Gerstein M: Genomic analysis of the hierarchical structure of regulatory networks. PNAS 2006, 103: 14724–14731. 10.1073/pnas.0508637103
 3.
Newman M, Girvan M: Finding and evaluating community structure in networks. Phys Rev E 2004, 69: 026113. 10.1103/PhysRevE.69.026113
 4.
ShenOrr S, Milo R, Mangan S, Alon U: Network motifs in the transcriptional regulation network of Escherichia coli . Nat Genet 2002, 31: 64–68. 10.1038/ng881
 5.
ResendisAntonio O, FreyreGonzález J, MenchacaMéndez R, GutiérrezRios R, MartinezAntonio A, ÁvilaSánchez C, ColladoVides J: Modular analysis of the transcriptional regulatory network of E. coli . Trends Genet 2005, 21: 16–20. 10.1016/j.tig.2004.11.010
 6.
Ingram P, Stumpf M, Stark J: Network motifs: structure does not determine function. BMC Genomics 2006, 7: 108. 10.1186/147121647108
 7.
BarJoseph Z, Gerber G, Lee T, Rinaldi N, Yoo J, Robert F, Gordon D, Fraenkel E, Jaakkola T, Young R, Gifford D: Computational discovery of gene modules and regulatory networks. Nat Biotech 2003, 21(11):1337–1342. 10.1038/nbt890
 8.
Kato M, Hata N, Banerjee N, Futcher B, Zhang M: Identifying combinatorial regulation of transcription factors and binding motifs. Genome Biol 2004, 5(8):1–13. 10.1186/gb200458r56
 9.
Chen K, Calzone L, CsikaszNagy A, Cross F, Novak B, Tyson J: Integrative Analysis of Cell Cycle Control in Budding Yeast. Mol Biol Cell 2004, 15: 3841–3862. 10.1091/mbc.E03110794
 10.
Albert R, Othmer H: The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila . J Theor Biol 2003, 223: 1–18. 10.1016/S00225193(03)000353
 11.
Chaves M, Albert R, Sontag E: Robustness and fragility of Boolean models for genetic regulatory networks. J Theor Biol 2005, 235: 431–449. 10.1016/j.jtbi.2005.01.023
 12.
Ma W, Lai L, Ouyang Q, Tang C: Robustness and modular design of the Drosophila segment polarity network. Mol Syst Biol 2006, 2():70. 10.1038/msb4100111
 13.
Irons D: Improving the efficiency of attractor cycle identification in Boolean networks. Physica D 2006, 217: 7–21. 10.1016/j.physd.2006.03.006
 14.
Grossniklaus U, Pearson R, Gehring W: The Drosophila sloppy paired locus encodes two proteins involved in segmentation that show homology with mammalian transcription factors. Genes & Dev 1992, 6: 1030–1051. 10.1101/gad.6.6.1030
Acknowledgements
We would like to thank C. Cannings, S. Bullock and K. Walters for critical reading. This work was supported by funding from BBSRC (grant ref: 02/B1/B/08370) and EPSRC (grant ref: EP/D003105/1).
Author information
Additional information
Authors' contributions
DI designed the method of identifying subsystems, coconceived the study and drafted the manuscript. NM coconceived the study and drafted the manuscript. All authors read and approved the final manuscript.
Electronic supplementary material
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Received
Accepted
Published
DOI
Keywords
 Attractor State
 Boolean Function
 Partial State
 Intersection Sequence
 Boolean Network Model