A compositional framework for Boolean networks

BioSystems xxx (xxxx) xxxx

Contents lists available at ScienceDirect

BioSystems journal homepage: www.elsevier.com/locate/biosystems

A compositional framework for Boolean networks H. Alkhudhayra,b, J. Stegglesa,

⁎

a b

School of Computing, Newcastle University, Newcastle upon Tyne, UK Faculty of Computing and Information Technology, King Abdulaziz University, Rabigh, Saudi Arabia

ARTICLE INFO

ABSTRACT

Keywords: Boolean network Model composition Behaviour preservation Interference

Boolean networks are a widely used qualitative approach for modelling and analysing biological systems. However, their application is restricted by the well-known state space explosion problem which means that modelling large-scale, realistic biological systems is challenging. In this paper we set out to facilitate the construction and analysis of large scale biological models by developing a formal framework for the composition of Boolean networks. The compositional approach we present is based on merging entities between Boolean networks using a binary Boolean operator and we formalise the preservation of behaviour under composition using a notion of compatibility. We investigate characterising compatibility in terms of the composed models by developing a trace alignment property. In particular, we use a formalisation of the interference that can occur in a composed model to define an extended trace alignment property that we show completely characterises compatibility.

1. Introduction Qualitative modelling techniques have become increasingly important in biology as the basis for developing formal techniques and tools for modelling, analysing and synthesizing biological systems (De Jong, 2002; Bartocci and Lió, 2016). Boolean networks (Kauffman, 1969, 1993) are a widely used qualitative approach based on representing the state of each biological entity (e.g. genes, proteins and other chemical signals) abstractly as a Boolean value, where 1 represents the entity is active and 0 that it is inactive. Entities then interact to regulate their states based on predefined next-state functions and the resulting dynamic behaviour gives rise to attractor cycles that can then be associated with biological phenomena. Entities can update their state either synchronously, where the state of all entities is updated simultaneously, or asynchronously, where entities update their state independently. Boolean networks have been shown to allow a range of interesting biological analysis to be performed and have been widely considered in the literature (for example, see Akutsu et al., 1999; Steggles et al., 2007; Saadatpour and Albert, 2013; Rosenblueth et al., 2014; Bartocci and Lió, 2016). It seems clear that they have an important role to play in advancing our understanding and engineering capability of complex biological systems. However, one important issue that limits the scalable application of Boolean networks is the well-known state space explosion problem which means that modelling and analysing largescale, realistic biological systems is challenging.

⁎

In this paper we set out to facilitate the construction and analysis of large scale biological models by developing a formal framework for the composition of (synchronous) Boolean networks. The compositional approach we present is based on composing Boolean networks by merging entities using a binary Boolean operator (there are 16 such binary Boolean operators to choose from, see for example Kurgalin and Borunov, 2018). We consider what it means to preserve the underlying behaviour of Boolean networks in a composed model and introduce the notion of compatibility to formalise this concept. It turns out that the compatibility property is problematic to check as it references the behaviour of the composed model and so we develop the alignment property which is able to predict whether the composition of Boolean networks is compatible based only on their individual behaviour. We show formally that the alignment property is sufficient to ensure compatibility when the underlying Boolean operator being used is idempotent and illustrate its application by presenting results about the compatibility of composing duplicate copies of a Boolean network. While the alignment property is interesting it turns out that is not a necessary property for compatibility and so does not completely characterise behaviour preservation. We therefore investigate extending the alignment property by considering the behavioural interference that can occur when two Boolean networks are composed. We formalise this interference using a labelled state graph approach and then use this to define an extended alignment property called weak alignment. We show formally that weak alignment completely characterises compatibility

Corresponding author. E-mail addresses: [email protected], [email protected] (H. Alkhudhayr), [email protected] (J. Steggles).

https://doi.org/10.1016/j.biosystems.2019.04.004 Received 20 March 2018; Received in revised form 4 February 2019; Accepted 8 April 2019 0303-2647/ © 2019 Elsevier B.V. All rights reserved.

Please cite this article as: H. Alkhudhayr and J. Steggles, BioSystems, https://doi.org/10.1016/j.biosystems.2019.04.004

BioSystems xxx (xxxx) xxxx

H. Alkhudhayr and J. Steggles

when the underlying Boolean operator being used is idempotent and therefore provides an important scalable means of checking behaviour preservation in composed models. The compositional framework developed here is supported by a prototype tool that automates the composition process and associated analysis. We note that an initial version of this work which focused on the sufficient property of alignment appears in Alkhudhayr and Steggles (2017). This paper is organized as follows. In Section 2 we provide a brief introduction to Boolean networks. In Section 3 we develop a compositional framework for Boolean networks and introduce the notion of compatibility to formalise the preservation of behaviour under composition. In Section 4 we investigate characterising compatibility and introduce the alignment property, a sufficient property for ensuring compatibility. In Section 5 we extend this alignment result to provide a complete characterisation of compatibility by modelling the interference that can occur within a composed model. Finally, in Section 6 we present some concluding remarks and discuss future work.

represents the state of entity gi . Note as a notational convenience we often use s1 … sn to represent a global state (s1, …, sn). When the current state of a Boolean network is clear from the context we allow gi to denote both the name of an entity and its corresponding current state. The state space of a Boolean network , denoted S , is = |G|. therefore the set of all possible global states S The state of a Boolean network can be updated either synchronously (Kauffman, 1993; Wuensche, 2004), where the state of all entities is updated simultaneously in a single update step, or asynchronously (Harvey and Bossomaier, 1997), where entities update their state independently. In the following we focus on the synchronous update semantics which has received considerable attention in the literature (see for example Kauffman, 1969, 1993; Akutsu et al., 1999; Wuensche, 2004; Schaub et al., in press; Banks and Steggles, 2012). Given two

S2 represent a (synchronous) update step states S1, S2 S , let S1 such that S2 is the state that results from simultaneously updating the state of each entity gi using its associated update function F(gi) and the appropriate neighbourhood of states from S1. As an example, consider the global state 011 for Ex1 (see Fig. 1), where entity g1 = 0, g2 = 1, Ex1

101 is an update step in and g3 = 1. Then 011 Ex1. from some initial state is The sequence of global states through S called a trace. Note that in the case of the synchronous update semantics such traces are deterministic and infinite. Given that the global state space is finite, this implies that a trace must eventually enter a cycle, known formally as an attractor cycle (Kauffman, 1993; Thieffry and Thomas, 1995). Attractor cycles are important biologically where they are seen as representing different biological states or functions (e.g. different cellular types such as proliferation, apoptosis and differentiation Huang and Ingber, 2000). We define a finite canonical representation for synchronous traces which specifies the infinite behaviour of a trace up to the first repeated state as follows. Let S0 S be a global state for . A trace is a list of global states σ(S0) = 〈S0, S1, …, Sn+1〉 such that:

2. Boolean networks Boolean networks (Kauffman, 1969, 1993) are a widely used qualitative modelling approach for biological control systems (see for example Akutsu et al., 1999; Steggles et al., 2007; Saadatpour and Albert, 2013; Rosenblueth et al., 2014; Bartocci and Lió, 2016). In this section we introduce the basic definitions for Boolean networks needed in the sequel and provide illustrative examples. A Boolean network consists of a set of regulatory entities G = {g1, …, gn } which can be in one of two possible states, either 1 representing the entity is active (e.g. a gene is expressed or a protein is present) or 0 representing the entity is inactive (e.g. a gene is not expressed or a protein is absent). The state of each entity is regulated by a subset of entities in the Boolean network and we refer to this subset as the neighbourhood of an entity (an entity may or may not be in its own neighbourhood). An entity updates its state by applying a logical nextstate function to the current states of the entities in its neighbourhood. We can define a Boolean network more formally as follows. Definition 1. A Boolean Network

is a tuple

Si + 1, for 0 ≤ i ≤ n; (ii) S0, …, Sn are unique states; and (iii) (i) Si Sn+1 = Si for some i ∈ {0, …, n}. ) = { (S ) | S S } therefore comThe set of all traces Tr( pletely characterizes the behaviour of a Boolean network under the (synchronous) update semantics. For example, in Ex1 (see Fig. 1) the

= (G, N , F ) where:

trace 011

(i) G = {g1, …, gn} is a non-empty, finite set of entities; (ii) N = (N (g1), …, N (gn)) is a tuple of neighbourhoods, such that N(gi) ⊆ G is the neighbourhood of gi; and (iii) F = (F (g1), …, F (gn)) is a tuple of next-state functions, such that the function F (gi ): |N (gi)| defines the next state of gi. As an example, consider the Boolean network Ex1 = (G Ex1, NEx1, FEx1) defined in Fig. 1. It consists of three entities GEx1 = {g1, g2, g3} with neighbourhoods NEx1(g1) = {g2}, NEx1(g2) = {g1}, and NEx1(g3) = {g1, g2}. The next-state functions FEx1 are defined equational in Fig. 1.(B), where we use [gi] to represent the next state of an entity gi. A global state of a Boolean network with n entities is represented by a tuple of Boolean states (s1, …, sn), where si

Fig. 1. Example of a Boolean network graph.

Ex1

Ex1

101

Ex1

010

Ex1

101

Ex1

… is denoted by

(011) = 011, 101, 010, 101 It can be seen that Ex1 has three attractors: two point attractors 〈000, 000〉 and 〈110, 110〉; and a cyclic attractor 〈101, 010, 101〉. The behaviour of a Boolean network can be concisely represented by a state graph in which the nodes are the global states and the edges are precisely the synchronous update steps allowed. We let

)=( , ) denote the state graph for a Boolean network under the synchronous trace semantics. As an example, consider the synchronous state graph SG( Ex1) for Ex1 presented in Fig. 1(C).

SG(

consisting of: (A) wiring diagram; (B) equational definition of next-state functions for

2

Ex1;

(C) synchronous state

BioSystems xxx (xxxx) xxxx

H. Alkhudhayr and J. Steggles

(i) If g1 ∉ N1(g1) and g2 ∉ N2(g2), then

(l1, …, lp, k1, …, kq) = F1 (g1)(l1, …, lp)

F2 (g 2)(k1, …, kq );

(ii) If g1 ∈ N1(g1) and g2 ∉ N2(g2), then

(g c , l1, …, lp, k1, …, kq) = F1 (g1)(g c , l1, …, lp)

F2 (g 2)(k1, …, kq);

(iii) If g1 ∉ N1(g1) and g2 ∈ N2(g2), then

(g c , l1, …, lp, k1, …, kq) = F1 (g1)(l1, …, lp)

F2 (g 2)(g c , k1, …, kq);

(iv) If g1 ∈ N1(g1) and g2 ∈ N2(g2), then Fig. 2. Pictorial representation of composing Boolean network by merging entities g1 entity gc.

1 1

and and g 2

2

to form a new 2 into a new

(g c , l1, …, lp, k1, …, kq) = F1 (g1)(g c , l1, …, lp)

In the sequel, we let gc denote the new entity created by merging g1 1 2 ( and g2 and assume has global states 1, 2, g , g ) (g c g11…gn1 g12…gm2 ) . As an example, consider composing Ex1 (Fig. 1) and Ex2 (Fig. 3) on entities g1 and g4 using conjunction ∧. The resulting Boolean network ( Ex1, Ex2, g1, g 4 ) is depicted in Fig. 4. We would like to be able to infer properties and behaviour of a composed model from the underlying Boolean networks used in the composition. Being able to do this would allow us to construct compositionally large Boolean networks with known properties and so help to address the limitations imposed by the state space explosion problem. The following definitions formalize the idea that the original behaviour of the underlying Boolean networks can be preserved in their composition. We begin by defining projection operators which are able to extract states and traces from a composed system.

3. Compositional framework In this section we begin to develop our formal framework for composing Boolean networks based on the idea of merging entities using a binary Boolean connective (for example, conjunction). We consider what it means for the behaviour of an individual Boolean network in a composed model to be preserved and formalise this by introducing a notion of compatibility. In the sequel, let 1 = (G1, N1, F1) and 2 = (G2, N2 , F2 ) be two Boolean networks such that G1 = {g1, g11, …, gn1} and G2 = {g 2, g12 , …, gm2 } are disjoint sets of entities, for some n, m . We let ⊙ represent an arbitrary binary Boolean operator chosen from the 16 possible binary Boolean operators available (see for example Kurgalin and Borunov, 2018). We formally define the composition of two Boolean networks 1 and 2 using ⊙ to merge the behaviour of entities (see Fig. 2) as follows.

( Definition 2. (Composition) Let 1, Boolean network constructed by composing g1 and g2 using ⊙ defined as follows:

1 2 ( Definition 3. (Projections) Let = 1, 2, g , g ) be the new 1 Boolean network constructed by composing 1 and 2 on entities g and g2 using ⊙. Let S = (g c g11…gn1 g12…gm2 ) be a global state in the composed system. Then we define the left 1: 1 and right 2: 2 projection operators by

g1, g 2) denote the and 2 on entities

2, 1

1 (S ) 1

1. Entities: The finite set of entities G = (G1/{g }) ∪ (G2/ {g2}) ∪ {gc}, where gc denotes the new entity created by merging g1 and g2 . 2. Neighbourhood: For any entity hi ∈ G, the neighbourhood N(hi) is defined as follows:

N (h i ) =

= (g c g11…gn1),

2 (S )

= (g c g12…gm2 )

We can extend the projection operators to traces by 1(

)=

1 (S1),

1 (S2),

… ,

2(

)=

2 (S1),

= S1, S2, …

2 (S2 ),

Tr( )

…

and let 1 (Tr( )) and 2 (Tr( )) represent the sets of projected traces derived by projecting each trace in Tr( ) . Note that projected traces may not be well-defined traces in their corresponding Boolean network, i.e. j (Tr( )) Tr( j) may hold, for j ∈ {1, 2}. We are interested in situations where composing two Boolean networks preserves their underlying trace behaviour and define a notion of compatibility to formalise this.

N1 (hi )[g1/ g c ], if hi G1 N2 (hi )[g2 / g c ], if hi G2 N1 (g1)[g1/ g c ] N2 (g 2)[g 2/ g c ], if hi = g c

where S[f/e] represents set S with all occurrences of element f replaced by e. 3. Functions: For any hi ∈ G, the next-state function F(hi) is defined:

F (h i ) =

F2 (g 2 )(g c , k1, …, kq).

1 2 = ( Definition 4. (Compatibility) Let 1, 2, g , g ) be the Boolean network resulting from composing 1 and 2 on entities 1 g1 and g2 using ⊙. Then we say that 1 and 2 are compatible on g 2 and g under ⊙ iff Tr( 1) 1 (Tr( )) and Tr( 2) 2 (Tr( )) .

F1 (hi ), if hi G1 F2 (hi ), if hi G2 , if hi = g c

Let N1(g1) = {l1, …, lp} if g1 ∉ N1(g1), otherwise let N1(g1) = {g1, l1, …, lp}. Similarly, let N2(g2) = {k1, …, kq} if g2 ∉ N2(g2), otherwise let c using four N2(g2) = {g2, k1, …, kq}. Then we define : |N (g )| cases as follows:

To illustrate the definition of compatibility consider composing and using conjunction ∧ to produce Ex1 Ex2 = ( Ex1, Ex2, g1, g4 ) (see Fig. 4). Then examples of projected traces in 2 (Tr( )) (assuming state order (gc g2 g3 g5)) will be Fig. 3. A second Boolean network example Ex2 containing the wiring diagram, next-state equations, and state graph.

3

BioSystems xxx (xxxx) xxxx

H. Alkhudhayr and J. Steggles

Fig. 4. Boolean network is (gc g2 g3 g5)). 2( 2(

(

Ex1,

Ex2,

g1, g4 ) resulting from the composition of

0100, 1011, 0100 ) = 00, 11, 00 1001, 0100, 1011, 0100 ) = 11, 00, 11

2( 2(

Ex1

and

Ex2

on entities g1 and g4 using conjunction ∧ (where the state order

and let g1 on g1 and g2 iff

0001, 0001 ) = 01, 01 1100, 1100 ) = 10, 10

Tr( It can be seen that and Ex1) 1 (Tr( )) Tr( Ex2) 2 (Tr( )) and so we have that Ex1 and Ex2 are compatible on g1 and g4 under ∧. Compatibility is important since it ensures that properties of the composed model can be inferred from the underlying models used in the composition which are normally much smaller and so simpler to analyse. In particular, if two Boolean networks are compatible then the behaviour they exhibit (such as attractors) must be present in the composed model. Since a composed model may extend the behaviour of its underlying models it does not allow the absence of behaviour to be verified. However, the work presented in Section 5 on interference graphs and in particular, Theorem 4 provides a starting point for addressing these types of negative properties. Given a Boolean operator ⊙ which is commutative and associative, we can show that the composition of Boolean networks under ⊙ is commutative and associative (see Alkhudhayr and Steggles, 2017 for more detail). This result is important as it means that the order in which multiple Boolean networks are composed under a commutative and associative Boolean operator does not affect the resulting model.

g1 (Tr(

1))

=

1

and g 2 g 2 (Tr(

2.

We say that

1

and

2

are aligned

2))

1 2 ( In the context of a composition = 1, 2, g , g ) we can define the merging of global states and traces as follows. Let 2 2 2 2 S1 = (g1 g11…gn1) 2 . Then we define 1 and S = (g g1 …gm ) 1 2 S S by merging the state of g1 with g2, that is S1 S 2 = (g1 g 2 g11…gn1 g12…gm2 ) . Let 1 = S11, S21, … Tr( 1) and 2 2 Tr( be two traces. Then we define 2 = S1 , S2 , … 2) 1 S12, S21 S22, … . Note that for any 1 Tr( 1 2 = S1 1) and Tr( Tr( ) . 2 2) we may have that 1 2 1 1 1 Lemma 1. Let 1 and 2 be Boolean networks with G1 = {g , g1 , …, gn } 2 , g 2 , …, g 2 } G = { g and 2 1 m . Let ⊙ be an idempotent binary Boolean operator 1 2 Tr( Tr( ( and let = 2) 1) and 2 1, 2, g , g ) . Let 1 Tr( ) ; and (ii) such that g1 ( 1) = g 2 ( 2 ) . Then we have: (i) 1 2 1( 1 2 ) = 1 and 2 ( 1 2) = 2 .

Tr( Tr( Proof. Let 1 = S1, S2, … 1) , 2 = T1, T2, … 2) , and suppose for any i we have Si = (s i s1i…sni) and Ti = (t i t1i…tmi ) . Then by the assumption g1 ( 1) = g 2 ( 2 ) we know si = ti and so since ⊙ is assumed to be idempotent we have

4. Compatibility and alignment

si

Compatibility is an important behavioural preservation property to consider when composing Boolean networks but unfortunately it is difficult to check since it refers to the full behaviour of the composed model. In this section we begin to address this issue by investigating how to infer compatibility without using the composed model. We formalise a trace alignment property which we show is sufficient for obtaining compatibility in a composed system. We use this result to show that duplicate Boolean networks are compatible under composition of corresponding entities. For any Boolean network with entities G = {g1, …, gn}, global state S = (s1…sn) S and any entity gi we define the state projection gi (S ) = si . Then gi ( ) denotes the projected trace of entity gi = S1, S2, … Tr( ) on trace defined by )) ={ gi ( )| Tr( )} . gi ( ) = gi (S1), gi (S2 ), … . We let gi (Tr( As an example, consider projecting the traces of Ex2 (Fig. 3) on g4 which gives

Then the proof for (i) and (ii) are straightforward based on (I) (see Alkhudhayr and Steggles, 2017 for an example of this proof based on using conjunction). □ We can now prove that alignment is a sufficient property for compatibility.

g4 (Tr(

Ex2))

1 2 ( Proof. Let = 1 and 2 are 1, 2, g , g ) and assume that aligned on g1 and g2. By the definition of compatibility (Definition 4) we have two properties to show:

Tr( (i) Tr( 1) we know by the 1) 1 (Tr( )) : For any 1 alignment assumption above that there exists 2 Tr( 2) such that Tr( ) . Then 2 g 1 ( 1) = g 2 ( 2 ) . It follows by Lemma 1.(i) that 1 by Lemma 1.(ii) we have 1 ( 1 2 ) = 1 and so 1 1 (Tr( )) as required. (ii) Tr( 2) 2 (Tr( )) : The proof follows along similar lines to (i) above. □

We can now define the property of alignment as follows. 1

and

2

(I)

1 Theorem 2. Let 1 and 1 2 be two Boolean networks with g 2 and g 2 . Let ⊙ be an idempotent binary Boolean operator. Then if 1 2 1 and 1 and 2 are aligned on g and g then 2 are compatible 1 2 on g and g under ⊙.

= { 0, 1, 0 , 0, 0 , 1, 1 , 1, 0, 1 }.

Definition 5. (Alignment) Let

t i = si = t i

be two Boolean networks 4

BioSystems xxx (xxxx) xxxx

H. Alkhudhayr and J. Steggles

5.1. Labelled state graphs

Fig. 5. Composing two duplicate copies of and g4 using conjunction.

Ex1

)=( , ) is defined to be the state graph for Recall that SG( )) denote the a Boolean network (see Section 2). Let Path(SG( ) . It is straightforward to set of all infinite paths in the state graph SG( see that under the synchronous update semantics traces and infinite paths are equivalent and in the sequel we use these terms interchangeably. Next we introduce the idea of a labelled state graph to allow the projection of a Boolean network's behaviour on a subset of its entities (similar to the approach used in the previous section for traces). Let be a Boolean network and let L be any non-empty set of labels. L is termed a state labelling function. Then any function : = S1, S2, … Path(SG( )) we define Given an infinite path ( ) = (S1), (S2 ), … and

on corresponding entities g1

The above result provides a means of ensuring compatibility holds without requiring the composed system to be considered. This is important since a composed model will be larger and so more affected by the state space explosion problem. Note that while alignment is a sufficient condition for compatibility it can be shown that it is not a necessary property for it. We consider this further in the next section where we extend the alignment property to address this. Two Boolean networks are said to be duplicates if they have the same structure and behaviour up to the renaming of entities (i.e. they are isomorphic). It turns out that duplicate Boolean networks are always aligned on corresponding entities (where corresponding is defined in the obvious way based on the underlying isomorphism) and therefore are compatible under an idempotent Boolean operator. As an illustration, consider the example presented in Fig. 5 based on composing two duplicate copies of Ex1 (Fig. 1) using conjunction. We can therefore use the alignment property to formally show that duplicate Boolean networks are always compatible when composed on corresponding entities using an idempotent Boolean operator.

(Path(SG(

))) = { ( ) |

Path(SG(

))}

We can now define a labelled state graph as follows. Definition 6. (Labelled State Graph) Let be a Boolean network and L be a state labelling function. Then we define the let : )) by labelled state graph (SG(

(SG(

)) = (

,

,

)

))) = (Path(SG( ))) . We define Path( (SG( 1 2 ( Let = 1 and 1, 2, g , g ) be the result of composing 1 2 2 on g and g using ⊙. Then using the projection functions and 2: 2 defined in Section 3 we can extract the relevant behaviour from corresponding to 1 and 2 by using the labelled state graphs 1 (SG( )) and 2 (SG( )) respectively. As an example, consider = ( Ex1, Ex2, g1, g4 ), the Boolean network resulting from the composition of Ex1 and Ex2 on entities g1 and g4 using conjunction (see Fig. 4). Then the labelled state graph 1 (SG( )) shown in Fig. 7 is the result of projecting from the relevant behaviour corresponding to Ex1. The equivalence between synchronous traces and infinite paths in the state graph means that we can define compatibility and alignment properties in terms of (labelled) state graphs. For example, it is straightforward to show that two Boolean networks 1 and 2 are compatible on g1 and g2 under ⊙ iff we have Path(SG( Path Path(SG( Path( 1(SG( ))) and 2)) 1)) ( 2 (SG( ))) , where

Theorem 3. Let 1 and 2 be two duplicate Boolean networks and let 2 g1 1 and g 2 be corresponding entities. Let ⊙ be an idempotent 1 2 binary Boolean operator. Then 1 and 2 are compatible on g and g under ⊙. Proof. Since 1 and 2 are duplicates it follows (assuming a corresponding state order) that Tr( by 1) = Tr( 2) . Thus Definition 5 we know that 1 and 2 are aligned on corresponding entities g1 and g2. Since ⊙ is assumed to be idempotent 1 we therefore have by Theorem 2 that 1 and 2 are compatible on g and g2 under ⊙, as required. □

5. Fully characterising compatibility

=

(

1,

2,

g1, g 2).

We also have that

The alignment property is a sufficient property for ensuring compatibility but is limited as it does not completely characterise the compatibility property. To illustrate this consider the counter example ( Ex1, Ex3, g1, g4 ) (see Fig. 6(B)), which is the Boolean network resulting from composing Ex3 (Fig. 6(A)) on g1 Ex1 (Fig. 1) and and g4 using conjunction. It can be seen that Ex3 are Ex1 and compatible on g1 and g4 under ∧. However, Ex3 do not Ex1 and Tr( align on g1 and g4 since there is a trace 00, 00 Ex1) whose projection g1 ( 00, 00 ) = 0, 0 is not in g4 (Tr( Ex3)) . Therefore, we have that the alignment property is a sufficient but not a necessary property for compatibility. In this section we extend the alignment property to provide a complete characterisation of compatibility. We do this by considering the interference that can occur to the behaviour of Boolean networks when they are composed. We develop a labelled state graph approach to formalise this interference and then use this to define a new alignment property referred to as weak alignment. We show that weak alignment completely characterises compatibility (i.e. it is a sufficient and necessary property for compatibility).

Path(

g1 (SG(

1)))

1

and

= Path(

2

align on g1 and g2 iff we have

g 2 (SG(

2))).

These definitions based on (labelled) state graphs provide a means of developing efficient algorithms for checking the compatibility and alignment properties, and we discuss this further in Section 6. 5.2. Modelling interference using state graphs When two Boolean networks are composed by merging entities it is possible for the behaviour of the underlying models to experience interference resulting in new behaviour. To illustrate this, consider ( Ex1, Ex2, g1, g 4 ) (see Fig. 4) and suppose the composed model is in global state 0111. Then the underlying Boolean networks will want Ex1

Ex2

01. In 101 and 01 to make the following transitions: 011 other words, entity g1 in Ex1 would like to transition from state 0 to state 1 but entity g4 in Ex2 would like to transition 0 to 0. The merged entity gc will therefore transition to 0 ∧ 1 =0 meaning that the underlying behaviour of Ex1 has been interfered with. This shows that when using conjunction interference will occur to the underlying behaviour of a merged entity whenever it wants to transition to 1 but its merged counterpart wants to transition to 0. We can formalise this

5

BioSystems xxx (xxxx) xxxx

H. Alkhudhayr and J. Steggles

Fig. 6. A counter example which shows that alignment is not a necessary property for compatibility. (A) The Boolean network ( Ex1, Ex3, g1, g4 ) resulting from the composition of Ex3 on entities g1 and g4 using ∧. Ex1 and

observation by defining an interference set Δ = {1} for conjunction. For each of the 16 possible different binary Boolean operators we can identify when interference can occur by considering their truth tables and identifying if the value of an argument is preserved after applying the operator. This results in a set of next state values which can be interfered with when using a Boolean operator (for example, for conjunction we have Δ = {1} and for disjunction Δ = {0}). There are four such possible interference sets {}, {0}, {1} and {0, 1}. Note that since not all of the Boolean operators are commutative the order (first or second) of the Boolean network being composed can affect the interference set and for this reason we normally provide an interference set for the first Boolean network Δ1 and for the second Δ2. The above idea is summarised in Table 1 where the interference sets for each possible binary Boolean operator is presented. In order to formalise the possible interference behaviour that can occur between composed Boolean networks we extend a Boolean network's state graph with additional edges representing interference. The idea is that whenever the entity to be merged transitions to a state s that can experience interference s ∈ Δ then we add another edge to the state graph to represent that the transition could instead go to ¬s. This situation is depicted in Fig. 8 based on the interference example presented above for Ex1 using conjunction. We formalise this idea by defining an interference state graph as

Ex3 ;

and (B) Boolean network

follows. be a Boolean network Definition 7. (Interference state graph) Let with entities g, g1, …, gn and let be an interference set. Then we ) for define the interference state graph SGg, ( on g with interference Δ by

SGg, (

)=(

,

g,

)

The extended edge relation

g,

is defined by

= {((s s1…sn ), ( ¬ s s1…sn)) | (s s1…sn )

g,

where

=

(s s1…sn ), s

}

To illustrate this definition, consider the interference state graph SGg1, ( Ex1) depicted in Fig. 9, where four new edges (given as dashed arrows) are needed to capture the interference Δ = {1} that can occur to the behaviour of Ex1 during a composition under conjunction. Interestingly, it turns out that the interference state graph is able to capture all the possible behaviours a Boolean network can have under composition. As an example, consider Fig. 10 which shows that all paths contained within the projected state graph ( 1 (SG( Ex1, Ex2, g1, g4 ))) are also contained in the interference state graph SGg1, ( Ex1) . This observation is formally proved in the following theorem.

Fig. 7. The labelled state graph 1 (SG( )) derived by applying the left projection 1 to the state graph SG( ) , where = ( Ex1, Ex2, g1, g4 ) (Fig. 4). Note that each node in 1 (SG( )) contains a bottom part representing the original global state and a top part representing the projected state (with state order (gc g2 g3)).

6

BioSystems xxx (xxxx) xxxx

H. Alkhudhayr and J. Steggles

Table 1 The interference sets Δ1 (first Boolean network) and Δ2 (second Boolean network) for each of the 16 binary Boolean operators. A

B

1 F

2 AND

3 AND!

4 A

5 !AND

6 B

7 XOR

8 OR

9 NOR

10 XNOR

11 !B

12 OR!

13 !A

14 !OR

15 NAND

16 T

0 0 1 1 Δ1 Δ2

0 1 0 1

0 0 0 0 {1} {1}

0 0 0 1 {1} {1}

0 0 1 0 {1} {0, 1}

0 0 1 1 {} {0, 1}

0 1 0 0 {0, 1} {1}

0 1 0 1 {0, 1} {}

0 1 1 0 {0, 1} {0, 1}

0 1 1 1 {0} {0}

1 0 0 0 {0, 1} {0, 1}

1 0 0 1 {0, 1} {0, 1}

1 0 1 0 {0, 1} {0, 1}

1 0 1 1 {0} {0, 1}

1 1 0 0 {0, 1} {0, 1}

1 1 0 1 {0, 1} {0}

1 1 1 0 {0, 1} {0, 1}

1 1 1 1 {0} {0}

(st t1i…tmi)

(t i+1 t1i+1…tmi+1)

2

(II)

Then it follows by Definition 2, (I) and (II) that

t i + 1 s1i + 1…sni+ 1 t1i + 1…tmi+ 1)

Si + 1 = (s i + 1 Fig. 8. Example based on Ex1 (Fig. 1) which illustrates how interference Δ= {1} (conjunction) can be modelled by adding additional edges (dashed edge) to a Boolean network's state graph.

We now have two cases to consider with respect to si+1 ⊙ ti+1: (1) Suppose si+1 ⊙ ti+1 = si+1 (i.e. there was no interference for 1). Then by the definition of projection we have 1 (Si + 1)

= (si + 1 s1i + 1…sni+ 1)

Then by (I) we know that 1 (Si )

1 (Si + 1)

(st s1i…sni)

1)) ;

and 2)) .

1 (Si )

g1, 1

1 (Si + 1)

In this section we use the interference state graph to extend the alignment property in order to fully characterise compatibility. The resulting property is referred to as weak alignment and we formally show that it is both sufficient and necessary for compatibility. Recall the entity state projection function gi: S defined for any (s1…sn) S by gi (s1…sn) = si (see Section 4). We can use this projection to extract from a state graph the behaviour of a single entity. As an example, consider Fig. 11 which depicts the labelled state graph Ex1. g1 (SG( Ex1)) representing the behaviour of entity g1 in We formally define weak alignment using the interference graph as follows.

in the interference graph SG g1, 1 ( , 1) . Let Si = (st where st is the state of the merged entity g1 ⊙ g . By definition of projection we know

= (st s1i…sni),and

2 (Si )

= (st t1i…tmi)

Suppose that

(st s1i…sni)

1

(si+1 s1i+1…sni+1)

( ¬ si+1 s1i+1…sni+1)

5.3. Weak alignment

s1i…sni t1i…tmi) 2

1 (Si )

1

g1, 1

to SG g1, 1 ( 1) as required. (ii) Follows along similar lines to (i) above. □

Proof. (i) Let = S1, S2, … Path(SG( )) . Then we need to show there Path(SG g1, 1 ( exists a path 1)) such that 1 ( ) = . It suffices to show that for any i there exists an edge 1

= ( ¬ si + 1 s1i + 1…sni+ 1)

By the definition of the interference set Δ1 this case can happen only when si+1 ∈ Δ1 and so by Definition 7 and (I), we must have added the edge

Theorem 4. Let 1 and 2 be two Boolean networks, and let 1 2 2 1 2 g1 = ( 1, g 2 . Let 1, 2, g , g ) , and let Δ and Δ be the first and second interference sets for ⊙. Then we have:

Path(SG g1, 1 ( Path(SG g 2, 2 (

1 (Si + 1)

and so the result follows since by Definition 7 we know all the edges in SG( 1) are contained in SG g1, 1 ( 1) . (2) Suppose si+1 ⊙ ti+1 ≠ si+1 (i.e. there was interference for 1). Then by the definition of projection we must have

Fig. 9. Example illustrating how the interference state graph SGg1, ( Ex1) is derived from the state graph SG( Ex1) by adding edges (dashed arrows) to model the interference on g1 that can occur during composition using conjunction (Δ = {1}).

(i) Path( 1(SG( ))) (ii) Path( 2 (SG( )))

1

(I)

Fig. 10. Example illustrating that the interference state graph SGg1, ( Ex1) contains all the behaviour in the projected state graph (where Δ = {1} represents the interference associated with conjunction). 7

1 (SG(

(

Ex1,

Ex2,

g1, g4 )))

BioSystems xxx (xxxx) xxxx

H. Alkhudhayr and J. Steggles

(iii) (iv)

1( 1

2)

1( 1

2)

= =

1

1

and and

2( 1

2)

2( 1

2)

2;

= =

2.

and

Proof. (i) Let 1 = S1, S2, … Path(SG( 1)) be an infinite path, and let 2)) such that g1 ( 1) = g 2 ( 2 ) . Let 2 = T1, T2, … Path(SG g 2, 2 ( 1, g 2 ) . In order to show that = ( , , g 1 2 Path(SG( )) it suffices to show that for each i , 1 2

Si

Ti

Si + 1

(I)

Ti + 1

We know by assumption that for any i

Si

Fig. 11. A labelled state graph g1(SG( Ex1)) based on projecting the state of a single entity g1 (where the entity order in the global state is (g1 g2 g3)).

1

Si + 1 and Ti

2

Ti + 1

g 2, 2

Si = (s i s1i…sni) , Si + 1 = (s i + 1 s1i + 1…sni+ 1) , Ti = (t i t1i…tmi ) , Let and Ti + 1 = (t i + 1 t1i + 1…tmi+ 1) . Then by definition of merging states we have Si

t i s1i…sni t1i…tmi)

Ti = (s i

t i + 1 s1i + 1…sni+ 1 t1i + 1…tmi+ 1)

Ti + 1 = (s i + 1

Si + 1

By assumption g1 ( 1) = idempotence of ⊙ we have

si

Ti

1)))

Path(

g2 (SG g 2, 2 (

2))),

Path(

g 2 (SG(

2)))

Path(

g1 (SG g1, 1 (

1)))

(s i

g1 (SG(

Ex1)))

Path(

g4 (SGg4,

(

and

s i+1

s i+1

g

2(

Then we have: (i) 1 2 (ii) 1 2

2 ) and

g 2 ( 2)

Path(SG( Path(SG(

= ( (

g

1(

1,

1,

2,

t i+1 = si+1

(s i + 1

that

ui + 1 s1i + 1…sni+ 1 t1i + 1…tmi+ 1)

ui + 1

t i+1 = si+1

s i+1 = s i+1

ui + 1

i+1

and so s ⊙ ti+1 = si+1 ⊙ ui+1 as required. (ii) Follows along similar lines to (i) above and so for brevity is omitted. (iii) We need to show that 1 ( 1 2 ) = 1 and 2 ( 1 2) = 2 . Let 1 = S1, S2, … Path(SG( 2)) . 1)) , 2 = T1, T2, … Path(SG g 2, 2 ( Then we have 1

2

= S1

For each i g1 ( 1) = have

si

g

T1, S2 2(

T2, …

, let Si = (s i s1i…sni),and Ti = (t i t1i…tmi) . Since we know i i 2 ) , it follows that s = t and so since ⊙ is idempotent we

t i = si = t i

Ti ) = Ti as required. Ti ) = Si and 2 (Si Then it follows that 1 (Si (iv) Follows along similar lines to (iii) above and so again is omitted for brevity.□ We can now formally prove that weak alignment is able to fully characterise compatibility (i.e. weak alignment is a sufficient and necessary property for compatibility).

1)

2,

t i s1i…sni t1i…tmi)

Given that we can infer si+1 = ¬ ui+1 it follows that

Ex3)))

1 Lemma 5. Let 1 and 2 be two Boolean networks such that g 1 2 and g 2 . Let ⊙ be a binary Boolean operator which is idempotent and let Δ1 and Δ2 be its associated interference sets. Let 1 Path(SG( 1)) , Path(SG( Path(SG g1, 1 ( and 2 2)) , 1)) , 1 2 Path(SG g 2, 2 ( 2)) be infinite paths such that

=

(ui + 1, t1i + 1, …, tmi+ 1)

as required. Case 2: Suppose ti+1 ≠ ui+1 (i.e. there was interference for 2 ). It follows by Definition 7 that ui+1 ∈ Δ2 and ti+1 = ¬ ui+1. Furthermore, since ⊙ is idempotent we must have ¬ui+1 ⊙ ui+1 = ¬ ui+1 (i.e. ¬ui+1 is the value that causes the interference). By assumption g1 ( 1) = g 2 ( 2 ) we know that si+1 = ti+1 and so since ⊙ is idempotent we have

(where Δ = {1} for conjunction) since interference on Ex3 provides the missing path 00, 00, …. It therefore follows that Ex3 Ex1 and weakly align on g1 and g4 under ∧. Recall from Section 4 that in the context of a composition 1 2 Tr( we can merge traces and ( 1 1) 1, 2, g , g ) 1 Tr( ) on g and g2 using the given binary Boolean operator 2 2 α1 ⊙ α2. Then this idea can be applied to paths in the obvious way. The following are some useful results that can be shown about merging paths that exhibit the weak alignment property.

g1 ( 1)

2

Therefore in order to prove (I) we have to show that si+1 ⊙ ti+1 = si+1 ⊙ ui+1 which we do by considering the following two cases. Case 1: Suppose ti+1 = ui+1 (i.e. there was no interference for 2 ). Then we must have

Weak alignment is a more expressive property than alignment and is able to capture the fundamental relationship that exists between compatible Boolean networks. To illustrate this, consider again the counter example presented in Fig. 6. In this example we have g1 (Tr( Ex1)) g4 (Tr( Ex3)) g4 (Tr( Ex3)) g1 (Tr( Ex1)) but and so Ex3 do not align on g1 and g4. However, as Fig. 12 Ex1 and shows we do have

Path(

we know that si = ti and so by

t i = si = t i

Then we know by definition of

Definition 8. (Weak alignment) Let 1 and 2 be Boolean 2 networks, and let g1 1 and g 2 . Let ⊙ be a binary Boolean operator, and let Δ1 and Δ2 be the first and second interference sets for 1 2 ⊙. Then we say that 1 and 2 are weakly aligned on g and g under ⊙ iff g1 (SG(

g 2 ( 2)

Now suppose

Fig. 12. State graph for Ex3 (where Ex1 and interference state graph for Path( g1 (SG Δ = {1} for conjunction) which show that ( Path(SGg4, ( Ex1))) Ex3)) (where the global states have entity order (g1 g2 g3) and (g4 g5)).

Path(

,

g1, g 2))) ; g1, g 2))) ; 8

BioSystems xxx (xxxx) xxxx

H. Alkhudhayr and J. Steggles

Theorem 6. Let 1 and 2 be two Boolean networks such that 2 g1 1 and g 2 , and let ⊙ be a binary Boolean operator which 1 2 is idempotent. Then 1 and 2 are compatible on g and g under ⊙ iff 1 2 1 and 2 are weakly aligned on g and g under ⊙.

To completely characterise compatibility, we considered formalising the interference that can occur to the behaviour of a Boolean network in a composition using a state graph approach. We used the resulting interference state graph to define the weak alignment property which we then showed completely characterised compatibility given an idempotent Boolean operator (i.e. it is sufficient and necessary for compatibility). This final result is very important since weak alignment allows us to check behaviour preservation (compatibility) without referencing the composed model and so helps avoid potentially limiting state space explosion issues. The compositional framework presented here is supported by a prototype tool that automates the composition process and associated behaviour preservation analysis.

1 1 2 = ( Proof. Let and Δ2 be the 1, 2, g , g ) and let Δ interference sets associated with ⊙. We have two cases to consider. 1 2 Case 1: Suppose 1 and 2 are compatible on g and g under ⊙. 1 Then we need to prove that 1 and 2 are weakly aligned on g and g2 under ⊙. By the definition of weak alignment (Definition 8) we need to show

Path(

g1 (SG(

1)))

Path(

g2 (SG g 2, 2 (

2)))

(1)

Path(

g 2 (SG(

2)))

Path(

g1 (SG g1, 1 (

1)))

(2)

Path(SG(

1)) .

To show (1) consider any of compatibility we know

Path(SG(

1))

Path(

1 (SG(

)

Path(SG g 2, 2 (

Since we must have g1 (

)=

g2 ( 2 (

Then by the assumption

The main focus of existing work on compositional techniques in Boolean networks is on the efficient identification of attractors. The composition of random Boolean networks is considered in Dubrova and Teslenko (2005) and an approach based on identifying independent components of a Boolean network whose composed behaviour can infer attractors is developed. In Tournier and Chaves (2013) a compositional approach for asynchronous Boolean networks is considered based on using a state graph approach. Various work has also considered extending the applicability of SAT based techniques by partitioning large models and then aggregating the results (see for example Guo et al., 2014; Hong et al., 2015). An algorithmic approach for efficient attractor identification in large-scale Boolean networks based on analysing subnetworks in a model is considered in Zhao et al. (2013), Cheng et al. (2012). An interesting range of research focusing on subnetwork embeddings and asynchronous behaviour preservation in Logical Regulatory Graphs (Thomas and d’Ari, 1990) can be found in the literature. In Siebert (2009) the identification of subnetworks that control the overall Boolean asynchronous behaviour of a model has been considered. By extending existing concepts of singular states and local interaction graphs they are able to identify the key structures influencing the behaviour of the model and so use this to predict its behaviour. The preservation of behaviour is also considered in Bernot and Tahi (2009) where they use the parameters associated with the output edges leaving an embedded network to formulate a necessary and sufficient condition for behaviour preservation. Further interesting work on the embedding of networks and modularity in this setting can be found in Mabrouki et al. (2011), Delaplace et al. (2012). In very recent work, a modular approach for Boolean Automata Networks (BANs) is presented in Perrot et al. (2018) where external inputs are added to the base model along with a notion of wiring. The initial results they present indicate a promising approach to investigating the composition and decomposition of BANs. The approach taken in this paper based on composing synchronous Boolean networks by using binary Boolean operators to merge entities and then characterising behavioural preservation using interference appears to differ significantly from previously published work.

)))

and so there must exist an infinite path = 1( ) . By Theorem 4 we know that 2(

6.1. Related work

Path(SG( )) such that

2))

g1 ( 1 (

)) =

g2 ( 2(

)) it follows that

))

Therefore (1) must hold as required. The proof of (2) follows along similar lines to (1) above and so for brevity is omitted. 1 2 Case 2: Suppose 1 and 2 are weakly aligned on g and g under 1 2 ⊙. Then we need to prove that 1 and 2 are compatible on g and g under ⊙. By the definition of compatibility (Definition 4) we need to show

Path(SG(

1))

Path(

1 (SG(

)))

(3)

Path(SG(

2))

Path(

2 (SG(

)))

(4)

To show (3) consider any path 1 Path(SG( 1)) . Then by assumption and the definition of weak alignment (Definition 8) there must exist a path 2 Path(SG g 2, 2 ( 2)) such that g1 ( 1)

=

g 2 ( 2)

It follows by Lemma 5.(i) That 1

2

Path(SG( ))

Since by Lemma 5.(ii) we know 1 ( 1 2 ) = 1 the result follows. The proof of (4) follows along similar lines to (3) above and so again is omitted for brevity. □ 6. Conclusions In this paper we set out to develop a compositional framework for Boolean networks in order to facilitate the construction and analysis of large scale models. This work was motivated by interesting interactions with the synthetic biology group at Newcastle1 and their search for formal tools and techniques to support their work on engineering biological systems. We have formally defined our compositional approach based on the idea of merging entities between models using binary Boolean operators. In order to formalise the preservation of a Boolean network's behaviour within a composed model we introduced the notion of compatibility. We formulated the alignment property which we showed was a sufficient condition for ensuring compatibility given an idempotent Boolean operator and used it to investigate the composition of duplicate models. 1

6.2. Future work The formal compositional framework developed in this paper provides the basis for a large range of future work. For example, we are currently working to generalise the compositional framework to allow multiple entities and multiple Boolean networks to be composed simultaneously. We have also begun to extend our compositional framework to multi-valued networks (see for example Rudell and Sangiovanni-Vincentelli, 1987; Thieffry and Thomas, 1995; Banks and Steggles, 2012), an extension of the Boolean network approach in which an entity's state is allowed to be a range of discrete values instead of just Boolean. In turn we aim to further strengthen the tool support with these results to ensure the practical application of the techniques

http://www.ncl.ac.uk/csbb/. 9

BioSystems xxx (xxxx) xxxx

H. Alkhudhayr and J. Steggles

developed is fully supported. Finally, we are interested in using our compositional framework as the basis for developing formal techniques for decomposing Boolean networks. The idea is to develop practical techniques and tools to allow the automatic decomposition of realistically large Boolean networks and so facilitate their analysis and investigation.

based on Boolean satisfiability for genetic regulatory networks. PLoS ONE 9 (4), e94258. https://doi.org/10.1371/journal.pone.0094258. Harvey, I., Bossomaier, T., 1997. Time out of joint: attractors in asynchronous random Boolean networks. Proceedings of the Fourth European Conference on Artificial Life. MIT Press, Cambridge, pp. 67–75. Hong, C., Hwang, J., Cho, K.-H., Shin, I., 2015. An efficient steady-state analysis method for large Boolean networks with high maximum node connectivity. PLoS ONE 10 (12), e0145734. https://doi.org/10.1371/journal.pone.0145734. Huang, S., Ingber, D.E., 2000. Shape-dependent control of cell growth, differentiation, and apoptosis: switching between attractors in cell regulatory networks. Exp. Cell Res. 261 (1), 91–103. Kauffman, S.A., 1969. Metabolic stability and epigenesis in randomly constructed genetic nets. J. Theor. Biol. 22 (3), 437–467. Kauffman, S.A., 1993. The Origins of Order: Self Organization and Selection in Evolution. Oxford University Press, USA. Kurgalin, S., Borunov, S., 2018. The Discrete Math Workbook: A Companion Manual for Practical Study, Texts in Computer Science. Springer. Mabrouki, M., Aiguier, M., Comet, J.-P., Gall, P., Richard, A., 2011. Embedding of biological regulatory networks and property preservation. Math. Comput. Sci. 5, 263–288. https://doi.org/10.1007/s11786-011-0092-3. Perrot, K., Perrotin, P., Sené, S., 2018. A framework for (de). In: Durand-Lose, J., Verlan, S. (Eds.), MCU 2018, LNCS 10881. Springer, pp. 121–136. Rosenblueth, D.A., Muñoz, S., Carrillo, M., Azpeitia, E., 2014. Inference of Boolean networks from gene interaction graphs using a sat solver. In: Dediu, A.-H., Martín-Vide, C., Truthe, B. (Eds.), International Conference on Algorithms for Computational Biology, AlCoB 2014, LNCS 8542. Springer, pp. 235–246. Rudell, R., Sangiovanni-Vincentelli, A., 1987. Multiple-valued minimization for pla optimization. IEEE Trans. Comput. Aided Des. CAD 6 (5), 727–750. Saadatpour, A., Albert, R., 2013. Boolean modeling of biological regulatory networks: a methodology tutorial. Methods 62 (1), 3–12. Schaub, M.A., Henzinger, T.A., Fisher, J., Qualitative networks: a symbolic approach to analyze biological signaling networks, BMC Systems Biology 1 (4). Siebert, H., 2009. Deriving behavior of Boolean bioregulatory networks from subnetwork dynamics. Math. Comput. Sci. 2, 421–442. https://doi.org/10.1007/s11786-0080064-4. Steggles, L.J., Banks, R., Shaw, O., Wipat, A., 2007. Qualitatively modelling and analysing genetic regulatory networks: a petri net approach. Bioinformatics 23 (3), 336–343. Thieffry, D., Thomas, R., 1995. Dynamical behaviour of biological regulatory networks – II. Immunity control in bacteriophage lambda. Bull. Math. Biol. 57 (2), 277–297. Thomas, R., d’Ari, R., 1990. Biological Feedback. CRC Press, Boca Raton. Tournier, L., Chaves, M., 2013. Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics. Automatica 49 (4), 884–893. Wuensche, A., 2004. Aggregation algorithm towards large-scale Boolean network analysis. In: Schlosser, G., Wagner, G.P. (Eds.), Modularity in Development and Evolution. University of Chicago Press, pp. 288–311 Ch. 13. Zhao, Y., Kim, J., Filippone, M., 2013. Compositional properties of random Boolean networks. IEEE Trans. Autom. Control 58, 1976–1985.

Acknowledgments We would like to thank Will Peckham for his work on developing tool support for our compositional framework. We also acknowledge the financial support provided by Faculty of Computing and Information Technology, King Abdulaziz University, Rabigh. Finally, we would like to thank the anonymous referees for their very helpful comments and suggestions. References Akutsu, T., Miyano, S., Kuhara, S., et al., 1999. Identification of genetic networks from a small number of gene expression patterns under the Boolean network model. Pac. Symp. Biocomput. 4, 17–28. Alkhudhayr, H., Steggles, J., 2017. A formal framework for composing qualitative models of biological systems. In: Martín-Vide, C., Neruda, R., Vega-Rodríguez, M. (Eds.), Theory and Practice of Natural Computing, TPNC2017, LNCS 10687. Springer, pp. 25–36. Banks, R., Steggles, L.J., 2012. An abstraction theory for qualitative models of biological systems. Theor. Comput. Sci. 431, 207–218. Bartocci, E., Lió, P., 2016. Computational modeling, formal analysis, and tools for systems biology. PLoS Comput. Biol. 12 (1), e1004591. Bernot, G., Tahi, F., 2009. Behaviour preservation of a biological regulatory network when embedded into a larger network. Fundam. Inform. 91, 463–485. https://doi. org/10.3233/FI-2009-0052. Cheng, D., Zhao, Y., Kim, J., Zhao, Y., 2012. Approximation of Boolean networks. Proc. of the 10th World Congress on Intelligent Control and Automation (WCICA) 2280–2285. https://doi.org/10.1109/WCICA.2012.6358254. De Jong, H., 2002. Modeling and simulation of genetic regulatory systems: a literature review. J. Comput. Biol. 9 (1), 67–103. Delaplace, F., Klaudel, H., Melliti, T., Sené, S., 2012. Analysis of modular organisation of interaction networks based on asymptotic dynamics. In: Gilbert, D., Heiner, M. (Eds.), Computational Methods in Systems Biology, LNCS 7605. Springer, pp. 148–165. Dubrova, E., Teslenko, M., 2005. Compositional properties of random Boolean networks. Phys. Rev. E 71, 056116. https://doi.org/10.1103/PhysRevE.71.056116. Guo, W., Yang, G., Wu, W., He, L., Sun, M., 2014. A parallel attractor-finding algorithm

10