L Tournier
email: laurent.tournier@curie.fr

M Chavesa

M Chaves

Uncovering operational interactions in genetic networks using asynchronous Boolean dynamics

Keywords: Boolean networks, Asynchronous transition graph, Model reduction, Apoptosis, NFκB signalling pathway

Biological networks of large dimensions, with their diagram of interactions, are often well represented by a Boolean model with a family of logical rules. The state space of a Boolean model is finite, and its asynchronous dynamics are fully described by a transition graph in the state space. In this context, a model reduction method will be developed for identifying the active or operational interactions responsible for a given dynamic behaviour. The first step in this procedure is the decomposition of the asynchronous transition graph into its strongly connected components, to obtain a "reduced" and hierarchically organized graph of transitions. The second step consists of the identification of a partial graph of interactions and a sub-family of logical rules that remain operational in a given region of the state space. This model reduction method and its usefulness are illustrated by an application to a model of programmed cell death. The method identifies two mechanisms used by the cell to respond to death-receptor stimulation and decide between the survival or apoptotic pathways.

A c c e p t e d m a n u s c r i p t 1. Introduction

Discrete and, in particular, Boolean models have been playing an increasingly important role in the study and analysis of complex biological systems [START_REF] Sánchez | A logical analysis of the Drosophila gap-gene system[END_REF][START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF][START_REF] Chaves | Robustness and fragility of boolean models for genetic regulatory networks[END_REF][START_REF] Saez-Rodriguez | A logical model provides insights into t cell receptor signaling[END_REF][START_REF] Fauré | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF][START_REF] Chaves | Regulation of apoptosis via the NFκB pathway: modeling and analysis[END_REF] (note that even though [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF] deals with piecewise linear differential models, a large part of the analysis of such systems pertains to the discrete framework). Based on the idea that each variable may take values only on a finite set, discrete models offer a very attractive framework for the systematic study of the dynamics of large systems, which may range from a few to hundreds of variables and their interactions.

The discrete modelling approach is highly relevant for many of the currently data acquisition techniques for signalling and genetic regulatory networks (microarrays, fluorescence markers, electrophoretic mobility shift assays, etc.) which involve more qualitative measurements. A discrete system may be expected to give a good idea of the system's dynamics from the available data (for example, on multistationary, stability, or oscillatory behavior). At the same time, Boolean networks provide a measure of the robustness of a system [START_REF] Chaves | Robustness and fragility of boolean models for genetic regulatory networks[END_REF], since this dynamical information is essentially independent of the system's parameters (such as kinetic rates, binding rates, or degradation constants) and depends only on the interconnection structure, and the logical function characterizing each node. For instance, the transition graph of the network indicates how much a given trajectory may be affected by perturbations, or whether the system is capable of maintaining a given dynamical behaviour despite fluctuations in the environment. Indeed, a major advantage of discrete and Boolean modelling is the possibility, for systems of reasonable size, of fully characterizing all qualitative dynamical trajectories of a particular network, based simply on the structure of links and interactions between nodes. This general characterization and "easier" handling of the state space, counterbalance the loss of detailed information on time evolution and (more realistic) continuous concentration changes.

The study of complex systems with many variables frequently raises questions concerning the possibility of simplifying or reducing the system in some way. To simplify the analysis and gain intuition, it is often useful to identify a smaller, easier to analyze, family of variables and interactions that still faithfully describe the original system and exhibit the same overall qualitative dynamics. Likewise, it is often of interest to find out whether different groups of variables are associated with different dynamics [START_REF] Schmidt | Identifying feedback mechanisms behind complex cell behavior[END_REF]. Another related question is whether all interactions operate at all times, or whether different groups of interactions become active or operational at different times, in response to a precise context [START_REF] Kholodenko | Untangling the wires: a strategy to trace functional interactions in signaling and gene networks[END_REF]. Similarly, find-

A c c e p t e d m a n u s c r i p t

ing interactions which prevent a given target function is useful from the point of view of therapeutical interventions, for instance [START_REF] Klamt | A methodology for the structural and functional analysis of signaling and regulatory networks[END_REF][START_REF] Saez-Rodriguez | A logical model provides insights into t cell receptor signaling[END_REF]. These are all challenging problems, and while some model reduction methods exist, they are generally aimed at special classes of systems (a survey of methods used for control systems can be found in [START_REF] Antoulas | A survey of model reduction methods for large-scale systems[END_REF]; a method for identifying the variables responsible for complex cell behavior was proposed in [START_REF] Schmidt | Identifying feedback mechanisms behind complex cell behavior[END_REF]). One of the objectives of this paper is to show that, to some extent, answers to these questions may be obtained through the discrete systems framework and some of its techniques.

In this context, a model reduction technique is proposed in this paper, for the analysis of asynchronous Boolean networks motivated by biological (namely, signalling or genetic) regulatory networks of intermediate dimension (e.g., 8-20 variables). The model reduction technique combines and adapts two methods: the classical decomposition of a graph into its strongly connected components [START_REF] Cormen | Introduction to algorithms[END_REF] and an identification algorithm described in [START_REF] Liang | REVEAL, a general reverse engineering algorithm for inference of genetic network architecture[END_REF][START_REF] Tournier | Etude et modélisation mathématique de réseaux de régulation génétique et métabolique[END_REF]. The first part of the model reduction technique involves the simplification of the asynchronous transition graph into its strongly connected components. These components are then organized into hierarchical levels, such that any given trajectory can only move into the next level in the hierarchy, but never into the previous level. A new "reduced" transition graph is then constructed which describes the transitions between the strongly connected components (Section 3). The second part of the model reduction involves the identification of the "operational" network (active interactions) that is responsible for the dynamics of the system from a given level in the hierarchy (Section 4).

Finally, in Section 5, biological knowledge is incorporated into the Boolean model by associating a (fixed) matrix of transition probabilities to the graph edges, a process which corresponds to generating a Markov chain. Several relevant quantities can then be computed, such as the expected times for convergence to a given attractor.

The methods proposed above are illustrated by an application to an apoptosis (or programmed cell death) network [START_REF] Schliemann | Modelling and experimental validation of TNFα induced pro-and antiapoptotic signalling[END_REF][START_REF] Chaves | Regulation of apoptosis via the NFκB pathway: modeling and analysis[END_REF], with n = 12, as described in Section 2.3. The dynamics of the network in response to death-receptor stimulation is studied, and two core groups of variables and pathways are identified: it is shown that these correspond to two mechanisms responsible for the decision between programmed cell death or cell survival. In addition, associating a transition probabilities' matrix to the apoptosis network allows us to estimate, among other quantities, the probability of cell survival or death upon stimulation of death receptors.

A c c e p t e d m a n u s c r i p t 2. Asynchronous Boolean models of gene regulatory networks

Discrete dynamical systems have been widely studied for decades, as they provide a good mathematical and algorithmic framework to model systems where variables are known or assumed to take values only in a finite set (as opposed to a continuum of values). In particular, the discrete framework has been often applied to model biological (for instance, genetic) regulatory networks [START_REF] Kauffman | The origins of order[END_REF][START_REF] Thomas | Biological feedback[END_REF]. A mathematical discrete model consists of a finite set of variables together with a family of activation functions, which describe the interactions among the variables. Both variables and activation functions are only allowed to take values in a finite set. A common example is that of Boolean variables, which take only the values 1 or 0 (ON or OFF). The interactions among variables can be described by a (finite) directed graph, called interaction graph. This graph, together with the family of activation functions, define the structure of a discrete system. Each variable will evolve according to a given rule, constructed from the interaction graph. In order to describe the dynamics of such a system, a fixed time unit is postulated, and the system is updated at integer multiples of this time unit (between two updates, the system's variables are assumed to be constant). To update the system one must define a strategy that determines the order in which the variables are evaluated over time.

Two updating strategies have been studied in the literature. The first one is the synchronous strategy, where all variables are simultaneously updated at each discrete instant (see [START_REF] Kauffman | The origins of order[END_REF] for an extensive study of this synchronous strategy; see also [START_REF] Shmulevich | Probabilistic boolean networks: a rule-based uncertainty model for gene regulatory networks[END_REF][START_REF] Wuensche | Basins of attraction in network dynamics: a conceptual framework for biomolecular networks[END_REF]). The dynamics resulting from the synchronous strategy presents some nice mathematical properties (mainly, the transition graph is deterministic) that allow one to simulate high-dimensional networks, in order to find statistically relevant types of dynamical behavior [START_REF] Kauffman | The origins of order[END_REF]. However, if one wants to model a given biological system in a more realistic manner, the synchronous updating strategy is quite a strong assumption, implying that all variables are produced or degraded at the same rates. Thus other approaches have been proposed, by developing asynchronous strategies, where the discrete variables are updated in a heterogeneous way over time. For example, the variables may be updated in a randomly generated order, or according to an order based on known properties (for instance, transcription or translation are slower processes than protein-protein binding) [START_REF] Chaves | Robustness and fragility of boolean models for genetic regulatory networks[END_REF]. Discrete networks with asynchronous updating orders are often called Thomas' networks [START_REF] Thomas | Biological feedback[END_REF][START_REF] Thomas | Multistationarity, the basis of cell differentiation and memory. I. Structural conditions of multistationarity and other nontrivial behavior[END_REF][START_REF] Thomas | Multistationarity, the basis of cell differentiation and memory. II. Logical analysis of regulatory networks in terms of feedback circuits[END_REF], and are much better suited to model the dynamical behavior of biological regulatory networks. Before giving more details on the type of asynchronicity that will be used in this paper, we first recall some basic definitions

A c c e p t e d m a n u s c r i p t

about the structure of discrete networks.

Structure of a Boolean network

In this paper, we will consider Boolean networks, where the variables (which represent, for instance, the level of expression of genes, or the level of concentration of different species, such as proteins) can take only two qualitative values. "0" represents a basal level (inhibition -or weak activation-of the transcription of a gene, or absence2 of a biochemical species) and "1" represents a high level (activation of the transcription of a gene, or presence of a species). We note here that there exist more general frameworks, where the variables can take more than two qualitative values (see for instance [START_REF] Sánchez | A logical analysis of the Drosophila gap-gene system[END_REF][START_REF] Chaouiya | Qualitative analysis of regulatory graphs: a computational tool based on a discrete formal framework[END_REF] or [START_REF] Siebert | Temporal constraints in the logical analysis of regulatory networks[END_REF]).

As most of the work on discrete or Boolean gene networks is based on the same mathematical objects (with slightly different definitions), the following part is only a brief summary of the main definitions and notations that will be used in the rest of this paper (for a detailed explanation of these definitions, one can refer to the extensive literature on discrete networks). Let us begin with the definition of the interaction graph, which is the core of the structure of a discrete model.

Definition 1 (Interaction Graph)

The interaction graph of a n-dimensional Boolean network is defined by G = (V, E), where V = {v 1 , . . . , v n } is the set of nodes (each node may represent a biological species) and E ⊂ V × V is the set of directed edges (representing the interactions between these species). The edge (v j , v i ) exists if node v j influences node v i (e.g. v j activates or inhibits v i ).

Each node v i ∈ V has a set of inputs (possibly empty), which are the nodes that influence its evolution:

I(v i ) = v j ∈ V | (v j , v i ) ∈ E ⊂ V.
For each v i ∈ V, the cardinality of the set I(v i ) (the number of its inputs) is often called the connectivity of node v i and is generally denoted by k i . In order to give a complete definition of the structure of the network, we now define the activation functions.

A c c e p t e d m a n u s c r i p t

Definition 2 The structure of a n-dimensional Boolean network is defined by an interaction graph G = (V, E) together with a collection F = { f i : i = 1, . . . , n}, of Boolean functions:

f i : {0, 1} k i -→ {0, 1},
where, for each i ∈ {1, . . . , n}, f i designates the activation function of node v i , and k i its number of inputs.

Let x i denote the Boolean variable associated with node v i . The updated value of x i , denoted by x i is therefore given by:

x i = f i x i 1 , . . . , x i k i , where: v i 1 , . . . , v i k i = I(v i ).
To illustrate these definitions, let us consider a simple 2-dimensional network, given by the following interaction graph and set of rules:

7654 0123 v 1 " " 7654 0123 v 2 g g b b f 1 (x 2 ) = x 2 , f 2 (x 1 , x 2 ) = xor(x 1 , x 2 ) = (x 1 ∧ x 2 ) ∨ (x 1 ∧ x 2 ).
The notations used in this example are classical in Boole's algebra: if x and y denote two Boolean variables, x denotes the negation of x, and x ∧ y, x ∨ y denote, respectively, the product (logical function and) and summation (logical function or) of x and y. The symbol xor denotes the exclusive or.

Remark 1

The three following points emphasize modeling choices, induced by these definitions, that will be made in the rest of the paper.

• Multiple edges are not allowed in the interaction graph. If a particular element of the network v i has two distinct influences on another element v j , they will be represented by a unique edge in G. Nevertheless, these two influences will appear in the logical rule of f j .

• In the study of discrete models of biochemical networks, the arrows of the interaction graph are usually signed, indicating whether the arrow represents an activation ( + → or →) or an inhibition ( -→ or ). It is to be noted that, for a general network given by Def. 2, it is not always possible to associate a sign to each arrow in an unequivocal manner. In the previous example, for instance, the arrows (v 1 , v 2 ) and (v 2 , v 2 ) cannot be signed as the interactions they represent are neither activations nor inhibitions. Nevertheless, in a Boolean network constructed from the description of a particular biological system, most of the interactions can be signed unequivocally. In the rest

A c c e p t e d m a n u s c r i p t

of the paper, activations and inhibitions will thus be labelled (respectively with '+' and '-') while all other interactions will stay unlabelled.

• In the above approach, connectivity is defined prior to the activation functions. This can lead to inaccuracies if the real connectivity (as defined in [START_REF] Tournier | Etude et modélisation mathématique de réseaux de régulation génétique et métabolique[END_REF]) of the function does not match its apparent connectivity. Thus, the function f

(x 1 , x 2 ) = (x 1 ∧ x 2 ) ∨ (x 1 ∧ x 2 )
has an apparent connectivity of 2, whereas its real connectivity is 1 (indeed f (x 1 , x 2 ) = x 1 ). In the terminology of [START_REF] Shmulevich | Probabilistic boolean networks: a rule-based uncertainty model for gene regulatory networks[END_REF], x 2 is called a fictitious (or non essential) variable for function f . This issue becomes particularly important if one wants to identify a network from given data. It will be addressed in more details in Section 4.

Synchronous vs asynchronous dynamics

Consider a n-dimensional network given by N = (V, E, F ) (see Def. 2). The state space of N is the set Ω = {0, 1} n whose cardinality is 2 n . As the state space is finite, one can represent the discrete dynamical behavior of the network with a finite directed graph, called transition graph. In order to define it properly, we need to assign an updating strategy for the network N. Two main approaches have been considered in the literature. The first one is to consider a synchronous update of all nodes at each time (see [START_REF] Kauffman | The origins of order[END_REF]). If the state of the network at time t is given by the Boolean vector: X(t) = (x 1 (t), . . . , x n (t)) ∈ Ω, then the next state X(t + 1) (also called the synchronous successor) is simply:

X(t + 1) = (x 1 (t + 1), . . . , x n (t + 1)) ∈ Ω,
where, for all i, x i (t + 1) = x i (t) (the updated value of x i (t)). In this case, the temporal evolution of the network is autonomous, in the sense that any vector X ∈ Ω has a (unique) successor F(X) = ( f 1 (X), . . . , f n (X)), and that successor is independent of time t. We can then construct a directed transition graph: its set of nodes is Ω and its set of directed edges is defined by the "successor" function. The main property of the synchronous graph is that it is deterministic, i.e. each state has a unique successor. In particular, this property implies that each connected component of the graph contains a unique attractor, and this attractor is either a cycle or a fixed point. More precisely, the connected components are in fact the basins of attraction of their attractor. The synchronous updating strategy is a very strong assumption, not very realistic if one wants to model the dynamical behavior of a given biological system.

A c c e p t e d m a n u s c r i p t

Actually, as discrete interactions are coarse-grained models of sometimes very complex biochemical processes (often implying several biochemical reactions), it is preferable to consider time-dependent, asynchronous updating strategies (historically introduced by R. Thomas [START_REF] Thomas | Biological feedback[END_REF]). In order to give a precise definition of an asynchronous transition graph, we first introduce the following notation:

• For each X ∈ Ω = {0, 1} n , F(X) ∈ Ω designates the synchronous successor of X: ∀i ∈ {1, . . . , n} , F i (X) = x i .
• For each X ∈ Ω and each i ∈ {1, . . . , n}, X i designates the vector:

X i = (x 1 , . . . , x i-1 , x i , x i+1 , . . . , x n ) ∈ Ω. • For each X ∈ Ω, let U(X) = v i ∈ V | x i x i ⊂ V . U(X)
designates the (possibly empty) set of nodes that can actually be updated when the system is in the state X.

• An element X * ∈ Ω is a steady state of the Boolean model if it satisfies:

U(X * ) = ∅, that is x * i = F i (X *
), for all i ∈ {1, . . . , n}. We also state the following assumptions, that we will suppose verified throughout this paper: Assumption 1 At each discrete time t, at most one node is updated (no update means the network is in a steady state).

Assumption 2 Each state X ∈ Ω such that X F(X) has exactly |U(X)| successors.

The first assumption forbids the simultaneous update of several nodes (which is reasonable from the biological point of view), whereas the second one implies that every possible update is taken into account (i.e. if at state X the node v i is liable to change, then that update must be present in the transition graph). With these assumptions, we can now define the (asynchronous) transition graph: Definition 3 (Asynchronous Transition Graph) The asynchronous transition graph of the network N = (V, E, F ) is the directed graph G = (V, E) where the set of nodes V is the state space Ω = {0, 1} n and the set of directed edges E is defined by:

E = X → X i | X ∈ Ω, v i ∈ U(X) .

A c c e p t e d m a n u s c r i p t

In the comparison between synchronous and asynchronous dynamics, a wellknown result is that, provided Assumptions 1 and 2 are satisfied, the asynchronous and synchronous steady states are the same. The proof follows immediately from the definition of steady state.

The non determinism of the asynchronous transition graph is a fundamental property. It allows to consider any possible trajectory implied by the structure of the network. If one wants to study one particular trajectory, then a particular path in the graph has to be chosen, which is equivalent to the choice of a particular updating strategy. Considering biological applications, such a choice will be based on the information available on the system. Unfortunately, the knowledge of the system is often incomplete. An advantage of this framework is the possibility to test and analyze different plausible sets of updating orders, as will be done later in this paper (it corresponds to the notion of priority classes developed in [START_REF] Fauré | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF]). The main advantage of the asynchronous graph is that it comprises all the possible choices in a finite structure, which allows to find general dynamical properties valid whatever the updating strategy. Obviously, although finite, the size of the transition graph grows exponentially with the dimension of the system (in the Boolean case, its size is exactly 2 n ). This limits the use of general graph algorithms to relatively low dimensional systems (on the order of n = 10-20), with respect to the synchronous case, where the dimension of the system under study can be higher [START_REF] Wuensche | Basins of attraction in network dynamics: a conceptual framework for biomolecular networks[END_REF].

Working example: an apoptosis signalling pathway

The model reduction method will be illustrated by application to an apoptosis network (Fig. 1). Apoptosis, or programmed cell death, is a physiological process which allows an organism to remove damaged or unwanted cells in a "clean" and natural way. The signalling pathways leading to apoptosis play fundamental roles in embryonic development and in adult organisms, by maintaining normal cellular homeostasis in organs and other cellular tissues [START_REF] Danial | Cell death: critical control points[END_REF]. Malfunctioning apoptotic pathways may lead to various diseases, such as cancer (in this case cells do not die, there is insufficient apoptosis), or immunodeficiency and infertility (in this case too many cells die, there is too much apoptosis) [START_REF] Danial | Cell death: critical control points[END_REF].

The apoptosis signalling pathway to be considered in this paper (Fig. 1) is based on the model presented in [START_REF] Chaves | Regulation of apoptosis via the NFκB pathway: modeling and analysis[END_REF], which is, in fact, a discrete version of a continuous model of apoptosis first developed in [START_REF] Schliemann | Modelling and experimental validation of TNFα induced pro-and antiapoptotic signalling[END_REF]. A brief description of the network is provided next, and the reader is referred to [START_REF] Chaves | Regulation of apoptosis via the NFκB pathway: modeling and analysis[END_REF][START_REF] Schliemann | Modelling and experimental validation of TNFα induced pro-and antiapoptotic signalling[END_REF] and references therein for more details. The network is composed essentially of a pro-apoptotic and an anti-apoptotic pathway, which are activated by the same signal: stimulation of death receptors by a factor such as Tumor Necrosis Factor α (denoted TNF in Fig. 1). The proapoptotic pathway is based on the model developed in [START_REF] Eissing | Bistability analysis of a caspase activation model for receptorinduced apoptosis[END_REF], and consists of a family of proteins called caspases, represented by active caspases 3 and 8 (resp., C3a and C8a) in Fig. 1. The caspases play the main role in apoptosis, as they cleave (or break into small pieces) the principal proteins in the cell, eventually leading to a "clean disposal" of the cell in response to an apoptotic signal. The anti-apoptotic pathway is based on the pioneering work of [START_REF] Hoffmann | The IkB-NFkB signaling module: temporal control and selective gene activation[END_REF] and on the models developed in [START_REF] Hoffmann | The IkB-NFkB signaling module: temporal control and selective gene activation[END_REF][START_REF] Lipniacki | Mathematical model of NFkB regulatory module[END_REF], and represents the Nuclear Factor κB (NFκB) signalling pathway. The links among the two pathways are as yet, not fully characterized. It is well known that the NFκB signalling pathway is responsible for activating transcription of both pro-and anti-apoptotic genes [START_REF] Perkins | Integrating cell-signalling pathways with NF-κB and IKK function[END_REF], and thus plays an important role in regulating apoptosis. There is also evidence to show that caspase 3 has an inhibitory influence on inhibitor of IkB kinases (IKK) [START_REF] Frelin | Inhibition of the NF-κB survival pathway via caspase-dependent cleavage of the IKK complex scaffold protein and NF-κB essential modulator NEMO[END_REF]. The components of the NFκB pathway are as follows (in biological terminology): Nuclear Factor κB in the cytoplasm (NFκB) and in the nucleus (NFκBnuc); inhibitor of NFκB (IκB); inhibitor of IκB kinases (IKKa); inhibitor of apoptosis proteins (IAP); caspase-8 and -10-associated RING proteins (CARP); a protein associated with inhibition of complex T2 (FLIP); and a protein regulating IKK activity (A20a).

Binding of TNF to a death receptor activates the anti-apoptotic pathway and, after a certain delay (upon formation of a second complex, denoted by T2), the

A c c e p t e d m a n u s c r i p t

pro-apoptotic pathway is also activated. The anti-apoptotic pathway activates synthesis of various proteins (IAP, CARP, FLIP) that will contribute to inhibit and regulate the caspases. Therefore, TNF stimulation triggers two opposite effects: activation and inhibition of caspases. The dynamics of the pro-and anti-apoptotic pathways, as well as the interconnections between them, will ultimately lead to a decision between cell death or cell survival. An abundance of active caspases (such as C3a) together with a low concentration of IAP typically leads to cell death. In contrast, a high concentration of IAP and a low level of active caspases typically characterizes a living cell (in this case, enough molecules IAP are present to down-regulate the level of active caspases). In the network represented in Fig. 1, and in particular in its corresponding Boolean model (Table 1), steady states corresponding to cell death should satisfy C3a = 1 and IAP = 0, while steady states corresponding to cell survival should satisfy C3a = 0 and IAP = 1.

The Boolean model in Table 1 has been slightly simplified from that in [START_REF] Chaves | Regulation of apoptosis via the NFκB pathway: modeling and analysis[END_REF], namely the mRNAs have been removed and only the corresponding proteins nodes are represented. This does not affect the overall dynamics, but reduces the number of variables to facilitate the use of the asynchronous algorithms. The analysis of the transition graph of the Boolean network will allow us to study the dynamics of the system, in particular the effect of the structure of the network in creating and/or maintaining a balance between the pro-and anti-apoptotic pathways, and ultimately the decision between death or survival.

Hierarchical organization of the asynchronous transition graph

In this section, a general methodology to analyze the asynchronous transition graph of a Boolean network is presented. This methodology is based on different algorithms that are classical in the field of graph theory (mainly the strongly connected components decomposition and the topological sort). One can refer to [START_REF] Cormen | Introduction to algorithms[END_REF] for a detailed analysis of these algorithms. Various implementations of strongly connected components decomposition and hierarchical organization, dedicated to biological networks, already exist in the literature (see, for instance, [START_REF] Fauré | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF][START_REF] Chaouiya | Qualitative analysis of regulatory graphs: a computational tool based on a discrete formal framework[END_REF]). We provide an alternative implementation of these algorithms, specifically designed to handle Boolean models of biological genetic regulatory networks (where a node has only two values). Our goal is to make an optimal use of the hierarchical organization of the transition graph in order to answer, algorithmically, some specific biological questions, related to local dynamical behaviors of the system. For instance, given any state, what regions of the phase space lead, or not, to it? What attractors are reachable, or not, from it? Moreover, in Section 4, we will also

A c c e p t e d m a n u s c r i p t

Node Boolean rule TNF TNF (input of the whole system)

T2 TNF ∧ FLIP IKKa TNF ∧ A20a ∧ C3a NFκB IκB NFκB nuc NFκB ∧ IκB IκB [TNF ∧ (NFκB nuc ∧ IKKa)] ∨ [TNF ∧ (NFκB nuc ∨ IKKa)] A20a TNF ∧ NFκB nuc IAP [TNF ∧ (NFκB nuc ∧ C3a)] ∨ [TNF ∧ (NFκB nuc ∨ C3a)] FLIP NFκB nuc C3a IAP ∧ C8a C8a CARP ∧ (C3a ∨ T2) CARP [TNF ∧ (NFκB nuc ∧ C3a)] ∨ [TNF ∧ (NFκB nuc ∨ C3a)]
Table 1: Boolean rules for the apoptosis network depicted in Fig. 1. See explanation of variables in the text. Note that the variable TNF can be considered an input of the system, as its activation function is TNF = TNF.

provide an algorithmic method, based on the hierarchical organization, to further analyse the system and identify a subnetwork of interactions responsible for such local behaviors.

SCC decomposition and hierarchical organization

The notion of hierarchical organization of a directed graph (or digraph) relies on the well known strongly connected components (SCC) decomposition algorithm. Let us first recall some basics about digraphs (see [START_REF] Cormen | Introduction to algorithms[END_REF] for more details). Let G = (V, E) be a digraph. Two vertices u, v ∈ V are mutually reachable (denoted u ∼ v) if and only if there exist two (directed) paths ρ and ρ such that ρ joins u to v and ρ joins v to u. This relation is clearly an equivalence relation on the set V of vertices. The strongly connected components of the digraph G are then defined as the elements of V/ ∼, that is to say the equivalence classes of the relation ∼.

In other words, a strongly connected component of G is a maximal set of vertices C ⊆ V such that for every pair u, v ∈ C, u and v are reachable from each other.

The SCC decomposition of a digraph G consists in computing the strongly connected components of G: C 1 , . . . , C p and then to compute the digraph G scc = (V scc , E scc ) defined as follows:

A c c e p t e d m a n u s c r i p t

• V scc = C 1 , . . . C p ,
• given 1 ≤ i, j ≤ p, the directed edge (C i , C j ) belongs to E scc if and only if there are u ∈ C i and v ∈ C j such that (u, v) ∈ E.

It can be easily proved (see [START_REF] Cormen | Introduction to algorithms[END_REF]) that the digraph G scc contains no (oriented) cycles.

It is called a dag (for directed acyclic graph). This is a key property of G scc , because every dag can be topologically sorted (see [START_REF] Cormen | Introduction to algorithms[END_REF], section 22.4). A topological sort of a dag can be viewed as a classification of its vertices in several hierarchical levels H 1 , H 2 , . . . , defined, by induction, as follows:

• H 1 = C q ∈ V scc : ∀C q ∈ V scc , (C q , C q ) E scc , • H i = C q ∈ V scc : ∀C q ∈ V scc , (C q , C q ) ∈ E scc ⇒ ∃ j < i, C q ∈ H j ,
so that the vertices of the first level H 1 are vertices with no predecessors, and the predecessors of vertices of level H i , i > 0, are contained in inferior levels H j with j < i (see also Fig. 2). The decomposition and hierarchical organization of a digraph G can be computed in linear time with respect to the number of vertices and edges of G [START_REF] Cormen | Introduction to algorithms[END_REF]. The main interest of this hierarchical organization, applied to the asynchronous transition graph of a Boolean network, is that, whatever path we choose in the graph (i.e. whatever updating order we choose for the variables), once the path leaves a hierarchical level H i , it cannot return to this level. So, any path will travel "down" the hierarchical levels:

H i 1 → H i 2 → . . . (with i 1 < i 2 < . . . ).
Due to this property, we can now give a precise definition of the term attractor for a Boolean network evolving according to an asynchronous strategy.

Definition 4 (Attractor) Let N be a Boolean network. An SCC c * ∈ V scc that has no successor in G scc is called an (asynchronous) attractor of N.

In graph theory, such SCCs are often called terminal SCCs. In other words, the asynchronous attractors of a Boolean network are the strongly connected components of the transition graph that cannot be escaped by the system, whatever the updating strategy. However, it should be noted that it is still possible to construct specific asynchronous updating strategies such that the system gets "stuck" in a non terminal SCC. Indeed, for any SCC that contains at least two states, it is obvious that we can find a particular strategy that allows the system to remain indefinitely in this component (by strong connectedness). Such intermediate SCCs will not be considered as attractors in this paper, as we seek general dynamical properties that are valid for all the choices of updating rules. The hierarchical for all γ ∈ P do 6:

A := A ∪ ATTR(γ) 7:
end for 8: end if 9: return A end for 8: end if 9: return R The functions predecessors and successors return, respectively, the -possibly empty-sets of immediate predecessors and successors of a node in the dag G scc . As explained in the text, if the node γ ∈ V scc belongs to a hierarchical level H i , then its predecessors (resp. its successors) can only lie in hierarchical levels H j with j < i (resp. with j > i).

Definition 5 Let c ∈ V scc be a SCC of the asynchronous transition graph. The sets A(c) and R(c) computed by Algorithm 1 are, respectively, the attraction set of c (i.e. the set of all SCCs that can lead to c) and the reachability set of c (i.e. the set of all SCCs that can be reached from c). If c is an attractor of the network (in other words

, if R(c) = {c}), the set A(c) is its basin of attraction.
The definition of attraction or reachability set should be understood as a weak notion, in the sense that the existence of a trajectory from a SCC c to a SCC c is sufficient to have c ∈ A(c ) (or c ∈ R(c)). On the contrary, Definition 4 of an attractor c * is strong, in the sense that no trajectories are allowed to move out of c * . These notions of strong and weak attraction are remniscent of those found in [START_REF] Casey | Piecewise linear models of genetic regulatory networks: equilibria and their stability[END_REF], in the context of equilibria of differential inclusions. In [START_REF] Casey | Piecewise linear models of genetic regulatory networks: equilibria and their stability[END_REF] a set

A c c e p t e d m a n u s c r i p t

of equilibria is said to be weakly asymptotically stable if there is at least one solution of the differential inclusion for which the set is asymptotically stable in the classical sense. The strong notion holds if the set is asymptotically stable in the classical sense for every solution of the differential inclusion. In a discrete graph, if a 1 , . . . , a r designate the attractors of the network, an element c ∈ A(a 1 ) may lead to attractor a 1 (i.e., there exists an updating strategy such that c leads to a 1 ). If one wants the basin of attraction of a 1 in a strong sense, that is, the set of SCCs that always lead to a 1 (whatever the updating order), then one has to compute the set:

A s (a 1 ) = A(a 1 )\ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ r i=2 A(a i ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,
which may in some cases be reduced to the singleton {a 1 }.

Application to the apoptosis network

The SCC decomposition and hierarchical organization were applied to the NFκB signalling pathway described in Section 2.3. The results presented here were obtained with codes implemented in Matlab. Following the matlab bgl3 library specifications, the graphs are represented with sparse matrices, which allow a quite efficient implementation.

We recall that the system under study is of dimension n = 12, and that one particular variable, TNF, is an input (i.e. its activation function is TNF = TNF). The state space is Ω = {0, 1} n , and the size of the asynchronous transition graph G is 2 n = 4096. The number of strongly connected components is p = 1472, therefore the size of the graph G scc is only 40% of the size of G. After the hierarchical organization of this graph, we found only 38 hierarchical levels, and 3 attractors. Fig. 2 represents a scheme of this graph with its main elements. The fact that TNF is an input implies that its value remains constant, whatever path is chosen. Mathematically, this means that G scc is the union of two disconnected (sub)graphs, denoted T 0 (where TNF = 0) and T 1 (where TNF = 1). The components T 0 and T 1 are two subsets of nodes of G scc that are completely separated (there exist no directed edge going from a SCC in T 0 to a SCC in T 1 , and vice versa).

We found three terminal SCCs in our graph, which means that the system has three different attractors. Using the SCC labels returned by the hierarchization, the SCCs 1 and 123 contain only one state (they are therefore steady states) whereas A c c e p t the third one, 49 is a more complex SCC with 56 states. Table 2 indicates the Boolean values taken by the variables within each attractor. The two steady states belong to T 0 (TNF is absent), the first one corresponds to survival of the cell (the caspases C3a and C8a are absent) and the second one corresponds to the triggering of apoptosis (with activation of the caspases). The complex attractor (SCC 49) belongs to T 1 (TNF is present). As we can see in Table 2, within this attractor the caspases are activated while NFκB, IκB and other factors oscillate. At first, this might seem to indicate that apoptosis will be the final outcome, but, as will be seen later, upon TNF removal, the cell may still choose either the survival or apoptotic steady states.

In the previous example, the computations of the asynchronous transition graph and of its hierarchical organization take only a few seconds (simulations made on a PC: Intel 2.13 GHz, RAM: 1.99 Go). The implementation of these graphs is based on sparse matrices, which makes the storage of data and the run of the different algorithms rather efficient (including the computation of attraction and reachability sets). However, as the size of the transition graph is 2 n , the time complexity of its construction grows exponentially. Moreover, the non determinism of the graph makes its analysis more difficult than in the synchronous case (in particular in the search of attractors, or attraction basins). At present, with not fully optimized codes, the computation time remains reasonable for dimensions n around 15. This is a major difference with the synchronous framework, where the determinism of the transition graph makes it possible to analyze much higher dimensional systems [START_REF] Wuensche | Basins of attraction in network dynamics: a conceptual framework for biomolecular networks[END_REF].

The hierarchical organization has the advantage of characterizing the dynamics for all possible updating orders. It allows for instance to detect "spurious" behaviors, that may appear when the network get apparently stuck in a non terminal SCC. It also comprises all the possible dynamical behaviors in a finite structure, avoiding generation of large numbers of simulations.

Identification of operational interactions

The family of SCCs and their transition graph describe a new state space (reduced from 2 n to p states), which characterizes the dynamics of a new, reduced, system. A second stage in our model reduction procedure involves the determination of a set of rules governing this new system. Starting from the simplified asynchronous transition graph G scc , the goal of this part is to reconstruct, as far as possible, the active, or operational interactions along time (the definition of an

A c c e p t e d m a n u s c r i p t

operational interaction can be found in the following). Recall that the SCCs are hierarchically organized, in such way that trajectories can go in one sense from one level to another. Thus, the G scc graph is particularly well suited for an identification process which consists, roughly, of finding all the operational interactions from a level l down to the terminal level (i.e. down to the possible attractors).

In general, this identification method can be applied to any, hierarchically organized, state space. It can be used to uncover groups of variables (and interactions) that are mainly responsible for the dynamical behavior of the system in a given region of the state space. More precisely, the method identifies groups of interactions responsible for the system's asymptotic behaviour, from a given level in the state space, or within a "self-contained" subgraph.

Synchronous identification algorithm

The identification technique that will be used consists in the Boolean identification algorithm REVEAL developed in [START_REF] Liang | REVEAL, a general reverse engineering algorithm for inference of genetic network architecture[END_REF]. The term identification must be understood in a precise sense, related to the Boolean structure of the networks under study. Basically, given a family of transition pairs {s i → s j , s i , s j ∈ Ω}, one wants to reconstruct the structure of a network, that is: its interaction graph G and (possibly) its set of Boolean rules F = { f i , i = 1 . . . n}. This algorithm is based on the fact that data (typically, time series issued from DNA microarrays), are supposed to be synchronous. More precisely, data consist in temporal sequences of Boolean vectors, X t t=0,1,... , that satisfy, for all t and i, X t+1 i = f i (X t ), where the f i are the functions that we seek to identify. Details about REVEAL can be found in [START_REF] Liang | REVEAL, a general reverse engineering algorithm for inference of genetic network architecture[END_REF][START_REF] Tournier | Etude et modélisation mathématique de réseaux de régulation génétique et métabolique[END_REF] and references therein. In order to describe our method, we summarize in the following some basic facts. Definition 6 Let f : {0, 1} n → {0, 1} and j ∈ {1, . . . , n}. The variable x j will be called fictitious in function f if, for all (x 1 , . . . , x j-1 , x j+1 , . . . ,

x n ) ∈ {0, 1} n-1 , f (x 1 , . . . , x j-1 , 0, x j+1 , . . . , x n ) = f (x 1 , . . . , x j-1 , 1, x j+1 , . . . , x n )
The set of indices j such that x j is non-fictitious in f is called the support of f and is denoted σ( f ).

The real connectivity of f , as evoked in Rem. 1, is in fact the size of σ( f ). We will also use the following notation. For 1 ≤ k ≤ n, let I = {i 1 , . . . , i k } ⊆ {1, . . . , n} be a set of indices and let Ī = { j 1 , . . . , j n-k } denote the set {1, . . . , n}\I. For each x = (x 1 , . . . , x n ) ∈ {0, 1} n , x I (respectively xĪ) designates the k-tuple (x i 1 , . . . , x i k )

A c c e p t e d m a n u s c r i p t

(respectively the (nk)-tuple (x j 1 , . . . , x j n-k )). If x I is fixed in {0, 1} k , f x I will designate the function:

f x I : {0, 1} n-k -→ {0, 1} (x j 1 , . . . , x j n-k ) -→ f (x 1 , . . . , x n ).
This lemma will be useful in the following:

Lemma 1 The support of f satisfies σ( f ) ⊆ I if and only if, for all x I ∈ {0, 1} k , the function f x I is constant.

The proof, that comes directly from Def. 6, is left to the reader.

The data used by REVEAL, given by a series of (input state,output) pairs, are first arranged under the form of a truth table, defined as a couple (In, Out) of q × n Boolean matrices such that the rows of In are all distinct (1 ≤ q ≤ 2 n ). The purpose of REVEAL is to find sets of Boolean functions that are consistent with (In, Out), as defined by: Definition 7 Let In and Out be two q × n Boolean matrices, such that all the rows of In are different. A Boolean function f j : {0, 1} n → {0, 1} is consistent with (In, Out) with respect to column j if, for all row R i (In) of In (i = 1 . . . q), we have f j (R i (In)) = Out i j . A function F = ( f 1 , . . . , f n ) : {0, 1} n → {0, 1} n is consistent with (In, Out) if for all j, f j is consistent with respect to column j.

Basically, the main loop of REVEAL consists in finding, for each node v j of the network, k-tuples of nodes (v i 1 , . . . , v i k ) that are possible sets of inputs for the logical rule of v j . Mathematically, the existence of such a tuple implies the existence of a function f j that is consistent with respect to column j and which satisfies σ( f j ) = {i 1 , . . . , i k }. To test whether the tuple (v i 1 , . . . , v i k ) is a possible set of inputs for v j , the following equality is checked:

H (C Out ( j), C In (i 1 ), . . . , C In (i k )) = H (C In (i 1 ), . . . , C In (i k )) ,
where C X (i) designates the i-th column of matrix X, and H designates the Shannon's entropy (a classical notion in information theory), defined by: Definition 8 Let M 1 , . . . , M r ∈ {0, 1} q , and let M designate the q × r Boolean matrix of columns M 1 , . . . , M r . The entropy of M is the quantity H(M) defined by:

H(M) = H(M 1 , . . . , M r ) = - x∈{0,1} r P x M log(P x M ),

A c c e p t e d m a n u s c r i p t

x i l , for all l = 1, . . . , k. Then, set f j (x) = Out i j (which, according to (iii), is independent of the choice of i). If, on the other hand, the pattern x I is not present in In, then we have a degree of freedom: we choose (x I ) ∈ {0, 1} and set f j (x) = (x I ). By construction, such a function f j is consistent with (In, Out) with respect to column j. Furthermore, each projection ( f j ) x I is constant which, according to Lemma 1, implies that σ( f j ) ⊆ I. This proves that (i)⇒(ii). The converse implication can be proved using analogous arguments.

The main loop of the algorithm is made for increasing k (from k = 0, detecting constant functions, to k = n). In order to avoid redundancies, the algorithm stops once it finds the minimal connectivity of each node, i.e. the minimal value of k that ensures the existence of a consistent function. Therefore, the functions identified by REVEAL are minimally consistent, in the sense that, among all consistent functions, they are the ones with a support of minimal size.

As a general reverse engineering algorithm, the purpose of REVEAL is to identify consistent functions "from scratch", that is, independently of any prior knowledge about the system under study. The problem that is addressed here is slightly different, as we are searching subsets of interactions included in an already known set of interactions. Therefore, given a collection of Boolean functions F = { f 1 , . . . , f n }, we force the algorithm to select only the minimally consistent functions g j whose support is included in σ( f j ) (see Alg. 2).

Algorithm 2 -Identification of minimally consistent sub-networks.

Input: (G, F ): Boolean network. In, Out: Boolean q × n matrices. Output: set of sub-networks minimally consistent with (In, Out).

1: mark all nodes {v 1 , . . . , v n } as untreated 2: for k = 0, 1, . . . , n do 3: for all untreated node v i do 4:

for all (i 1 , . . . , i k ) ⊆ σ( f i ) do 5: if H (C Out (i), C In (i 1 ), . . . , C In (i k )) = H (C In (i 1 ), . . . , C In (i k )) then 6: → (v i 1 , . . . , v i k ) are inputs of node v i 7:
→ mark node v i as treated

A c c e p t e d m a n u s c r i p t

For the original REVEAL algorithm, the existence of a minimally consistent function is guaranteed, but this is not the case for Algorithm 2. In Section 4.2, we propose a way to construct truth tables (In, Out) that will ensure the existence of solutions. Regarding uniqueness, the proof of Prop. 1 shows that if a solution exists, it is not necessarily unique if the size of the table q is < 2 n .

Algorithmic search for operational interactions

In order to properly define the algorithmic search for operational interactions, let us first recall the following definitions, classical in graph theory.

Definition 9 Let G = (V, E) be a directed graph. • For V ⊂ V, the subgraph of G, induced by V , is the graph G(V ) = (V , E(V ))
, where E(V ) is the subset of E containing only the directed edges whose head and tail both belong to V .

•

For E ⊂ E, the graph G = (V, E ) is called a partial graph of G (G is said to be included in G, denoted by G ⊂ G).
Let N = (G, F ) be a (given) n-dimensional Boolean network, and let G = (Ω, E), G scc = (V scc , E scc ) denote, respectively, its asynchronous transition graph and its (hierarchically organized) SCC decomposition. We will also make use of the synchronous transition graph of the system, which will be designated by G sync = (Ω, E sync ) (recall that each state s ∈ Ω has a unique successor in that graph, which is the state s = (s 1 , . . . , s n ) where s i = f i (s)). If c ∈ V scc denotes a SCC of G, Algorithm 1 computes the reachability set of c, R(c). This set contains all SCCs (and thus all states) that are reachable by the system starting from c, whatever the updating order of the variables. Therefore, it is possible to construct the subgraph of G scc (respectively, the subgraph of G), induced by R(c), that contains all possible SCC trajectories starting from c (resp., all possible state trajectories starting from any state in c). Let G scc (c) (resp., G(c)) denote this subgraph. Using Assumptions 1 and 2, it is straightforward to deduce, from G(c), the subgraph G sync (c) of G sync induced by R(c). This subgraph G sync (c) is then organized under the form of a truth table (In, Out) (see Algorithm 3 -step 1). The application of Algorithm 2 on these synchronous data (step 2 of Alg. 3) allows identification of all interactions which are sufficient to reproduce all asynchronous trajectories issued from c. These interactions are called operational and can be formally defined as follows: • An edge e ∈ E is a non-operational interaction associated with c if the asynchronous subgraph Ĝ(c), generated by Ĝ = (V, E \ {e}) (starting from c), satisfies: Ĝ(c) = G(c).

• An edge e ∈ E is an operational interaction associated with c if e is not non-operational.

Let E c denote the set of operational interactions, the interaction graph G c = (V, E c ) will be called operational graph.

In other words, a minimal family of operational interactions (for the transition graph generated by a set of states c) contains as small a set as possible of the original interactions that still generates the original graph G(c). The term minimal refers to the fact that, thanks to REVEAL, the supports of the identified functions are of minimal size.

From the construction of the table (In, Out), it is obvious that Alg. 2 effectively returns a solution (in the worst case, it will return the original functions f 1 , . . . , f n , that are obviously consistent with (In, Out)). However, as already evoked, the returned solution is not necessarily unique. Nevertheless, for the apoptosis network treated in this paper, the operational graph is in fact unique for any SCC c (see Section 4.3). As the reachability set R(c) captures all possible asynchronous successors of c, it is easy to see that, if we successively compute the operational graphs along a particular SCC trajectory: (c 1 , c 2 , . . . , c l ) (where c l is an attractor of the system), then we have the following inclusions:

G ⊃ G c 1 ⊃ G c 2 ⊃ • • • ⊃ G c l .
(Symbol ⊂ designates the graph inclusion defined in Def. 9). In other words, this means that it is possible to visualize, along a trajectory, at which step an interaction (represented by a directed edge of G) may become non-operational. Ultimately, the final graph G c l comprises interactions that remain always operational through the whole trajectory (up to the attractor).

A related notion is that of minimal cut sets (MCS) for logical interaction graphs, suggested and used in [START_REF] Klamt | A methodology for the structural and functional analysis of signaling and regulatory networks[END_REF]. A MCS has been defined with respect to a

A c c e p t e d m a n u s c r i p t

Algorithm 3 -Identification of asynchronous operational interactions.

Step 1: Construction of synchronous data from a subgraph Γ of the asynchronous transition graph G.

Input: Γ = (V(Γ), E(Γ)): subgraph of G. Output: In, Out: Boolean matrices of size q × n (partial truth tables).

1: q := 1 2: for all state X ∈ V(Γ) do 3:

In(q, :) := X /* fill in the q-th row of In */ for all asynchronous successors Y of X (in Γ) do 6:

find i such that Y = X i 7:

S (i) := not(X(i))

8:

end for

9:

Out(q, :) := S , q := q + 1 10: end for Step 2: Identification of operational interactions from a SCC c. certain target function or response, and by analogy with the analysis of (continuous, stoichiometric) metabolic networks. A MCS is a minimal set of reactions whose removal will prevent the target response. The notion of operational interaction defined here differs from that of MCS. Mainly, the operational graph identifies a subnetwork responsible for a certain dynamical behaviour, that is comprised in a certain region of the transition graph. The reduction of the whole state space to this particular region allows to ensure all non-operational interactions to be removed. In contrast, MCS are not defined in terms of transition graphs.

Application to the apoptosis network

In this section, Algorithm 3 is applied to the apoptosis NFκB signalling pathway. This system is interesting because it exhibits two global dynamical properties that are often studied in systems biology. The first one is the multistationarity,

A c c e p t e d m a n u s c r i p t

i.e. the coexistence of several attractors. According to Fig. 2, this happens when TNF is absent (the two attractors are steady states). The second one is the presence of oscillations, which appear when TNF is activated. These two general properties are often related to the interaction graph, and in particular to the presence of positive and negative feedback loops. Two famous conjectures, stated in the eighties by R. Thomas [START_REF] Thomas | Biological feedback[END_REF], have been proved in different mathematical frameworks (see for instance [START_REF] Gouzé | Positive and negative circuits in dynamical systems[END_REF][START_REF] Soulé | Graphic requirements for multistationarity[END_REF][START_REF] Remy | Graphic requirements for multistability and attractive cycles in a Boolean dynamical framework[END_REF]). The first one states that the presence of positive feedback loops is a necessary condition for multistationarity. The second one states that the presence of negative feedback loops is a necessary condition for the presence of oscillations (in continuous frameworks, oscillations may be damped, and can then be related to the biological concept of homeostasis). In the following, Alg. 3 is used to identify which feedback loops of Fig. 1 are effectively responsible for the presence of oscillations and for multistationarity.

According to the analysis presented in Section 3.2, the state space of the apoptosis NFκB system can be separated in two regions, T 0 and T 1 , according to the value of the input TNF. Within region T 1 (TNF = 1), we found a unique attractor, a 0 (see Fig. 2). It contains 56 states, and the variables that can oscillate within a 0 are indicated in Table 2. As it is the only attractor present in T 1 , we know that all trajectories starting from any state of T 1 will eventually reach it and remain in it for all subsequent times. When applying Algorithm 3 to a 0 , we obtain the operational graph G a 0 depicted in Figure 3. . The isolated variables have a fixed value, that can be found in Table 2.

This interaction graph is included in the initial one (Fig. 1), and comprises all interactions that remain active in a 0 . As we can see in Fig. 3, six of the twelve

A c c e p t e d m a n u s c r i p t

variables are isolated, that is, they keep a constant value over time. Their values can be found in Table 2 (in particular, the caspases C3a and C8a are activated). Among the six remaining variables, only three are involved in feedback loops: IκB, NFκB and NFκB nuc (the three other variables are affected by those three, but do not affect them). One can see that both feedback loops: NFκB nuc IκB and IκB → NFκB → NFκB nuc → IκB are negative, as they contain an odd number of inhibitions. Thomas' conjectures [START_REF] Gouzé | Positive and negative circuits in dynamical systems[END_REF][START_REF] Soulé | Graphic requirements for multistationarity[END_REF][START_REF] Remy | Graphic requirements for multistability and attractive cycles in a Boolean dynamical framework[END_REF] lead to the following conclusions: (i) there cannot be more than one steady state because the graph has no positive loop, and (ii) oscillations are possible because the graph has negative loops. This is indeed what is observed in the system. As these loops are the only ones that are active in attractor a 0 , it can be infered that they are the ones responsible for the presence of oscillations in the system.

Let us now consider the region T 0 , where TNF = 0. As we can see in Figure 2, this region contains two attractors, corresponding to two steady states. The first steady state (a L ) corresponds to the "survival" of the cell, with the inhibition of the caspases, whereas the second one (a D ) corresponds to the triggering of apoptosis, with activation of the caspases. In order to find the feedback loop that is responsible for the coexistence of these two steady states, we have to find a SCC c * that can lead to both a L and a D , i.e.: c * ∈ A(a L ) ∩ A(a D ) (where A(a L ) and A(a D ) are the attraction sets of a L and a D ). There are several possible choices for c * . As we only want to identify the loop responsible for the coexistence of the two steady states, we choose c * to be the one with the highest hierarchical level. The operational graph G c * obtained is depicted in Figure 4.

Note that this graph contains only one feedback loop: C3a C8a, which is a classical double activation system, involving the caspases. From Thomas' conjectures, there can be no oscillations, only multistationarity. This is indeed the case, as there are exactly two steady states: C3a = C8a = 1 (apoptosis) and C3a = C8a = 0 (survival).

These two examples show that it is algorithmically possible, in the Boolean framework, to identify, in a complex interaction graph such as the one in Fig. 1 (comprising multiple positive and negative loops), subsets of loops which are effectively responsible for those two dynamical behaviors.

Each of the two operational graphs represents two different biological scenarios. The first operational graph (Fig. 3), represents the interactions that remain with respect to Boolean rules (Table 1).

active after a sufficiently long time interval of TNF stimulation has elapsed. That is, immediately upon TNF stimulation, the system responds and evolves towards a certain configuration; once this has been reached, most interactions have been "stabilized" or achieved a natural balance. At this stage only remain active the cycle corresponding to NFκB-activated transcription of IκB, and the subsequent inhibition of NFκB. In this case, oscillations in these two variables may be observed [START_REF] Hoffmann | The IkB-NFkB signaling module: temporal control and selective gene activation[END_REF], but all the other variables are constant. In particular, the model predicts that inhibitor of apoptosis proteins (IAP) is present at a low concentration, but the caspases at high concentration so, from the biological point of view, it may be expected that a very long stimulation with TNF will eventually lead to apoptosis. However, at this stage, the system may still "reverse" its apoptotic decision if TNF is shut down: the system will then leave the configuration shown in Fig. 3, and postpone the survival/death decision. In the absence of TNF, the system will evolve towards the configuration shown in Fig. 4, where only the positive feedback cycle representing the caspase cascase [START_REF] Eissing | Bistability analysis of a caspase activation model for receptorinduced apoptosis[END_REF] remains functional. From the biological point of view this means that, depending on the state of the system when TNF was shut down (or its initial state), the cell may decide between survival (attractor a L ) or initiating apoptosis (attractor a D ).

Probabilistic analysis of asynchronous dynamics

The asynchronous transition graph of a given network is a very general object, in the sense that it contains information on all the possible trajectories of the system: that is, for each state, the graph indicates every possible successor. However,

A c c e p t e d m a n u s c r i p t

the graph does not contain any indication on which of the successors is more probable at a given time. Probabilistic approaches to study trajectories of Boolean networks have been used before, under different modelling ideas. For instance, [START_REF] Shmulevich | Probabilistic boolean networks: a rule-based uncertainty model for gene regulatory networks[END_REF] introduces Probabilistic Boolean Networks (PBNs), where several families of logical rules are available (X = F r (X), r ∈ {1, . . . , R}). At each updating instant, the family R to be used is chosen according to a probability distribution, but once F r is chosen the updating scheme is synchronous. In [START_REF] Chaves | Robustness and fragility of boolean models for genetic regulatory networks[END_REF][START_REF] Fauré | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF][START_REF] Gonzalez | GINsim: a software suite for the qualitative modelling, simulation and analysis of regulatory networks[END_REF], asynchronous updating schemes for Boolean models (with a single family of logical rules) were developed, based on the nature of the biological process (for example, proteins are always updated at a higher frequency than mRNAs). In [START_REF] Fauré | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF] the variables are divided into priority classes, and each priority class is updated at its own frequency. The idea of asynchronous updates was coupled to PBNs in [START_REF] Faryabi | Optimal intervention in asynchronous genetic regulatory networks[END_REF].

Here, we will use the concept of priority classes, to construct a matrix of transition probabilities for the graph G, thus incorporating further quantitative biological knowledge into discrete dynamics' models. Application of our method to the apoptosis/NFκB network leads to estimation, among other quantities, of the probability of cell survival or apoptosis upon stimulation of death receptors. Such quantities can be compared with experimental data (see also numerical results in [START_REF] Chaves | Regulation of apoptosis via the NFκB pathway: modeling and analysis[END_REF] and references therein).

Construction of a probabilistic transition graph

In order to confront the general discrete asynchronous dynamics, given by the transition graph G (see Definition 3) with biological experiments, one has to be a little bit more specific on the updating strategies. Previous algorithms and results are essentially qualitative, as they are mainly attached to the structure, without considering updating strategies at all. If one wants to represent the different choices of updating strategies, in a more quantitative manner, one way is to introduce probabilities in the transition graph G, and to consider a choice of the updating order of the variables (that is, the choice of a particular trajectory in G) as the choice of a trajectory in a stochastic process.

G can be implemented by means of its adjacency matrix, that will be denoted A(G). Recall that each state of the system is a Boolean vector X ∈ Ω = {0, 1} n . Such a vector X is unequivocally associated with a unique integer s(X) ∈ {1, . . . , 2 n } (X is the binary decomposition of s(X)), so that in the following, "states" (i.e. elements of Ω) will be represented by integers lying in {1, . . . , 2 n }. The adjacency matrix A(G) is defined as follows:

A(G) = a i, j 1≤i, j≤2 n , with a i j = 1 if "state j" is a successor of "state i", a i j = 0 otherwise.

A c c e p t e d m a n u s c r i p t

The size of A(G) is 2 n ×2 n and, in general, A(G) is implemented as a sparse matrix. As a first example, we propose a naive construction of the transition probabilities p i j , where all asynchronous successors of a state have the same probability. The uniform distribution of probabilities of asynchronous successors means that we make no a priori assumption about the updating rules (beyond of course Assumptions 1 and 2). Define the matrix P(G) as the 2 n × 2 n real matrix:

P(G) = p i, j 1≤i, j≤2 n , with: ∀1 ≤ i, j ≤ 2 n , p i j = a i j 2 n k=1 a ik = a i j N + (i) ,
where N + (i) = 2 n k=1 a ik designates the number of directed edges of G that leave i (in other words, N + (i) is the number of asynchronous successors of state i), and p i j ∈ [0, 1] is the probability for the system to go from state i to state j. By construction, matrices P(G) and A(G) share the same sparsity pattern. Furthermore, it is easy to see that P(G) is a stochastic matrix, i.e. all its elements p i j are nonnegative, and the sum of the elements of each row, 2 n j=1 p i j is equal to 1. From the construction of P(G), we can now consider the (asynchronous) dynamics on G as a discrete time Markov chain, over the state space {1, . . . , 2 n }. In particular, if P(X 0 = i 0 ) designates the probability for the system to be initially in the state i 0 ∈ {1, . . . , 2 n }, then the probability that a trajectory follows the path p = (i 0 , i 1 , . . . , i q ) (where the i j are elements of {1, . . . , 2 n }) is equal to:

P(p) = P(X 0 = i 0 ) q-1 j=0 p i j ,i j+1 .
A similar matrix of probabilities can be constructed for the transition graph G scc , which describes transitions among strongly connected components (SCCs). As explained in Section 3.1, a SCC c ∈ {1, . . . , p} is, by definition, a subset of the state space Ω. In the following, S (c) denotes the subset of {1, . . . , 2 n } of the states belonging to c, and |S (c)| denotes the size of S (c). The matrix P scc is a square real matrix of size p × p, and its entries are given by:

p scc c,c = 1 |S (c)| i∈S (c) j∈S (c ) p i j . (1)
In words, this definition means that to compute the probability of the transition from a SCC c to a SCC c , we sum all transitions from any state i in S (c) to any state j in S (c ), and we divide this sum by the number of states in c, in order to obtain an average transition probability from c to c . Moreover, we prove the following:

A c c e p t e d m a n u s c r i p t

Proposition 2 The matrix P scc is stochastic.

Proof.

The nonnegativity of P scc elements is obvious. We prove here that the sum of the elements of its rows is equal to 1. Let c ∈ {1, . . . , p}. We have:

p c =1 p scc c,c = p c =1 1 |S (c)| i∈S (c) j∈S (c ) p i j = 1 |S (c)| p c =1 i∈S (c) j∈S (c ) a i j N + (i) = 1 |S (c)| i∈S (c) 1 N + (i) p c =1 j∈S (c ) a i j .
As the SCCs of a directed graph form a partition of its set of vertices, the quantity p c =1 j∈S (c ) a i j is equal to the number of edges that leave i (in graph G). By definition, it is equal to N + (i). This leads to:

p c =1 p scc c,c = 1 |S (c)| i∈S (c) N + (i) N + (i) = |S (c)| |S (c)| = 1.
The main advantage of considering matrix P scc instead of matrix P(G) is that it satisfies some useful properties. For instance, each ergodic class of P scc (roughly, an ergodic class is a set of non transient SCCs, see [START_REF] Berman | Nonnegative matrices in the mathematical sciences[END_REF] for a precise definition) contains only one element. The corresponding Markov chain is then called absorbing and its absorbing elements are in fact the attractors of the system. For absorbing chains, we can use the following well known result (see, e.g. [START_REF] Grinstead | Introduction to probability[END_REF]): there exists a p × p permutation matrix P 2 such that

P 2 P scc P t 2 = I r 0 R Q , (2)
where r ∈ {1, . . . , p} is the number of attractors of the system, I r denotes the r × r identity matrix, and Q is a (pr) × (pr) lower triangular matrix that satisfies:

lim n→∞ Q n = 0. (3)

A c c e p t e d m a n u s c r i p t

The form (2) is often called the canonical form of P scc . From (3), we can define the (pr) × (pr) matrix N = (I -Q) -1 (often called fundamental matrix). Its entry n cc gives the expected number of times that the process is in the transient SCC c if it started somewhere in the transient SCC c. In particular, it can be used to compute the vector t = N1 (where 1 designates the column vector of whose entries are 1). Given a transient SCC c, the entry t c of t gives the expected number of steps before the chain reaches an attractor, given that it started somewhere in c. As a consequence, Markov chains provide an efficient mathematical framework, in which useful global parameters can be computed. These parameters provide important biological knowledge about the system, as they contribute to further characterize qualitative dynamical properties, that are robust with respect to the topology of the network. (See next Sections 5.2 and 5.3.)

Towards a biological probabilistic graph

The matrix of transition probabilities, P(G), associated with the asynchronous graph in Section 5.1 was based on a uniform probability distribution. The probabilities of transition can instead be computed according to biological data, by using the notion of priority classes [START_REF] Fauré | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF][START_REF] Gonzalez | GINsim: a software suite for the qualitative modelling, simulation and analysis of regulatory networks[END_REF]. For many biological networks, partial knowledge of the parameters is often available. For example, the relative rates of two reactions are known (e.g., the rate of formation of protein A is larger than that of protein B). This biological knowledge can be straigthforwardly incorporated into the transition graph, by stipulating an updating strategy such as an updating order among all variables. Roughly, the idea is to group the variables into several groups, called priority classes, and assign a weight to each of these groups: higher weights denote a more probable transition. A similar idea was used in [START_REF] Chaves | Robustness and fragility of boolean models for genetic regulatory networks[END_REF], where two classes were considered, one for proteins and another for mRNAs. The updating order stipulated that proteins were always updated first and mRNAs next. More generally, to implement the notion of priority classes, consider ρ classes C 1 , . . ., C ρ ⊂ V and their respective weights,

W 1 > W 2 > • • • > W ρ .
Recall that, by assumption 1, only one variable is updated at a given instant, so each transition in G involves a 0/1 switch in exactly one variable. Consider an edge i → j in G, and suppose that v s ∈ V is the node that switches. Then associate with edge i → j the value:

w i j = W r , if v s ∈ C r

A c c e p t e d m a n u s c r i p t

that is, if the node v s updated in the transition from state i to state j belongs to class C r . If no transition from i to j is possible, then set w i j = 0. Then define a new transition matrix, where each p i j represents a weighted average:

P bio (G) = p i, j 1≤i, j≤2 n , with: ∀1 ≤ i, j ≤ 2 n , p i j = w i j a i j 2 n k=1 w ik . (4)
As before, a corresponding matrix, P scc bio , can be constructed for the graph G scc . The probability of transition between two SCCs c and c is given by Eq. ( 1), where the p i j are replaced by the transition probabilities computed in (4). Again, it is easy to check that P scc bio represents an absorbing Markov chain process.

Application to the apoptosis network

In order to illustrate the probabilistic approach, and the type of results it provides, an application to the apoptosis network is next described. In particular, the results show how Tumor Necrosis Factor (TNF) influences the choice of the system between the two possible steady states: the "survival" of the cell (attractor a L , with inhibition of the caspases) and the triggering of apoptosis (attractor a D , with activation of the caspases). Experimentally, it is observed that a cell irreversibly enters the apoptotic pathway once a certain threshold in caspase activation has been reached (see [START_REF] Rehm | Single-cell fluorescence resonance energy transfer analysis demonstrates that caspase activation during apoptosis is a rapid process[END_REF]). In turn, caspase activation is observed to depend on the duration of TNF stimulation, as well as on TNF concentration (typically, the caspase activation threshold is reached faster for higher TNF concentrations). In [START_REF] Chaves | Regulation of apoptosis via the NFκB pathway: modeling and analysis[END_REF], this property was statistically observed by computing many different trajectories (of a continuous, piecewise linear system), with different updating strategies. Within the framework presented in this paper, the probability that the cell follows the "survival" or "apoptosis" pathway can be computed directly from the matrices P scc (or P scc bio ), without performing large numbers of simulations. Recall that a trajectory of the system will converge towards a L or a D in the absence of TNF. Technically, the shutdown of death receptor stimulation is represented in our system by switching the input variable TNF from 1 to 0. For each state X in T 1 (where TNF is equal to 1), a successor X s is computed in T 0 (where the variable TNF is 0, and all other variables stay unchanged). Then, using the matrix P scc and its canonical form (2), the probabilities to reach attractors a L and a D from initial state X s are computed (Figure 5).

Contrary to [START_REF] Chaves | Regulation of apoptosis via the NFκB pathway: modeling and analysis[END_REF], where time is represented by a continuous variable, in the asynchronous graphs there is no direct measure of time. However, the hierarchical levels of the graph G scc do give an indication of time progression. Indeed, consider two states along any given trajectory, X 1 , X 2 ∈ T 1 where X 1 belongs to a hierarchical level H i and X 2 belongs to H j , with j > i. Then we can say that TNF stimulation has been longer for X 2 than for X 1 . The x-axis of Figure 5, which represents the hierarchical levels, is thus a relative measure of the duration of TNF. Following the procedure developed in Section 5.2, a more realistic matrix P scc bio can be constructed. Based on the parameters reported in [START_REF] Eissing | Bistability analysis of a caspase activation model for receptorinduced apoptosis[END_REF][START_REF] Schliemann | Modelling and experimental validation of TNFα induced pro-and antiapoptotic signalling[END_REF] (and references therein), and on the results from Section 4.3, three priority classes were established, based on biological function groups (Table 3). The first class represents the negative feedback loops related to the oscillatory behaviour: NFκB, NFκB nuc , IκB, and IKKa. The second class represents the caspases' positive loop and its activation: C3a, C8a, and T 2 . Finally, the third class contains the remaining proteins (A20a, FLIP, IAP, and CARP) whose transcription is regulated by NFκB nuc .

A c c e p t e d m a n u s c r i p t

Setting up these classes allowed us to compare the contribution of each family of operational interactions to the final survival/apoptosis decision. In a first study, we computed the probabilities of apoptosis and survival, for all possible combinations of w i ∈ {1, 5} (i = 1, 2, 3). This captures the effect of the relative updating frequencies (results are summarized in Table 4). Some examples are shown in Figure 5, which depicts the average probabilities of reaching attractors a L and a D , after 50 updating steps with TNF=0,

P L (H i ) = 1 N 1 i N 1 i k=1 P 50 bio (X 0 k , a L ), P D (H i ) = 1 N 1 i N 1 i k=1 P 50 bio (X 0 k , a D ),
where N 1 i denotes the number of states in T 1 ∩ H i , and X 0 k ∈ T 0 denotes the successor of X k ∈ T 1 after TNF is switched off (computed as explained above). After 50 steps, the system is likely to have reached an attractor with (practically) 100% probability. Thus, P L (H i ) ≈ 1 -P D (H i ). In all cases studied, we observe that the probability of survival is larger than that of apoptosis, for almost all levels H i , i = 1, . . . , 37. This situation is reversed for the last level, H 38 , where the probability of apoptosis is always larger than that of survival (see Table 4, column P L (H 38 )). It is striking that this property appears in all cases, suggesting the existence of a threshold state configuration which is independent of the updating strategy, and constitutes a robust feature of the system. In other words, if one wants to promote apoptosis, TNF should be sustained long enough for the system to reach attractor a 0 , i.e. apoptotic oscillations.

More generally, we observe that the probability of apoptosis (or survival) can be modulated within a certain range of values. At the last level, the probability of apoptosis is above 50%. In the levels H 1 to H 37 , the probability of reaching a L dominates the probability of reaching a D . Nevertheless, P D (H i ), for i = 1, . . . , 37, is maximal (between 30 and 40%) when w 1 > w 2 , that is, when NFκB-IκB pathway interactions are more frequent than caspase cascade interactions. (Figure 5d). In fact, it is interesting to note that, as the frequency of NFκB-IκB inter-A c c e p t e d m a n u s c r i p t

Conclusion

A method for model reduction of Boolean networks has been developed, based on the hierarchical decomposition of asynchronous graphs. The first aspect to be analyzed is the decomposition of the state space of the n-dimensional Boolean network into strongly connected components (SCCs), and the construction of the graph of transitions among them. The SCCs can be viewed as the "new states" of a "new" reduced system, since very often the number of SCC is less than or equal to the number of states of the original system. The second aspect in the model reduction method is the reconstruction of the Boolean rules that represent the graph of transitions among the SCCs. An identification algorithm (known as REVEAL) was adapted and used to determine a family of Boolean rules that describe the dynamics represented by a (sub-)graph of transitions.

More generally, the model reduction method uses the structure of interactions to isolate and identify smaller subsystems (or groups of variables and interactions) responsible for a given qualitative dynamical behaviour. This is a particularly relevant characterization for biological systems where experimental data consists (mostly) of qualitative measurements. The techniques described here are based on the fact that, for a Boolean system, the whole state space can, in principle, be enumerated; this introduces one other limitation to the method, on the size n of the network that can be computationally managed. Networks of intermediate size (up to n = 15) are easily computed. For larger networks, one may still use this method, by first isolating more basic modules. Each module would then be separately reduced and treated as one "node" with its Boolean rule.

As illustrated with the apoptosis example, model reduction using the asynchronous graph decomposition is a powerful potential source of valuable knowledge on a system. All the possible qualitative trajectories of the system are characterized, as well as their robustness to environmental perturbations. It is possible to identify the mechanism (in the form of smaller groups of variables and interactions) which is responsible for a given asymptotic behaviour of the system, for instance, the existence of oscillatory dynamics or (multi-)stability. Finally, the asynchronous transition graph can be naturally associated with a matrix of transition probabilities. Biological knowledge on the system's kinetics can thus be incorporated to obtain a more quantitative description of the system. These are also useful tools to test hypothesis and generate predictions concerning the structure of interconnections and the importance of each variable to the overall dynamics.

Figure 1 :

 1 Figure 1: Interaction graph of the simplified model of regulation of apoptosis via the NFκB pathway. As noted in Remark 1, some edges do not have a fixed sign (influence of TNF on IκB, IAP and CARP).

A c c e p t e d m a n u

  s c r i p t organization of the transition graph allows the formulation of simple algorithmic definitions of attraction and reachability sets. Algorithm 1 -Computation of attraction and reachability sets in G scc . Input: c ∈ {1, . . . , p} (a SCC of the asynchronous transition graph). Output: A(c), R(c): attraction and reachability sets of c. 1: A(c) := ATTR(c) 2: R(c) := REACH(c) where ATTR and REACH are two simple recursive functions: 1: ATTR(c): 2: A := {c} 3: P := predecessors(c) 4: if P ∅ then 5:

1 :

 1 REACH(c): 2: R := {c} 3: S := successors(c) 4: if S ∅ then 5: for all γ ∈ S do 6: R := R ∪ REACH(γ) 7:

Figure 2 :

 2 Figure 2: Scheme of the main elements of the hierarchical graph G scc for the apoptosis network.The vertical line separates the two connected components (generated by the input TNF), and the horizontal lines separates the different hierarchical levels. The SCCs are designated by their integer index (between 1 and p = 1472). The only SCCs that are represented here are the roots (SCCs that belong to the first hierarchical level) and the attractors (in bold characters).

A c c e p t e d m a n u s c r i p t Definition 10

 10 Consider an n-dimensional Boolean network with interaction graph G = (V, E). Let G = (Ω, E) be its asynchronous state transition graph. Let c ⊂ Ω and let R(c) denote the set of states reachable from any state in c. The subgraph of G induced by R(c) is designated by G(c).

4 :S

 4 := X /* will contain the synchronous successor of X */ 5:

1 :

 1 compute reachability set R(c) /* use Algorithm 1 */ 2: compute subgraph G scc (c) of G scc induced by R(c) 3: compute corresponding subgraph G(c) of G 4: compute synchronous truth table (In, Out) from G(c) /* use step 1 */ 5: compute operational graphs /* use Alg. 2 */

Figure 3 :

 3 Figure3: Operational graph of the NFκB pathway in attractor a 0 (TNF = 1). The isolated variables have a fixed value, that can be found in Table2.

Figure 4 :

 4 Figure4: Operational graph of the NFκB pathway just before the choice between attractors a L (survival) and a D (apoptosis) (i.e. in the region where TNF= 0). Interactions have been signed with respect to Boolean rules (Table1).

Figure 5 :

 5 Figure 5: The curves represent the system's response to TNF switch off, after the system has reached a certain hierarchical level H i (x-axis). The dashed lines with points represent the average probability for the system to reach attractor a L (the "survival" steady state), starting from a state in hierarchical level H i . The solid lines with '+' signs represent the probability to reach attractor a D (the "apoptosis" steady state), starting from H i .The weights used in each case were: (a) w 1 = w 2 = w 3 = 1; (b) w 3 = 1, w 1 = w 2 = 5; (c) w 1 = w 3 = 1, w 2 = 5; and (d) w 2 = w 3 = 1, w 1 = 5.

A c c e p t e d m a n u s c r i p t

Table 2 :

 2 Boolean patterns of the three attractors. Attractor a L is a steady state that corresponds to cell survival (or "life"). Attractor a D is another steady state that corresponds to the triggering of apoptosis (or "death"). Attractor a O is a SCC which contains 56 states. The symbol * means that the corresponding variable has no fixed value in the attractor and may oscillate between 0 and 1.

Table 3 :

 3 Priority classes.

generally, a given species is considered to be absent if its concentration is lower than a given threshold, present otherwise

see http://www.stanford.edu/ ˜dgleich/programs/matlab_bgl/

A c c e p t e d m a n u s c r i p t

where log designates the logarithm to the base 2, and P x M ∈ [0, 1] is the proportion of rows of M that are equal to the boolean vector x.

To justify the correctness of the algorithm, we prove the following Proposition 1 Let I = {i 1 , . . . , i k } ⊆ {1, . . . , n} and j ∈ {1, . . . , n}. The two assertions are equivalent:

(ii) There exists a function f j that is consistent with (In, Out) with respect to column j and such that σ( f j ) ⊆ {i 1 , . . . , i k }.

Proof.

Let M (resp. M ) designate the matrix formed by the columns C In (i 1 ), . . . , C In (i k ) (resp. by the columns C Out ( j), C In (i 1 ), . . . , C In (i k )). Definition 8 yields:

Moreover, we have:

As this is a sum of nonnegative terms, it is equivalent to:

which both imply that either P 0,x M or P 1,x M is equal to zero in each term of the sum, i.e. for each x ∈ {0, 1} k . In other words, this means that (i) is equivalent to the following assertion:

denotes the i-th row of In). Then there exists α x ∈ {0, 1} such that Out i j = α x for all i ∈ A x .

Suppose (i) (or equivalently (iii)) is satisfied, let us construct a function f j that satisfies (ii). To do that, take any Boolean vector x ∈ {0, 1} n . We have to consider two different cases. First, suppose there exists a row R i (In) of In such that In i,i l =

A c c e p t e d m a n u s c r i p t

Weights Survival probability (w 1 , w 2 , w 3 ) P L (H 38 ) actions becomes larger, also the probability of apoptosis increases (compare cases w 1 = 1, 3, 5, 10 with w 2 = w 3 = 1). In contrast, for the cases with w 1 = 1 (the NFκB-IκB interactions happen with smaller frequency), the probability of survival is significantly higher (larger than 80%) (Figure 5a-c). This suggests that the decision between cell survival and apoptosis strongly depends on the dynamics of the anti-apoptotic NFκB-IκB pathway: high "turnover rates" (that is, the interactions in this pathway happen at a faster rate) increase the probability of apoptosis. The contribution of this Markov chain approach to the analysis of biological networks is thus two-fold: first, the common traits of the curves obtained with different matrices P scc suggest global qualitative dynamical properties, which are robust with respect to the structure of the system, that is to say, independent of the choice of the updating order of the variables. Second, the quantitative aspects capture the variability and possible operating range of the network. Further applications of this technique include hypothesis testing and validation. For instance, one can easily analyze the effect of new interactions in the system's structure; similarly, the relative impact of two proteins can be studied by comparing the response of the system with different priority classes and updating strategies.