Identifying dynamical modules from genetic regulatory systems: applications to the segment polarity network
 David J Irons^{1}Email author and
 Nicholas AM Monk^{1, 2}
DOI: 10.1186/147121058413
© Irons and Monk; licensee BioMed Central Ltd. 2007
Received: 24 May 2007
Accepted: 25 October 2007
Published: 25 October 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.
Results: Identifying subsystems and regulatory interactions between subsystems
 1.
Conserved across some set of attractors,
 2.
Distinguishable from the stable dynamics for the rest of the system.
(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.
 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
 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).

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}_{0}^{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})${x}_{1}^{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
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.
 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.
Intersection sequences in the example Boolean network model
Intersection sequences  Intersects at ... 

P _{3}  $\u2102$_{3} = {A_{1}, A_{2}, A_{3}} 
P _{4}  $\u2102$_{4} = {A_{3}, A_{4}} 
P _{5}  $\u2102$_{5} = {A_{1}, A_{2}} 
P _{6}  $\u2102$_{6} = {A_{1}, A_{2}} 
P _{7}  $\u2102$_{7} = {A_{1}} 
P _{8}  $\u2102$_{8} = {A_{2}} 
P _{9}  $\u2102$_{9} = {A_{3}} 
P _{10}  $\u2102$_{10} = {A_{4}} 
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
 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
 1.
$\u2102\cap \mathbb{D}$ ≥ 2
 2.
N ∪ M = V (the set of all nodes)
Partition sequences in the example Boolean network model
Partition sequences  Core to ...  Exclusive to ...  Independently Oscillating with ... 

P _{1}  {A_{1}, A_{2}, A_{3}}  ×  × 
P _{2}  {A_{1}, A_{2}}  ×  × 
P _{3}  ×  {A_{1}, A_{2}, A_{3}}  P _{5} 
P _{4}  {A_{3}, A_{4}}  {A_{3}, A_{4}}  × 
P _{5}  ×  ×  P_{3} and P_{6} 
P _{6}  ×  ×  P _{5} 
P _{7}  {A_{1}}  {A_{1}}  × 
P _{8}  {A_{2}}  {A_{2}}  × 
P _{9}  {A_{3}}  {A_{3}}  × 
P _{10}  {A_{4}}  {A_{4}}  × 
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.
 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.
Unique components of partition sequences in the example Boolean network model
Partition sequence  Unique Component 

P _{1}  S _{1} 
P _{2}  S _{2} 
P _{3}  S _{3} 
P _{4}  S _{4} 
P _{5}  S _{5} 
P _{6}  N/a 
P _{7}  N/a 
P _{8}  N/a 
P _{9}  N/a 
P _{10}  S _{6} 
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 }.
 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 }.
Regulation of subsystems in the example Boolean network model
Subsystem  Regulation set 

S _{1}  $\mathbb{S}$ = {S_{1}} 
$\mathbb{S}$ = {S_{2}}  
S _{2}  $\mathbb{S}$ = {S_{2}} 
S _{3}  $\mathbb{S}$ = {S_{1}, S_{3}} 
$\mathbb{S}$ = {S_{2}, S_{3}}  
S _{4}  $\mathbb{S}$ = {S_{4}} 
S _{5}  $\mathbb{S}$ = {S_{2}, S_{5}} 
S _{6}  $\mathbb{S}$ = {S_{4}, S_{6}} 
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.

S_{ x }occurs in an attractor A ⇒ S_{ y }occurs in an attractor A
 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.
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
Other subsystems for the segment polarity network model
Subsystem  States 

S _{D 1}  wg_{1} = 1, WG_{1} = 1 
S _{D 2}  wg_{1} = 0, WG_{1} = 0 
S _{E 1}  wg_{2} = 1, WG_{2} = 1 
S _{E 2}  Wg_{2} = 0, WG_{2} = 0 
S _{F 1}  PTC_{1} = 1 
S _{F 2}  PTC_{1} = 0 
S _{G 1}  PTC_{2} = 1 
S _{G 2}  PTC_{2} = 0 
S _{H 1}  CIA_{1} = 1, ptc_{1} = 1 
S _{H 2}  CIA_{1} = 0, ptc_{1} = 0 
S _{I 1}  CIA_{2} = 1, ptc_{2} = 1 
S _{I 2}  CIA_{2} = 0, ptc_{2} = 0 
S _{J 1}  CIR_{1} = 1, CIR_{2} = 1 
S _{J 2}  CIR_{1} = 0, CIR_{2} = 0 
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.
Regulation of subsystems in the segment polarity network model
Subsystem  Regulation set 

S _{ A }  $\mathbb{S}$ = {S_{ A }, S_{B 1}, S_{C 1}} 
$\mathbb{S}$ = {S_{ A }, S_{B 2}, S_{C 1}}  
$\mathbb{S}$ = {S_{ A }, S_{B 1}, S_{C 2}}  
$\mathbb{S}$ = {S_{ A }, S_{B 2}, S_{C 2}}  
S _{B 1}  $\mathbb{S}$ = {S_{ A }, S_{B 1}} 
S _{B 2}  $\mathbb{S}$ = {S_{ A }, S_{B 2}, S_{E 2}} 
$\mathbb{S}$ = {S_{ A }, S_{B 2}, S_{C 2}, S_{F 1}, S_{G 1}}  
S _{C 1}  $\mathbb{S}$ = {S_{ A }, S_{C 1}} 
S _{C 2}  $\mathbb{S}$ = {S_{ A }, S_{C 2}, S_{D 2}} 
$\mathbb{S}$ = {S_{ A }, S_{B 2}, S_{C 2}, S_{F 1}, S_{G 1}}  
S _{D 1}  $\mathbb{S}$ = {S_{D 1}, S_{H 1}, S_{J 2}} 
S _{D 2}  $\mathbb{S}$ = {S_{D 2}} 
$\mathbb{S}$ = {S_{B 1}}  
$\mathbb{S}$ = {S_{E 1}}  
$\mathbb{S}$ = {S_{ A }, S_{C 2}, S_{F 1}}  
S _{E 1}  $\mathbb{S}$ = {S_{E 1}, S_{I 1}, S_{J 2}} 
S _{E 2}  $\mathbb{S}$ = {S_{E 2}} 
$\mathbb{S}$ = {S_{C 1}}  
$\mathbb{S}$ = {S_{D 1}}  
$\mathbb{S}$ = {S_{ A }, S_{B 2}, S_{G 1}}  
S _{F 1}  $\mathbb{S}$ = {S_{B 2}, S_{C 1}} 
$\mathbb{S}$ = {S_{ A }, S_{C 2}, S_{F 1}}  
S _{F 2}  $\mathbb{S}$ = {S_{F 2}, S_{H 2}} 
$\mathbb{S}$ = {S_{B 1}, S_{C 1}}  
S _{G 1}  $\mathbb{S}$ = {S_{B 1}, S_{C 2}} 
$\mathbb{S}$ = {S_{ A }, S_{B 2}, S_{G 1}}  
S _{G 2}  $\mathbb{S}$ = {S_{G 2}, S_{I 2}} 
$\mathbb{S}$ = {S_{B 1}, S_{C 1}}  
S _{H 1}  $\mathbb{S}$ = {S_{B 2}, S_{C 1}} 
S _{H 2}  $\mathbb{S}$ = {S_{B 1}} 
$\mathbb{S}$ = {S_{E 1}}  
$\mathbb{S}$ = {S_{ A }, S_{C 2}, S_{F 1}}  
S _{I 1}  $\mathbb{S}$ = {S_{B 1}, S_{C 2}} 
S _{I 2}  $\mathbb{S}$ = {S_{C 1}} 
$\mathbb{S}$ = {S_{D 1}}  
$\mathbb{S}$ = {S_{ A }, S_{B 2}, S_{G 1}}  
S _{J 1}  $\mathbb{S}$ = {S_{ A }, S_{B 2}, S_{C 2}, S_{F 1}, S_{G 1}} 
S _{J 2}  $\mathbb{S}$ = {S_{B 1}} 
$\mathbb{S}$ = {S_{C 1}}  
$\mathbb{S}$ = {S_{F 2}, S_{G 2}} 

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

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

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

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

= S_{C 2}∪ S_{G 1}${S}_{C2}^{\ast}$
Robustness of attractors in the segment polarity network model
Attractor  Basin of attraction  Robustness 

A _{1}  0.01 %  0.67 
A _{2}  98.63 %  1.0 
A _{3}  1.27 %  1.0 
A _{4}  0.00035 %  0.60 
A _{5}  0.00005 %  0.58 
A _{6}  0.04 %  0.69 
A _{7}  0.0003 %  0.60 
A _{8}  0.00015 %  0.58 
A _{9}  0.04 %  0.69 
A _{10}  0.01 %  0.67 
Robustness of subsystems in the segment polarity network model
Subsystem  Robustness  

Global  Internal  External  
S _{ A }  1.0  1.0  1.0 
S _{B 1}  1.0  1.0  1.0 
S _{B 2}  0.78  0.56  0.85 
S _{C 1}  1.0  1.0  1.0 
S _{C 2}  0.78  0.56  0.85 
S _{D 1}  0.63  0.5  0.63 
S _{D 2}  0.99  0.88  1.0 
S _{E 1}  0.63  0.5  0.63 
S _{E 2}  0.99  0.88  1.0 
S _{F 1}  0.76  0.71  0.76 
S _{F 2}  0.97  0.33  0.99 
S _{G 1}  0.76  0.71  0.76 
S _{G 2}  0.97  0.33  0.99 
S _{H 1}  0.73  1.0  0.72 
S _{H 2}  1.0  1.0  1.0 
S _{I 1}  0.73  1.0  0.72 
S _{I 2}  1.0  1.0  1.0 
S _{J 1}  1.0  1.0  1.0 
S _{J 2}  1.0  1.0  1.0 
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.
 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))).

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

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

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).
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.

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.
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.

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

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

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).

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}
 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}}
 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.
 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
 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

triggers the occurrence of ${y}_{i}^{M}$ in the attractor state z_{ i }∈ A.${x}_{i}^{N}$
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 }

∈ U_{ j }(for j = 0, ..., p  1)${x}_{j}^{{N}_{j}}$
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.
Declarations
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).
Authors’ Affiliations
References
 Hartwell L, Hopfield J, Leibler S, Murray A: From molecular to modular cell biology. Nature (supp) 1999, 402: C47C52. 10.1038/35011540View Article
 Yu H, Gerstein M: Genomic analysis of the hierarchical structure of regulatory networks. PNAS 2006, 103: 14724–14731. 10.1073/pnas.0508637103PubMed CentralView ArticlePubMed
 Newman M, Girvan M: Finding and evaluating community structure in networks. Phys Rev E 2004, 69: 026113. 10.1103/PhysRevE.69.026113View Article
 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/ng881View ArticlePubMed
 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.010View ArticlePubMed
 Ingram P, Stumpf M, Stark J: Network motifs: structure does not determine function. BMC Genomics 2006, 7: 108. 10.1186/147121647108PubMed CentralView ArticlePubMed
 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/nbt890View Article
 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/gb200458r56View Article
 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.E03110794PubMed CentralView ArticlePubMed
 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)000353View ArticlePubMed
 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.023View ArticlePubMed
 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/msb4100111PubMed CentralView ArticlePubMed
 Irons D: Improving the efficiency of attractor cycle identification in Boolean networks. Physica D 2006, 217: 7–21. 10.1016/j.physd.2006.03.006View Article
 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.1030View Article
Copyright
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.