Franck Delaplace
email: franck.delaplace@ibisc.univ-evry.fr.

Hanna ⋆⋆

Klaudel

Tarek Melliti

Sylvain Sené

Analysis of modular organisation of interaction networks based on asymptotic dynamics ⋆

Keywords: modularity, interaction networks, discrete dynamics, equilibria

This paper investigates questions related to modularity in biological interaction networks. We develop a discrete theoretical framework based on the analysis of the asymptotic dynamics of biological interaction networks. More precisely, we exhibit formal conditions under which agents of interaction networks can be grouped into modules, forming a modular organisation. Our main result is that the conventional decomposition into strongly connected components fulfills the formal conditions of being a modular organisation. We also propose a modular and incremental algorithm for an efficient equilibria computation. Furthermore, we point out that our framework enables a finer analysis providing a decomposition in elementary modules, possibly smaller than strongly connected components.

Introduction

The analysis of the relations between the structure of a biological system and the related biological functions that identify specific states describing particular behaviours is among the most challenging problems [START_REF] Monod | Chance and necessity: an essay on the natural philosophy of modern biology[END_REF] at the frontier of theoretical computer science and biology. Let us introduce an illustration of these structure/function relations. On the one hand, gene regulation may be structured into a directed graph, called the interaction graph, from which a dynamics is computed. On the other hand, the attractors (i.e., stable configurations and/or sustained oscillations) of such a dynamics identify the functions of the system. For instance, for the bacteriophage λ, the reciprocal regulations between genes Cro and cI induce two biological functions, namely the lysis and the lysogeny [START_REF] Delbrück | Viruses[END_REF][START_REF] Lederberg | Genetic studies of lysogenicity in Escherichia coli[END_REF][START_REF] Jacob | Genetic regulatory mechanisms in the synthesis of proteins[END_REF], each corresponding to a distinct attractor.

Generally, studying complex biological functions relies on their decomposition into sub-functions identifying some basic behaviours. Each subfunction is supported by a part of the structure. In the context of gene regulation, this part corresponds to a sub-graph of the interaction graph. Thus, the whole system can be viewed in a modular way, where modularity establishes the link between the parts and their related sub-functions. The module composition refers to a structural composition as well as a dynamical one.

In the literature, methods related to module discovery in interaction networks are generally based on both the analysis of the network structures (a field close to graph theory) and the study of their associated dynamics [START_REF] Qi | Modularity and dynamics of cellular networks[END_REF]. Structural analysis identifies sub-networks with specific topological properties motivated either by a correspondence between topology and functionality [START_REF] Gagneur | Modular decomposition of protein-protein interaction networks[END_REF][START_REF] Chaouiya | A modular, qualitative modeling of regulatory networks using Petri nets[END_REF] or by the existence of statistical biases with respect to random networks [START_REF] Rives | Modular organization of cellular networks[END_REF]. Specific topologies like cliques [START_REF] Spirin | Protein complexes and functional modules in molecular networks[END_REF], or more generally strongly connected components (SCCs) are commonly used to reveal modules by structural analysis. Particular motifs [START_REF] Milo | Network motifs: simple building blocks of complex networks[END_REF][START_REF] Alon | Biological networks: The tinkerer as an engineer[END_REF] may also be interpreted as modules viewed as basic components. They represent over-represented biological sub-networks with respect to random ones. Moreover, dynamical analysis lays on the hypothesis that expression profiles provide insights on the relations between regulators, modules being possibly revealed from correlations between the expressions of biological agents. For instance, using yeast gene expression data, the authors of [START_REF] Bar-Joseph | Computational discovery of gene modules and regulatory networks[END_REF][START_REF] Segal | Module networks: discovering regulatory modules and their condition specific regulators from gene expression data[END_REF] inferred modules from co-regulated genes and the condition under which the regulation occurs. As a consequence, the discovery of a modular organisation in biological interaction networks is closely related to the influence of agents on one another and needs to investigate their expression dynamics [START_REF] Thieffry | The modularity of biological regulatory networks[END_REF][START_REF] Han | Understanding biological functions through molecular networks[END_REF].

In [START_REF] Siebert | Dynamical and structural modularity of discrete regulatory networks[END_REF][START_REF] Delaplace | Discrete causal model view of biological networks[END_REF][START_REF] Demongeot | Attraction basins as gauges of robustness against boundary conditions in biological complex systems[END_REF], the authors point out the need to relate structure and function to deal with modular organisations. The objective of this article is to define formally the notion of modular organisation as a list of modules together with a composition operation so that the dynamics of the module composition meets the global dynamics of the system. Indeed, modularity is somehow related to an invariance property of module asymptotic dynamics against regulatory perturbations of other modules [START_REF] Bernot | Behaviour preservation of a biological regulatory network when embedded into a larger network[END_REF], which supports the idea of viewing the global dynamics as the composition of the module's dynamics. Thus, using a discrete model of biological interaction networks [START_REF] Thomas | Boolean formalisation of genetic control circuits[END_REF][START_REF] Thomas | On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations[END_REF], we propose an approach that analyses the conditions of module formation and characterises the relations between the global behaviour of a network and the local behaviours of its components. We show under which conditions interaction networks can be divided into modules. As main results, we propose a modular and incremental algorithm to compute equilibria and we show that the conventional structural decomposition into strongly connected components fulfils the formal conditions of being a modular organisation. Furthermore, we show that our framework enables a finer analysis providing a decomposition in elementary modules, possibly smaller than strongly connected components. The paper is structured as follows: First, Section 2 introduces the main definitions and notations used throughout the paper. Section 3 presents the central notion of a modular organisation of a network along with its structural and dynamical properties. Section 4 defines elementary modular organisation and the conditions leading to obtain it. Some concluding remarks and perspectives are provided in Section 5.

The interaction network and its associated dynamics

This section introduces the discrete based asynchronous dynamics, modelling the dynamics of biological networks.

Relation. First, we introduce basic notations. Let ⇀ ⊆ S ×S be a binary relation on a set S, given s, s ′ ∈ S and S ′ ⊆ S, we denote by s ⇀ s ′ the fact that (s, s ′ ) ∈ ⇀, by (s ⇀) {s ′ | s ⇀ s ′ } the image of s by ⇀, and by (S ′ ⇀) its generalisation to the state set S ′ . Similarly, we denote by (⇀ s) and (⇀ S ′ ) the corresponding preimages. The composition of two binary relations will be denoted by ⇀ • ⇀ ′ and the reflexive and transitive closure by ⇀ * = i∈IN ⇀ i , with ⇀ 0 as the identity relation.

States and operations on states. Given a set A = {a 1 , . . . , a n } of agents of interest, each a i ∈ A has a local state, denoted by s a i , taking values in some nonempty finite set S a i . In the examples, all S a i are Boolean sets {0, 1}, but the proposed framework is not restricted to it. A state of A (or a configuration) is defined as a vector s ∈ S associating to each a i ∈ A a value in S a i , where S S a 1 × . . . × S an is the set of all possible states. For any X ⊆ A and s ∈ S, we denote by s| X the restriction of s to the agents in X, and by s| X the completion of s by all the values of agents in X; these notations extend to sets of states naturally. For example, the completion of the state s a 2 = 0 by the set of agents {a 1 , a 3 } is s a 2 | {a 1 ,a3 } = {000, 001, 100, 101}, and the restriction of s = 101, s ∈ S a 1 × S a 2 × S a 3 , on {a 2 } is 0. The X-equivalence defines an equivalence relation on states with regard to the state restriction on the agent set X. Definition 1. Two states s 1 , s 2 ∈ S are said to be X-equivalent and denoted by s 1 ∼ X s 2 , for some X ⊆ A, if and only if s 1 | X = s 2 | X , i.e., if they cannot be distinguished in S| X .

Evolution and asynchronous dynamics. An evolution is a relation on states ⇀. Each s ⇀ s ′ is a transition meaning that s evolves to s ′ by ⇀. Thus, the global evolution of η can be represented by a directed graph G = (S, ⇀) called the state graph. In this work, we pay particular attention to local evolutions, since each agent a ∈ A has its own evolution ⇀ a . The collection of all these local evolutions results in the asynchronous view of the global evolution of η, i.e., ⇀= a∈A ⇀ a . Definition 2. The asynchronous dynamics (or dynamics for short) of a network η is the triple A, S, (⇀ a ) a∈A , where A is a set of agents, S is a set of states, and for each a ∈ A, ⇀ a ⊆ S × S is a total or empty relation characterising the evolution of agent a such that for any s ⇀ a s ′ , either s = s ′ or s differs from s ′ only on the a-th component.

Interaction network and interaction graph. We are now in a position to introduce formally the interaction network as a family of functions η = {η a } a∈A , such that each η a : S → S a defines the next state η a (s) with respect to the asynchronous evolution of a from s. Network η allows to deduce a directed interaction graph G (A, -→) such that a i -→ a j if a i occurs in the definition of η a j . When ⇀ a is empty for some a (i.e., the local state of a remains invariant), then a plays the role of an input, which means that no other agents of A influence it (i.e., there are no arcs towards a in G), see Figure 1.

Orbit and equilibrium. Given a set S ′ ⊆ S, we introduce the following notions:

an orbit of S ′ , Ω(S ′ ), is the set of states comprising S ′ and all the states reachable from S ′ by ⇀; -an equilibrium e ∈ S is a state reachable infinitely often by ⇀; Ψ(S ′ ) denotes the set of equilibria reachable from S ′ ; -an attractor is a set of equilibria E ⊆ S such that ∀e ∈ E : Ψ({e}) = E.

In a state graph, an attractor is the set of states comprised in one terminal strongly connected component 3 that can be of two kinds:

η =        ηa 1 (s) = sa 1 ηa 2 (s) = sa 1 ∨ sa 3 ηa 3 (s) =
a 2 a 2 a 2 a 2 a 2 a 2 a 2 a 2 a 3 a 3 a 3 a 3 a 3 a 3 a 3 a 3 a 4 a 4 a 4 a 4 a 4 a 4 a 4 a 4
Fig. 1. An interaction network η (top left), its graphical representation (top right), and the state graph G of η composed of two disconnected components (bottom). G represents the dynamics of η for each state s ∈ {0, 1} 4 , in which by convention, self loops are omitted, stable states are depicted in gray while limit sets are in black.

• a stable state is a singleton E ⊆ S;

• a limit set is an attractor E such that |E| > 1.
Moreover, the restriction of ⇀ to X ⊆ A is defined as: ⇀ X = a∈X ⇀ a . Orbits and equilibria are determined by two operators having two arguments, an agent set and a state set. Definition 3. The orbit operator Ω and the equilibrium operator Ψ , are defined as follows for X ⊆ A and S ′ ⊆ S:

-Ω X (S ′ ) = (S ′ ⇀ * X ); -Ψ X (S ′ ) = {s ∈ Ω X (S ′ ) | ∀s ′ ∈ S : s ⇀ * X s ′ =⇒ s ′ ⇀ * X s}.
The equilibrium operator Ψ X is idempotent, upper-continuous and monotone (see Proposition 2 in Appendix). The example in Figure 1 illustrates the dynamics of an interaction network. It is defined by η, each η a being the local transition function of agent a. Given a state s, the evolution s ⇀ a s ′ means that s ′ is obtained by applying

η a to s, i.e., s ⇀ a s ′ s ′ a = η a (s) ∧ ∀a ′ ∈ A \ {a} : s ′ a ′ = s a ′ . Let us remark that: -the orbit of {1111} is Ω({1111}) = {1111, 1101, 1100};
-1100 is an equilibrium, as well as 0101 and 0010 are; -the set of equilibria reachable from 0000 is {0xyz | x, y, z ∈ {0, 1}}; -two attractors exist: a stable state {1100} and a limit set Ψ({0000}).

Regulation. The regulation is a sub-relation of the interaction specifying a dynamics-based dependence between two agents. Agent a k regulates

η = ηa 1 (s) = sa 1 ∧ sa 2 ηa 2 (s) = sa 1 ∨ 1 a1 a2 00 01 10 11 a 1 a 2 a 2 a1 a2
Fig. 2. An interaction network, its state graph and the corresponding regulation graph.

agent a ℓ , if at least one modification of a state of a ℓ requires a modification of a state of a k .

Definition 4. an interaction, a k -→ a ℓ is a regulation if and only if there exist two states s, s ′ ∈ S such that (s

∼ A\{a k } s ′ ) ∧ ((s ⇀ a ℓ ) ≁ a ℓ (s ′ ⇀ a ℓ )). By extension, given X i , X j ⊆ A, X i -→ X j if and only if ∃a k ∈ X i , ∃a ℓ ∈ X j : a k -→ a ℓ .
It may arise that the graph of interaction differs from the graph of regulation because the interaction depends on the syntactic definition of a network whereas the regulation relies on a property of the dynamics. In Figure 1, the sets of regulators of agents a 1 , a 2 , a 3 and a 4 are respectively {a 1 }, {a 1 , a 3 }, {a 2 } and {a 3 }. Notice also that there are the following relations on sets of agents: {a 1 , a 2 } -→ {a 3 , a 4 }, {a 3 } -→ {a 2 , a 4 } and {a 1 } -→ {a 2 , a 3 } -→ {a 4 }. Another example of a regulation graph is given in Figure 2 (right). It shows that interaction (a 1 , a 2 ) in the interaction network is actually not a regulation because no modification of a 1 influences the state of a 2 . All other interactions are effective, meaning that the underlying regulation graph contains all interactions from the network but (a 1 , a 2 ).

Composition of equilibria

In this section, a relation between the equilibria of an interaction network and that of its parts is presented. It allows to consider a modular view of the system in which each part is seen as a module, i.e., a subset of agents. It means that modules, which influence each other, reveal the underlying biological functions materialised by their equilibria.

Modular organisation

Our objective is to find a decomposition of the set of agents A into modules, i.e., a partition4 of A, together with a composition operator ⊘ of the Fig. 3. An interaction network and its state graph.

η =    ηa 1 (s) = ¬sa 1 ∨ sa 2 ηa 2 (s) = 1 ηa 3 (s) = sa 2 ∨
equilibria of these modules, allowing to retrieve the global equilibria of the network. Finding an adequate operator ⊘ is a challenging question. Basically, a modular organisation (X 1 , . . . , X m ) should satisfy the following equation characterizing the composition of the module equilibria:

Ψ X 1 ⊘ . . . ⊘ Ψ Xm = Ψ m i=1 X i . (1)
One can easily see that, in general, taking ⊘ = ∪ for example is not a solution. If we consider the interaction network in Figure 3 and a partition into two parts {a 1 , a 2 } and {a 3 }, then the corresponding sets of equilibria are respectively {110, 111} and {000, 001, 100, 101, 011, 111}, while the set of global equilibria is {111}, which is not the union of the previous ones. However, one may see that, for the same parts, the computation of the equilibria of {a 3 } from the equilibria of {a 1 , a 2 }, gives the expected property

Ψ {a 3 } • Ψ {a 1 ,a 2 } = Ψ {a 1 ,a 2 ,a 3 } , whereas Ψ {a 1 ,a 2 } • Ψ {a 3 } = Ψ {a 1 ,a 2 ,a 3 } .
This suggests that the order in which parts are taken into account plays an important role in the definition of the composition operator. Unfortunately, in general none of the usual operators such as ∪, ∩ and • can be used as the modular composition operator as one can check in the following network: {η a 1 (s) = s a 1 ∧ ¬s a 3 , η a 2 (s) = s a 2 ∧ ¬s a 3 , η a 3 (s) = s a 2 }, while a modular decomposition exists: {a 1 } followed by {a 2 , a 3 }. Thus, we will focus on an ordered partition π = (X 1 , . . . , X m ) of A, i.e., a partition of A provided with a strict total order and represented by a sequence, called a modular organisation, preserving [START_REF] Monod | Chance and necessity: an essay on the natural philosophy of modern biology[END_REF]. Furthermore, we would like to be able to "fold" contiguous modules in π in order to deal with them as with a single module 5 , while preserving the result of the composition of equilibria. As a consequence, we require a modular organisation to support folding and to be such that the composition operator ⊘ is associative according to the order in π. In order to form a modular organisation, the modules and their order in π should satisfy some conditions related to their dynamics. Intuitively, two disjoint sets of agents X i and X j , i < j, can be modules in π, either if they do not regulate each other, or if X i regulates X j . In both cases, we can remark that the equilibria of X i should embed the asymptotic evolution of X j , which leads to encompass the equilibria of X j in the equilibria of X i . These conditions are expressed by the modularity relation (M -relation).

Definition 5. The M -relation ⊆ P(A) × P(A) is defined as:

X i X j ∀S ′ ⊆ S : (Ψ X i • Ψ X i ∪X j (S ′ ))⇀ X j ⊆ (Ψ X i • Ψ X i ∪X j (S ′ )).
Some fundamental properties of the M -relation can be found in Proposition 3 of the Appendix. In this context, a modular organisation can be defined as follows.

Definition 6. A modular organisation (X 1 , . . . , X m ) is an ordered partition of A such that for all 1 < i ≤ m:

( i-1 j=1 X j ) X i .
From Definition 6, Proposition 1 states that being a modular organisation is preserved by any folding of its contiguous parts 6 .

Proposition 1. Let π = (X 1 , . . . , X m ) be a modular organisation. For all

1 ≤ i ≤ j ≤ m, (X 1 , . . . , X i-1 , j k=i X k , X j+1 , . . . , X m ) is a modular organisation.
In the literature [START_REF] Gagneur | Modular decomposition of protein-protein interaction networks[END_REF][START_REF] Milo | Network motifs: simple building blocks of complex networks[END_REF][START_REF] Alon | Biological networks: The tinkerer as an engineer[END_REF], modules are frequently assimilated to SCCs of interaction networks. Although these works focus on structural arguments only, it turns out that they are compatible with Definition 6. Indeed, any topological order 7 of SCCs is actually a modular organisation. Notice that, a topological order on the quotient graph of SCCs always exist since the graph is acyclic. For instance, ({a 1 }, {a 2 , a 3 }, {a 4 }) is a modular organisation of the interaction network presented in Figure 1. In what follows, we present an approach addressing formally this aspect. As a result, we show that, in particular, the structural decomposition in SCCs makes sense and may be improved by a deeper analysis leading to the decomposition of SCCs in elementary modules (see Section 4), potentially smaller than those coming from SCCs.

Regulation and modularity relation

The regulation and the M -relation are related, as shown below.

Lemma 1. For any X i , X j subsets of A:

¬(X j -→ X i ) =⇒ X i X j . Xi Xi Xi Xi Xi Xi Xi Xi Xi Xi Xi Xi Xi Xi S ′ ∪ = ΨX i (S ′ ) Xj Xi Xi Xi Xi Xi Xi Xi Xi Xi Xi Xi Xi Xi Xi S ′ Xj Xj Xj [ ΨX j ]⇀ ↽X i (S ′ ) (a) (b)
Fig. 4. Successive steps leading to the definition of the composition operator ⊘. The SCCs defined for ⇀X i are in gray (a). The two terminal SCCs at the bottom correspond to attractors of ⇀X i (a). The equilibria ΨX i ∪X j (S ′ ) is computed from ⇀X j (bold Xj arrows) on these attractors (b).

According to Definition 6, Theorem 1 provides a connection between structural properties of a regulation graph and the corresponding modular organisations (possibly reduced to a single module).

Theorem 1. Any topological order of the SCC quotient graph of a regulation graph is a modular organisation.

Composition operator

In this section, we present the successive steps leading to the definition of the composition operator ⊘. From (1), ⊘ is a binary operator that applies on the equilibria of parts X i and X j of π, with i < j. Thus, its definition is based on the attractors of ⇀ X i which correspond to terminal nodes (terminal SCCs) of the SCC quotient graph of ⇀ X i , namely G/ ⇋ X i , where ⇋ X i is the equivalence relation identifying states belonging to the same SCC and defined as s 4.a). Moreover, for all S ′ ⊆ S, we denote by:

⇋ X i s ′ (s ⇀ * X i s ′ ) ∧ (s ′ ⇀ * X i s). For any S ′ ⊆ S, an attractor of ⇀ X i coincides with [s] ⇋ X i ⊆ Ψ X i (S ′ ) (see Figure
-[S ′ ] ⇋ X i = {[s] ⇋ X i |s ∈ S ′ } the set of equivalence classes of ⇋ X i in S ′ ; -[s] ⇋ X i [⇀ X j ] ⇋ X i [s ′ ] ⇋ X i ∃s ∈ [s] ⇋ X i , ∃s ′ ∈ [s ′ ] ⇋ X i : s ⇀ X j s ′
the evolution by agents of X j on these equivalence classes. We define an operator [ Ψ X j ] ⇋ X i , similar to the equilibria operator, computing the set of equilibria of [⇀ X j ] ⇋ X i in G/ ⇋ X i (see Figure 4.b and Section 3.4 for the algorithm) as follows:

[ Ψ X j ] ⇋ X i (S ′ ) {[s] ⇋ X i ∈ [S ′ ] ⇋ X i | (([s] ⇋ X i [⇀ X j ] * ⇋ X i ) ⊆ [S ′ ] ⇋ X i ) ∧ ∀[s ′ ] ⇋ X i ∈ [S] ⇋ X i : [s] ⇋ X i [⇀ X j ] * ⇋ X i [s ′ ] ⇋ X i =⇒ [s ′ ] ⇋ X i [⇀ X j ] * ⇋ X i [s] ⇋ X i }.
Operator ⊘ is thus defined as:

Ψ X i ⊘ Ψ X j Flat • [ Ψ X j ] ⇋ X i • Ψ X i , (2)
where, for any set E ⊆ P(S), Flat(E) = e∈E e flattens the set. Thus, if applied to a set of attractors, Flat gives the underlying set of equilibria. For example, Flat({{00, 01}, {10}}) = {00, 01, 10}. As a result, one can see that ⊘ does compute the set of states belonging to the attractors of X i which are also the equilibria of ⇀ X j . Lemma 2 below shows that the global equilibria of X i ∪ X j are obtained using the composition Ψ X i ⊘ Ψ X j .

Lemma 2. For all X i , X j disjoint subsets of A:

X i X j =⇒ ∀S ′ ⊆ S : Ψ X i ∪X j (S ′ ) = (Ψ X i ⊘ Ψ X j ) • Ω X i ∪X j (S ′ ).
As a main result, the computation of global equilibria of an interaction network can be obtained modularly (Theorem 2).

Theorem 2. Let A ′ = m i=1 X i ⊆ A be a set of agents, if (X 1 , . . . , X m
) is a modular organisation then we have:

Ψ A ′ = (Ψ X 1 ⊘ . . . ⊘ Ψ Xm ) • Ω A ′ .

Modular and incremental computation of equilibria

Modularity allows an incremental and efficient computation of the equilibria avoiding the generation of the complete state space S. The algorithm presented in Figure 5 is an application of Theorem 2 introducing incremental processing based the fact that the evolution of an agent only depends on the current states of its regulators. Indeed, the equilibria of a set of agents X in a state space S| X∪Y where X and Y are disjoint sets and X ∪ Y contains the regulators R X = (-→ X) of X, can be computed from the restriction to X ∪ R X completed by Y \ R X . In other words, the equilibria computation operator and the completion operator commute for an appropriate selection of sets of agents, as shown by the following lemma.

Lemma 3. Ψ X (S| X∪Y ) = Ψ X S| X∪R X Y \R X , with R X = (-→ X), R X ⊆ X ∪ Y and X ∩ Y = ∅.
Each step of the algorithm is seen as the computation of equilibria for the following modular organisation (X, X i ) where X = X 1 ∪ . . . ∪ X i-1 is the folding of the modules preceding X i in the initial modular organisation and corresponds to the following equation:

Ψ X∪X i (S) = Flat Ψ X i ⇋ X Ψ X (S| X∪R X ) | X i \R X . ( 3
)
The algorithm of Figure 5 is divided in two parts: part i corresponding to the completion of states by X i , and part ii corresponding to the computation of the attractors taking into account ⇀ X i . In part i the attractors are duplicated by completing their state values while preserving the structure of attractors. part ii is divided into two subparts. part ii.1 computes the quotient graph with attractors as vertices and the quotiented evolution by X i as arcs. Notice that some attractors are removed during this step. Indeed, from some states belonging to attractors, ⇀ X i may reach states located outside the attractors set. By definition of the M -relation (Definition 5), they cannot be considered as equilibria and are not included. The remaining set of attractors, called the core, is used in part ii.2 to compute the equilibria as the terminal SCCs of the core graph. Since TermSCCs returns a set of "attractors of attractors", they are finally flattened to retrieve the structure of a set of equilibria.

The complexity of the algorithm in the product of the maximal number of equilibria by the number of agents is exponential in general. It is however linear for acyclic regulatory graphs. Let α i be the number of equilibria computed at step i, and N the set of newly introduced agents. The computation time is bounded by k

• m i=1 α i • 2 |N | , for k ∈ IN.
The exponential time corresponds to the computation of the completion and of the quotient graph. |N | is bounded by β, the number of agents in the greatest SCC, under the assumption that all modular organisations are subdivisions of topological orders of the SCC quotient graphs. Indeed, in the modular organisation the regulators of X i always precede X i unless they are also regulated by X i . Hence, the computation time is bounded by k • 2 β • m • α where α stands for the maximal number of equilibria for all steps, leading to a complexity in O(2 β • m • α). In the worst case, notably corresponding to a regulatory graph reduced to a single SCC, the complexity is the same as the brute-force algorithm computing equilibria from the whole state graph (i.e., π = {A}, β = |A|). The algorithm is more efficient in practice. In particular, for networks whose interaction graph is acyclic, each module corresponds to a single agent (β = 1) leading to a complexity in O(|A| • α) with α ≤ 2 ∆ 0 where ∆ 0 is the number of all input agents having no regulators but possibly themselves. Hence, for regulatory path-graphs, the algorithm is linear in the number of agents because ∆ 0 = 1.

Elementary modular organisation

Informally, a module is elementary if it is not separable, i.e., if the equilibria of each of its agents depend entirely on the equilibria of all the others. For instance, consider negative circuits that lead to asymptotic sustained oscillations [START_REF] Thomas | On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations[END_REF]. In such regulation patterns, the equilibria of an agent cannot be encompassed into that of the others because, in order to reach its own equilibria, each agent evolves from the equilibria of all the others. In this context, a modular organisation provided by some topological orders of the SCC quotient graph (see Theorem 1) does not always provide an elementary decomposition. Figure 6 depicts an interaction network η composed of three agents a 1 , a 2 and a 3 with the associated strongly connected regulation graph. Its underlying state graph shows that the global dynamics of η leads to two attractors, stable state {111} and limit set {000, 100, 101, 001}. It is easy to see that {a 2 } {a 1 , a 3 } and that this M -relation is (obviously) preserved by folding, because there are only two modules. Hence, ordered partition ({a 2 }, {a 1 , a 3 }) is a modular organisation of η. However, separa- bility is not possible in general as illustrated in Figure 7. Indeed, starting from modular organisation π = ({a 1 }, {a 2 , a 3 }) obtained from the SCCs, the separation of {a 2 , a 3 } should lead to one of the ordered partitions π ′ = ({a 1 }, {a 2 }, {a 3 }) and π ′′ = ({a 1 }, {a 3 }, {a 2 }). The condition for π ′ to be a modular organisation is that {a 1 , a 2 } {a 3 }, i.e., the evolution by a 3 from the equilibria of {a 1 , a 2 } has to be included in the equilibria of {a 1 , a 2 }. We can observe that the attractors for {a 1 , a 2 } are {001, 101} and {010, 110}, while the evolution by a 3 from either 101 or 110 leaves the attractors of {a 1 , a 2 }, which means that π ′ is not a modular organisation. As a consequence, {a 2 , a 3 } cannot be separated. Indeed, here, although agents a 2 and a 3 are together M -related and thus can be separated a priori, they cannot be in the context of agent a 1 . The same reasoning applies for π ′′ . Hence, the separation condition of a module X i in π is not local to this module but depends on the module "context", that is the global equilibria (i.e., Ψ i-1 k=1 X k ) of the modules that precede. Deciding the separability of X i into X 1 i and X 2 i implies checking two conditions:

η =    ηa 1 (s) = (sa 1 ∧ sa 2 ) ∨ ¬sa 3 ηa 2 (s) = sa 1 ∧ sa 2 ∧ sa 3 ηa 3 (s) = sa 1
i-1 k=1 X k X 1 i and i-1 k=1 X k ∪ X 1 i X 2 i .
Of course, the complexity of the underlying computation is exponential in the size of π and also depends on the position of X i in π. Nevertheless, brute-force computation may be used in practice for small interaction networks (of about 15 agents). A more efficient method allowing to go beyond this limitation is, for the moment, an open question.

η =    ηa 1 (s) = ¬sa 1 ηa 2 (s) = (sa 1 ∧ sa 2 ) ∨ (sa 1 ∧ ¬sa 3 ) ∨ (sa 2 ∧ ¬sa 3 ) ηa 3 (s) = (sa 1 ∧ sa 2 ) ∨ (¬sa 1 ∧ sa 3 ) ∨ (sa 2 ∧ sa 3 ) a1 a2 a3 000

Conclusion

We developed a formal framework for the analysis of the modularity in interaction networks assuming asymptotic dynamics of modules and enabling their composition. We exhibited modularity conditions governing the composition of modules and an efficient computation method such that the global equilibria of interaction networks are obtained from the local ones, leading to an efficient algorithm. Moreover, we confirmed that usual assumptions identifying modules with SCCs have a strong motivation coming from theory. The next step should be identifying a characteristic property for finding elementary modular organisations. Then, since this work provides a rigorous setting for studying other questions around modularity. For example, a success factor of synthetic biology is to ensure the safety of modular design [START_REF] Purnick | The second wave of synthetic biology: from modules to systems[END_REF] which in our theoretical framework is guaranteed by construction through the concept of modular organisation. Also, questions related to robustness and evolution could be tackled thanks to the modular knowledge of interaction networks.

Proposition 3 (Properties of the modularity relation). Let S ′ be a subset of S. For all X i , X j subsets of A, we have the following properties:

1. X i X j ⇐⇒ Ψ X i ∪X j (S ′ ) ⊆ (Ψ X i • Ω X i ∪X j (S ′ )); 2. X i X j ⇐⇒ Ψ X i ∪X j (S ′ ) = Ψ X i • Ψ X i ∪X j (S ′ ). Proof. 1. (⇒) Let S ′ ⊆ S, s ∈ Ψ X i ∪X j (S ′ ) and s / ∈ Ψ X i (Ω X i ∪X j (S ′ )), and s ′ ∈ Ψ X i (Ω X i ∪X j ({s})). By definition of equilibrium, s ′ ⇀ * X i ∪X j s. Now, we have: -∀s ′′ ∈ (s ′ ⇀ X i ) : s ′′ ∈ Ψ X i • Ψ X i ∪X j ({s ′ })
, by definition of equilibria; -∀s ′′ ∈ (s ′ ⇀ X j : s ′′ ∈ Ψ X i • Ψ X i ∪X j ({s ′ }), by definition of . As a consequence, s ′ ⇀ * X i ∪X j s, which leads to a contradiction. (⇐) Let S ′ ⊆ S, s ∈ Ψ X i ∪X j (S ′ ), and Ψ X i ∪X j (S ′ ) ⊆ Ψ X i (Ω X i ∪X j (S ′ )). Then, by hypothesis, we have:

∀s ′′ ∈ (s ⇀ X j ) : s ′′ ∈ Ψ X i ∪X j (S ′ ) ∧ s ′′ ∈ Ψ X i (Ω X i ∪X j (S ′ )), which means that s, s ′′ ∈ Ψ X i • Ψ X i ∪X j (S ′ ). Thus, X i X j .

2. From Proposition 3.1 and since Ψ X i ∪X j (S ′ ) = Ω X i ∪X j (Ψ X i ∪X j (S ′ )).

Proof (Proposition 1). Let π = (X 1 , . . . , X i-1 , X i , X i+1 , . . . , X m ) be a modular organisation and let X = i-1 k=1 X k . We want to show that (X 1 , . . . , X i-1 , X i ∪ X i+1 , . . . , X m ) is a modular organisation. By definition 6, we have:

(X ∪ X i ) X i+1 (4)
and:

X X i . ( 5
)
We want to show that: (4) ∧ (5) =⇒ X (X i ∪ X i+1 ). First, by Proposition 3.2, we can write:

Ψ X∪X i • Ψ X∪X i ∪X i+1 = Ψ X∪X i ∪X i+1 by (4), (6)
Ψ X • Ψ X∪X i = Ψ X∪X i by (5). (7)
Thus:

Ψ X • Ψ X∪X i ∪X i+1 = Ψ X • (Ψ X∪X i • Ψ X∪X i ∪X i+1 ) by (6) = (Ψ X • Ψ X∪X i ) • Ψ X∪X i ∪X i+1 = Ψ X∪X i • Ψ X∪X i ∪X i+1 by (7)
= Ψ X∪X i ∪X i+1 by [START_REF] Gagneur | Modular decomposition of protein-protein interaction networks[END_REF], which is the expected result. From Proposition 3.2, we can deduce that: X (X i ∪ X i+1 ). Iteratively, we show that X j k=i X k . As a result, (X 1 , . . . , X i-1 , j k=i X k , X j+1 , . . . , X m ) is a modular organisation.

the following property holds by definition of the evolution: ∀s 1 , s 2 ∈ S| X∪Y : (s

1 ⇀ X s 2 ) ⇐⇒ (s 1 | X∪R X ⇀ X s 2 | X∪R X ∧ s 1 | Y \R X = s 2 | Y \R X ).
This property extends to the transitive closure by induction:

∀s 1 , s 2 ∈ S| X∪Y : (s 1 ⇀ * X s 2 ) ⇐⇒ (s 1 | X∪R X ⇀ * X s 2 | X∪R X ∧ s 1 | Y \R X = s 2 | Y \R X ). ( 8
)
First we prove that:

∀s ∈ S| X∪Y : s ∈ Ψ X (S| X∪Y ) ⇐⇒ s| X∪R X ∈ Ψ X (S| X∪R X ).

Let s ∈ Ψ X (S| X∪Y ). By definition of the equilibrium operator (Definition 3), s complies with the following equivalent property: ⇐⇒ ∀s ′ ∈ S| X∪Y : s ⇀ * X s ′ =⇒ s ′ ⇀ * X s. By application of Equation 8, we derive:

⇐⇒ ∀s ′ ∈ S| X∪Y : (s| Y \R X = s ′ | Y \R X ) =⇒ (s| X∪R X ⇀ * X s ′ | X∪R X =⇒ s ′ | X∪R X ⇀ * X s| X∪R X ) . Since s ⇀ *
X s ′ insures that s| X∪R X = s ′ | X∪R X , by Equation 8, we simplify and obtain:

⇐⇒ ∀s ′ ∈ S| X∪Y : s| X∪R X ⇀ * X s ′ | X∪R X =⇒ s ′ | X∪R X ⇀ * X s| X∪R X .
By definition of the equilibrium operator, this equivalently leads to: ⇐⇒ s| X∪R X ∈ Ψ X (S| X∪R X ). Now, by definition of the completion, any state s complies with the following property: s ∈ (s| X∪R X )| Y \R X . Since s| X∪R X ∈ Ψ X (S| X∪R X ), we deduce that s ∈ Ψ X (S| X∪R X )| Y \R X . Hence, we conclude that:

∀s ∈ S| X∪Y : s ∈ Ψ X (S| X∪Y ) ⇐⇒ s ∈ Ψ X (S| X∪R X )| Y \R X .

¬sa 2 ηa 4 (

 4

2 /Fig. 5 .

 25 Fig. 5. Algorithm of modular and incremental computation of equilibria

3 Fig. 6 .

 36 Fig.6. A separable network with modular organisation ({a2}, {a1, a3}).

Recall that, in a terminal strongly connected component each path starting from a vertex of the component remains in this component.

A partition of a set A is a set of nonempty disjoint subsets of A which covers A.

The folding of modules corresponds to the union of these modules.

Proofs are in Appendix.

A topological order is a total order obtained by topological sorting.

⋆ This work is supported by the project synbiotic of French National Agency for Research, anr blan-0307-01.

Appendix

Proposition 2 (Properties of Ψ ). Let A, S, (⇀ a ) a∈A be an asynchronous dynamics and X ⊆ A be a subset of agents. Ψ X has the following properties, for all sets S ′ , S ′′ subsets of states: a. Idempotency: Ψ X • Ψ X (S ′ ) = Ψ X (Ψ X (S ′ )) = Ψ X (S ′ ); b. Upper-continuity: Ψ X (S ′ ∪ S ′′ ) = Ψ X (S ′ ) ∪ Ψ X (S ′′ ); c. Monotony (order-preserving) : S ′ ⊆ S ′′ =⇒ Ψ X (S ′ ) ⊆ Ψ X (S ′′ ).

Proof (Proposition 2). Let Eq X (s) be the predicate meaning that s is an equilibrium for ⇀ X . a. By expanding Ψ X (Ψ X (S ′ )), we have:

X is upper-continuous on the lattice of state sets, we have:

c. An upper-continuous function is monotone.

Proof (Lemma 1). By Definition 4, for any X i , X j ⊆ A and for any s, s ′ ∈ S, we have ¬(

This property is obviously preserved at equilibria. Indeed, for any s, s ′ ∈ S, we have ¬(

Thus, the restrictions Ψ X i (s) and Ψ X i (s ′ ) to X i are identical. Then, the evolution by X j from the equilibria of X i remains in the equilibria of X i . Hence, we get that ¬(X j -→ X i ) =⇒ X i X j .

Proof (Theorem 1). Observe that, in the SCC quotient graph G of a regulation graph, X i -→ X j always implies that ¬(X j -→ X i ), because of the acyclicity of G. Thus, folding contiguous modules with respect to any topological order preserves the absence of regulation. As a consequence, if (X 1 , . . . , X m ) is a topological order of G, for all i, j ∈ IN such that 1 ≤ i ≤ j ≤ m, we have ¬(X j -→ X i ), and by Lemma 1, X i X j .

Similarly, an evolution by X i ∪ X j from an attractor of X i remain in the same attractor except potentially with evolutions by X j . We have then, for all s, s ′ ∈ Ψ Xi∪X j (S ′ ):

Now, since both s and s ′ belong to attractors of X i ∪ X j (by hypothesis), if there exists an evolution by X j from [s] ⇋ X i to [s ′ ] ⇋ X i , there exists obviously another path labelled by X j from [s ′ ] ⇋ X i to [s] ⇋ X i . Hence, for all s, s ′ ∈ Ψ X i ∪X j (S ′ ), we have:

As a result, we have 2) and Proposition 3.2, we can write:

Hence, for all s ∈ Ψ X∪Y (S ′ ), we have: s ∈ Ψ X ⊘ Ψ Y • Ω X∪Y (S ′ ), which corresponds to the following inclusion:

Proof (Theorem 2). This proof is made directly by induction on the modular organisation, using Definition 6 and Lemma 2. Since π = (X 1 , . . . , X m ) is a modular organisation, it is folding preserving and m-1 i=1 X i X m . Then, using (2) and Lemma 2, we have:

As a result, we obtain: Ψ A ′ (S ′ ) = (Ψ X 1 ⊘ . . . ⊘ Ψ Xm ) • Ω A ′ (S ′ ), which is the expected result.

Proof (Lemma 3). The evolution is governed by the values of the regulators and the evolution concerns the states of agents of X only. Under the assumptions that X and Y are two disjoint sets and that R X ⊆ X ∪ Y ,


Claire Chauvin

PIANO Reference Manual Python Interface for Assimilation in NemO

1 First steps with PIANO PIANO intends to provide a user-friendly interface to a full run of assimilation. For the moment, PIANO is only about one full cycle of variationnal assimilation. Its aims is to manage runs and communications between the model NEMO and the Incremental minimization of the cost function available in NEMOVAR. In the following, we make a distinction between the namelist -the file that contains all the definition of all the parameters needed by a run -, and a namelist blocka section of namelist defined by its name &nam_block. Recall that a NEMO run needs a namelist called namelist in the outer directory, and NEMOVAR asks for the namelist namelist.nemovar in the inner directory. These two namelists contain an exhaustive list of all the existing namelist blocks, needed by a run of NEMO (or NEMOTAM / NEMOVAR). This example contains an exhaustive list of all the arguments defined in file descr. Some are not implemented yet (parallel architecture not managed). The following parameters of final namelist are computed from ndays : nitend, nstock, nwrite, nstockfl, nwritfl and nitiaufin. Parameter grid stands for the configuration name : piano loads the namelists defined in config/grid. Four directories have to be declared in the description file, namely tdir, exe_dir, obs_dir, and forcing_dir. If needed, the name of the restart file, and the background normalization file, can be provided. The two last parameters can be usefull to : if the user wants to run only one executable, the parameter run_solo has to be specified. Also, it is possible to ask for piano to start the assimilation cycle from a previous one, stopped at iteration ndone.

Description of the directory piano

Files namelist_nemo and namelist_nemovar contains a subset of the full namelist of a NEMO/NEMOVAR/NEMOTAM run. Indeed, they contain only the modifications made by the user to the default namelists located in config/grid.

For instance :

> cat exp/GYRE-ENACT/namelist_nemo &nam_asminc ln_bkgwri = .true. nitbkg = 0 nitdin = 0 &namtam ln_trjwri = .true. &namobs ln_t3d = .true. ln_s3d = .true. ln_profb = .true. profbfiles='profb_01_fdbk_0000.nc' piano merges namelist_nemo files loaded from exp/GYRE-ENACT and config/GYRE_Z31, and generates the namelist file that will be used for the NEMO run. In this example, the experiment asks for NEMO to assimilate an increment from a NE-MOVAR run (ln_bkgwri), to generate the new trajectory (ln_trjwri), and to write the model equivalent of T and S profile observations (parameters set to true in &namobs and name of the observation file).

> cat exp/GYRE-ENACT/namelist_nemovar &namobs ln_prf=.true. ln_t3d=.true. ln_s3d=.true.

&namalg noutmax = 8 noutit = 1 nvarex = 1
Here the experiment parameters specify the kind of assimilated observations, the number of outer loop, of the current loop (noutit that can eventually be updated by ndone), the kind of assimilation (4D-var in this case, since nvarex=1).

Files namelist_nemo and namelist_nemovar in the exp directory can be empty : in this case, the namelist used by the run is the one defined in the asso-ciated configuration directory (for instance, GYRE_Z31 in the case of GYRE-ENACT experiment).

Running piano 1.2.1 Running piano -general features

In directory piano, run the command >python piano.py exp_name where exp_name is the name of an experiment. Depending on your experiment, several directories are created in tdir. If a full assimilation cycle is asked, then an outer and an inner directory are created, namely :

> cd $tdir inner_GYRE-ENACT_mac_intel_1x1 outer_GYRE-ENACT_mac_intel_1x1 outer_tam_GYRE-ENACT_mac_intel_1x1
If tangent and adjoint model tests are done (parameters ln_tst=ln_tst_tam = .true. in namelist_nemovar, for instance) then the output is located in the outer_tam directory. Else, if run_solo=nemo (resp. nemovar), only the outer (resp. inner) directory is created, and modified. If a full assimilation cycle is done (run_solo deactivated), then both inner and outer directories are created (or overwritten).

A run is described by two parameters of description file : ndone and run_solo : • run_solo=full_assim (default) : if ndone not defined clean all directories. Else set noutit=ndone.

• run_solo=nemo : If ndone not
In the three last cases, the NEMOVAR namelist is initialized. For nemovar or full_assim, check that nvarex=1 if you want a 4D-Var assimilation.

Running piano for Tangent Linear Hypothesis (HLT) test

HLT stands for Tangent Linear Hypothesis test. The user needs first t define its experiments through namelist_nemo and name-list_nemovar and descr in the experiment. (see NEMOTAM User's Manual Reference) More particularly, blocks namhlt and tl_tamtrj. The HLT test is activated if run_solo is defined as nemohlt.

For instance :

> cat exp/GYRE-HLT/descr compile=mac_g95 tdir=MY_PATH exe_dir=MY_EXE_PATH npx=1 npy=1 grid=GYRE_Z31 run_solo=nemohlt nittrjfrq_hlt=75 increments_hlt=MY_INCREMENT_PATH restart_file=MY_RESTART_PATH #restart_hlt=MY_2ND_RESTART_PATH #ndone=1 time_dir=MY_TIME_PATH common_dir=MY_COMMON_PATH
The increments δx is read from an increments file increments_hlt or computed from two restart files restart_file and restart_hlt. The choice is done through the namelist parameter ln_incdx from namhlt namelist block. The response time frequency for HLT is controlled with nittrjfrq_hlt (Note : it will orverwrite nittrjfrq when computing M(x + δx) and nittrjfrq_tan when computing L(δx). If ndone is equal to 1, the driver assumes that M(x) is already available in outer

Running piano for tangent test

The tangent test refers to the analysis of the tangent behaviour with respect to the direct module when the perturbation evolves as p n δX. If p n has a geometric progression, PIANO is able to perfom the loops on n and concatenate the final results into a single file tan_diag.global_0000 . To do so, the user must provide the element to compute p n = p 0 r -n :

-the initial value p 0 (tlm_start) -the common ratio r (tlm_fact) -the number of loops on n (tlm_loop)

The tangent test is activated if run_solo is defined as nemotam and ln_tst_tan is set to .true. (plus ln_tst_tan_cpd for instance) in namtst block namelist. For instance : The Forcing files are usually time dependent data. They will be linked or copy to the outer and inner directories (and if needed test directories). The directory is identified in the descr file under the name time_dir. The forcing files ( surface boundary conditions, ice cover, fluxes, surface wind) are defined in the NEMO namelist (namelists &namsbc_xxx, activated by the logical ln_xxx in namelist block &namsbc). Once the time_dir is specified in descr, the loading of all the files activated in these namelist blocks should be transparent. Note that for the moment only the forcing files for ORCA2 has been implemented : the use of other files (by then other parts of the namelist) has not been tested. Since information on files is given in the namelist, note that it is possible to load different forcing files for NEMO and NEMOVAR.

> cat exp/GYRE-XXX

Data files (common_dir)

We call Data files, files related to the configuration and general forcing file. The are usualyy used for realistic configuration such as ORCA2 and their file names are hardcoded into NEMO. The directory is identified in the descr file under the name common_dir. For the moment, they are automatically loaded for the ORCA2 grid, by mean of the method init_real in module init_run.py. def init_orca(descrlist,common_files): grid=descrlist.get_value('grid') if grid=="ORCA2_Z31": forcing_files.append("data_1m_potential_temperature_nomask.nc") forcing_files.append("data_1m_salinity_nomask.nc") forcing_files.append("dist.coast.nc") forcing_files.append("bathy_level.nc") forcing_files.append("bathy_meter.nc") forcing_files.append("ahmcoef") forcing_files.append("coordinates.nc") forcing_files.append("ssh.nc") forcing_files.append("geothermal_heating.nc") forcing_files.append("slaReference.nc") Note 1 : treatment of "slaReference.nc" does not rely on a specific configuration.

Note 2 : For the moment, piano works only for one cycle of assimilation : no specific treatment is done with respect to the date management, and so the storage is done in a very simple way. This question will be addressed when the implementation of several cycles of assimilation is done.

Observation files (obs_dir)

It is possible to load observation files of different format, as specified in the namelist &namobs, but the communication between NEMO and NEMOVAR is done only with feedback files.

Note that observations such as sla needs an additional set of data (that can be called slaReference.nc) to work.

3 Specific features

NEMOVAR tests

By default feedback files needed by nemovar (inner loop) is the output of a direct simulation (outer loop). In a inner loop test context it is possible to specify profiles and/or sla feedback files (key word : profile_fdbk_file and sla_fdbk_file respectively) when the key word nemovar is used for run_solo experiment. Then, your descr may contain the following lines : run_solo=nemovar profile_fdbk_file=$MY_PATH_TO_PROF_FB_FILE.nc sla_fdbk_file=$MY_PATH_TO_SLA_FB_FILE.nc

Change Assimilation method between outer loops (nvarex_list)

For testing purpose we offer the possibility to change the nvarex value (means 3DFGAT or 4DVAR assimilation method for nvar equals 0 or 1 respectively). It can be done adding the parameter nvarex_list and the corresponding list in your descr file. Parameter nvarex is a list a maximum noutmax (maximum number of outer loops) elements. Each element is separated by a comma (,). If the number elements is less than noutmax, nvarex is set to the default value specified in name-list_nemovar (or by default in the default configuration).

Then, if you want to perform 4 outer loops in your assimilation process using 3DFGAT method for the first and the second set of inner loops and then 4DVAR for the thrids and the fourth sets, your descr may contain the following lines : nvarex_list=0,0,1,1 4 Installation New Configuration

Input files

The management of input files for NEMO/NEMOVAR depends on the namelist, and on convention name defined in the code. There are several kind of input files :

Observation files : their name can either be defined by default in case of feedback files for profiles, or sla observations, or be present as arguments in the namelist block namobs.

Forcing files : surface boundary condition, ice cover, defined in the namelists namsbc_xxx : Their names are defined in the corresponding namelist block, when the flag in nam_sbc is activated. These kind of files are thus automatically treated in PIANO.

Data files (for instance the file slaReferences.nc needed when there is some SLA observation)

Files that are not SBC files, or that are about specific data, need a specific treatment in PIANO. Examples are available for ORCA2.

Available options in 'descr' file

Tables 1 and2 give a list of available options for the descr file.

Références

  The software PIANO is available in the VODA development repository in directory UTIL :>cd ~/VODA/UTIL/python/pianoThis directory contains all the python sources (files *.py), and two directories called exp and config

	>ls config/
	GYRE_Z31/
	ORCA2_Z31/
	POMME_Z46/
	SQB_Z11/
	Each directory contains two default namelists
	>ls GYRE_Z31/
	namelist_nemo
	namelist_nemovar
	>ls
	block.py	! Python class defining one namelist block
	blocklist.py	! Python class defining the whole namelist
	char_handling.py ! Extract strings from a files
	comments.py	! Python class defining a comment in the namelist
	config/	! Directory containing the available NEMO configuration
	descrlist.py	! Python class defining the description file
	exp/	! Directory containing the user experiments
	init_run.py	! Method for initialization of PIANO
	piano.py	! Main
	param.py	! Python class defining a parameter (name, value, comment)
	ps_handling.py	! Manages processes
	run_main.py	! Methods that manage the run of one executable
	hlt_handling.py	! Method to compute RMS, coeff. corr of lin.tgt hyp test
	npz2plot	! Method to convert hlt_handling output to plot
	tlmdiagmrg.py	! Method to merge outputs files from tangent tests
	Details on the python code architecture are given in section ??. The directory
	config contains the following configurations :

  The directory exp contains examples of user experiments :

	> cat exp/GYRE-ENACT/descr
	#
	# Description file for the Experiment
	#
	ndays=6	! to compute nitend and nwrite
	grid=GYRE_Z31	! Configuration
	compile=mac_intel	! Compiler
	npx=1	! Parallel run:
	npy=1	!	not implemented
	tdir=/tmp	! Contains output directories
	exe_dir=/tmp	! Contains nemofcm_build dir
	time_dir=/path_to_forcing_dir	! For forcing files
	common_dir=/path_to_common_dir	! For data files
	obs_dir=/path_to_obs_dir	! Observation files
	restart_file=/path_to/restart.nc ! If ln_rstart
	#bkgnorm_file=/path_to/background.normalization.nc
	#run_solo=nemovar	! Run only one executable
	#ndone=1	! Start assim cycle from a previous one
	#trj_dir=/path_to/tam_trajectories! When only nemovar run
	>ls exp/
	GYRE-ENACT
	GYRE-SGL
	ORCA2
	POMME
	An experiment is described by three files :
	>ls GYRE-ENACT/
	descr
	namelist_nemo
	namelist_nemovar
	The description file descr contains all the parameters needed by a run that are
	not in the namelists :

  defined in descr or ndone=1 : Delete outer directory. Outer loop counter noutit=ndone. • run_solo=nemotam : Delete outer_tam directory. Executable model_tam.exe is used to run tests on tangent and adjoint model of NEMO : check that logical for such tests are activated in namelist namtst : ln_tst=.true., and at least one of the following flags has to be set to .true. : If ndone not defined in descr or ndone=1 : Delete inner directory. Outer loop counter noutit=ndone.

	ln_tst_cpd_tam
	ln_tst_stp_tam
	ln_tst_tan
	• run_solo=nemovar :

TABLE 1 -

 1 Available options in descr file. Part 1/2


Pascal Cotret

Guy Gogniat

Jean-Philippe Diguet

Jérémie Crenne
email: jeremie.crenne@lirmm.fr

LIGHTWEIGHT RECONFIGURATION SECURITY SERVICES FOR AXI-BASED MPSOCS

Nowadays, security is a key constraint in MPSoC development as many critical and secret information can be stored and manipulated within these systems. Addressing the protection issue in an efficient way is challenging as information can leak from many points. However one strategic component of a bus-based MPSoC is the communication architecture as all information that an attacker could try to extract or modify would be visible on the bus. Thus monitoring and controlling communications allows an efficient protection of the whole system. Attacks can be detected and discarded before system corruption. In this work, we propose a lightweight solution to dynamically update hardware firewall enhancements which secure data exchanges in a busbased MPSoC. It provides a standalone security solution for AXI-based embedded systems where no user intervention is required for security mechanisms update. An FPGA implementation demonstrates an area overhead of around 11% for the adaptive version of the hardware firewall compared to the static one.

INTRODUCTION

Embedded systems are facing an increasing number of threats as attackers' motivation is raising every day. Our devices contain many sensitive information (passwords, private information) that needs to be protected from software and hardware attacks. Reconfigurable technologies such as FPGAs can be a good candidate to build trusted devices as they embed processors, memories and application-specific IPs in a single chip with moderate development costs. They also offer several interesting features to build high-performance high-security systems due to their intrinsic performance and adaptive properties. When dealing with logic attacks (e.g. targeting the external memory through code/data corruption) main existing solutions are based on software countermeasures. However, relying the system security on software-only solutions may not be adapted for high constrained embedded systems. We believe adding some hardware mechanisms strongly increases systems security at a very low cost. Monitoring and controlling the communication architecture in a bus-based MPSoC is particularly well adapted to enhance the security of a system as all data are exposed to its structure. Thus if an attacker tries to extract, destroy or modify some data through a logical attack he will have to use the communication architecture. Therefore, enhancing the communication architecture with protection mechanisms is a crucial point. However, when building such a solution several key points need to be addressed: what happens when an attack is detected? How the system should behave? Is it required to increase the protection level of the system? In this paper, we address these questions and propose security update features to provide designers dynamic security levels. The paper is organized as follows. Section 2 presents related work. Section 3 summarizes a static solution defined in [START_REF] Cotret | Distributed security for communications and memories in a multiprocessor architecture[END_REF]. Section 4 describes our contribution and Sections 5 and 6 propose several results and analysis. Section 7 highlights main perspectives.

RELATED WORK

In the literature, several studies have addressed the security of embedded systems [START_REF] Ravi | Security in embedded systems: Design challenges[END_REF]. At the communication level, these systems can be protected either by software or hardware mechanisms. Software solutions generally do not require additional hardware but offer low efficiency in terms of latency which can be critical for applications where reactivity is essential to fend off attacks. From an hardware point of view, several solutions have been proposed depending on the communication architecture technology: network-on-chip (NoC) or bus. Regarding NoC-based architectures, Evain et al. [START_REF] Diguet | Noccentric security of reconfigurable soc[END_REF] propose a solution where security controls are done in each network interface in a distributed manner. A management unit gathers all information from network interfaces according to a security policy. Fiorin et al. [START_REF] Fiorin | A data protection unit for noc-based architectures[END_REF] [5][6] propose an alternative to this work by adding probes within the interface structure to refine the protection mechanisms. These probes can block incoming traffic according to parameters stored in an embedded memory. A security manager collects information from individual security-enhanced interfaces to detect any collision or error in the traffic. They also provide a adaptive implementation by adding a shadow memory (which behaves like a buffer mechanism) to avoid a temporary state of the security enhancement during its reconfiguration. For bus-based communication architectures, the work of Coburn et al. [START_REF] Coburn | Seca: security-enhanced communication architecture[END_REF] which is similar to Fiorin's work (without updates) is based on SEI (Security Enforcement Interface) implemented in each interface between an IP and the bus. Each SEI computes information from the data handled by the IP and sends it to a global manager (SEM, Security Enforcement Module) which aims to check SEI data.

This work focuses on an FPGA implementation of a adaptive version of the work presented in [START_REF] Cotret | Distributed security for communications and memories in a multiprocessor architecture[END_REF]. It is inspired by [START_REF] Fiorin | A monitoring system for nocs[END_REF] with an implementation specific to the AXI bus standard with a low area/latency trade-off. Moreover, several scenarios and a case study are discussed to demonstrate all the features provided in this work.

The key contributions of this work include:

• Demonstration of adaptive firewall enhancements.

• Design of a set of security policies.

• Analysis of attack scenarios and case study implementation.

STATIC SECURITY ENHANCEMENTS WITHIN A BUS-BASED MPSOC SYSTEM

This work is based on static security enhancements proposed in [START_REF] Cotret | Distributed security for communications and memories in a multiprocessor architecture[END_REF] and [START_REF] Cotret | Bus-based mpsoc security through communication protection: A latencyefficient alternative[END_REF]. It provides a low-latency solution based on hardware firewalls embedded in each IP bus interface providing protection against read/write access and format disruptions (Local Firewalls). The firewall connected to the external memory controller (Cryptographic Firewall) adds flexible cryptographic services to protect the memory with confidentiality and/or authentication (Figure 1).

System bus

Local Firewall Proc.

Crypto Firewall External Memory

FPGA (trusted) untrusted

Mem. ctrl

Firewall Interface

Security Builder

Crypto Module BRAM_addr Corr. security policies contain cryptographic information (such as keys), read/write access and format rules for a given address space. Then, these parameters are sent to the Checking Module which compares the system bus parameters with the values extracted from the security policy; at this step, if cryptographic operations are needed, the Crypto Module manages the encryption/decryption and authentication tasks with a dedicated BRAM for cryptographic information storage (this is not the main point of this work). Once Checking Module completes its operations, a check out signal is sent to the Firewall Interface to confirm or not the data validity (security policies are verified or not). Finally, Firewall Interface provides the final data and manages synchronization tasks in order to fit with the output bus interface. Using this method, the system is protected against logical attacks aiming to tamper the external memory and the external bus without encrypting the whole memory which has a strong impact in terms of latency [START_REF] Cotret | Distributed security for communications and memories in a multiprocessor architecture[END_REF]. Unfortunately, this is a static solution where security policies cannot be changed. Security services update could be done with partial reconfiguration of firewalls or a complete reconfiguration of the FPGA. However these solutions would negatively impact the reaction time of the system in case of an attack and potentially compromise the security. This work focuses on the implementation of a solution based on Block RAM real-time modifications without downloading any new bitstream.

SECURITY SERVICES UPDATE REQUIREMENTS

When an attack event is detected, BRAM contents must be updated with new security policies in order to keep a safe execution environment for the target MPSoC. This work proposes an architecture for this purpose as detailed in Figure 2. All the components are connected through an AXI-Lite bus (also known as Security bus) and managed by a trustworthy processor (program stored in a trusted ROM) which stores important events in a log file readable by the processor running the main application (timestamps, attack events. . . ). Each firewall has a custom bus connection with the monitoring IP. While log timer and monitoring IP are used for detecting and reporting an attack, AXI-Lite BRAM controller connections are used for security policies update.

The first task accomplished in the architecture presented in Figure 2 is the monitoring of attack events by a dedicated IP in the left hand side of the figure. Once an attack is detected, register values are update and an interruption routine is launched for security update on the trustworthy processor depending on the main register value (representative of all attack events). On an attack event, the system security must be updated in order to avoid malicious data leakage. One very important issue during firewalls update is the data availability while switching between two security policies. For this purpose, a mechanism based on the handshake property of the AXI protocol is implemented in each Firewall Interface to block outcoming traffic while updating security rules associated with the current firewall. Two classes of components are defined based on their ability to manipulate confidential information. Critical IPs (for instance, ciphering algorithms implementation) must not reveal any information when an attack is detected. Extracting keys and/or signatures would be a major threat for the system. In that case as soon as an attack is detected critical IPs are isolated from the system (known as "error mode" or quarantine).

For non-critical IPs, an intermediate protection layer is authorized where reading accesses are still allowed but no writing ones. This feature aims for example to allow a backup of data before IP isolation.

In case of an attack event (detected through interruption routines launched by the monitoring IP), the trustworthy processor saves the current security policy of the attacked firewall in a dedicated on-chip memory and applies a higher security level according to the two schemes previously defined. For instance, if the current protection level of an IP is a "readonly mode", the next one could be an "error mode", equivalent to a quarantine feature.

In the most critical case (i.e. an attack event is detected even if the IP being monitored is under the "error mode" protection configuration), it is assumed that a reboot of the complete system is required. Therefore, the initial bitstream and security policies configuration are reloaded. The system restarts from an initial state.

Switching between the different modes

When an IP firewall must be updated, Security Policies parameters stored in Block RAMs have to be modified (assuming that the whole IP address space is covered by one or more SP, the only component to be updated is the BRAM containing SP values). This work considers only the update of read/write rights (fields 1 in Figure 3). Local Firewall SP is stored on a single 32-bit block while Cryptographic Firewall read/write information is stored in the first 32-bit word in a 32-bit word) in a BRAM takes one clock cycle. Therefore, the update of N Security Policies in one firewall is done in N clock cycles. This time has no real impact on firewall data analysis because data is blocked as soon as an attack is detected (previously described in this work).

IMPLEMENTATIONS RESULTS

Implementations were done using Virtex-6 Xilinx FPGA technology (model xc6vlx240t1156-1). This device has around 240,000 logic cells and 15 Mb of Block RAM. First, different implementation options for firewalls are studied.

Then, this work focuses on generic scenarios that serve as a basis for further case study analysis.

Area

In Table 2, two implementation options are considered: static Local Firewall (based on results of [START_REF] Cotret | Distributed security for communications and memories in a multiprocessor architecture[END_REF]) and the adaptive version. The adaptive version takes into account two features presented in this work: the run-time update and the error mode option. The area overhead due to these enhancements is low (around 11%). The update is based on simple logic elements and the use of bus properties, there is no need to store data in a buffer mechanism while updating security policies of a firewall.

Slices

Latency

Update latency of security policies depends on the number of policies to be updated. The most critical step in this process is the computation of the new security policy configuration done by the trustworthy processor. For a basic implementation of this process, the new configuration is computed in 148 clock cycles. When more than one attack has to be reported (N firewalls detect an attack), they are first blocked according to their order in the main register of the Monitoring IP. The trustworthy processor computes all the up-to-date security policies. Finally, each firewall is released as soon as its update is completed. Firewalls being updated are not affected because they are not able to transmit any information (handshake signals are controlled by the update logic). This method allows any updated firewall to be used by the main application of the MPSoC system. The update latency depends on the firewall location in the update queue: the first firewall will be modified in N cycles, while the firewall placed in location #k in the queue is updated in k(N ) cycles because it must wait for the first k -1 firewalls to be updated.

CASE STUDY

The MPSoC case study considers 2 Microblaze processors, a 64KB shared Block RAM, an image processing IP and an external memory. Each processor has code and data sections in an external memory. This architecture is set with mixed cryptographic options and access rights to get all the options (integrity only, read/write, confidentiality and integrity. . . ). 4 Local Firewalls and 1 Cryptographic Firewall are needed for the protection of this case study. Three options are considered: the MPSoC without firewalls, the static firewallenhanced MPSoC (based on the results of [START_REF] Cotret | Distributed security for communications and memories in a multiprocessor architecture[END_REF]) and the adaptive version with firewalls. Area results are summarized in Table 3. The firewall-enhanced case study has a quite high overhead: this is mainly due to the cryptographic module embedded in the firewall attached to the external memory controller. The logic added for update purposes implies a quite high area overhead of around 40% compared to the static firewall implementation: this is mainly due to the trustworthy processor, the monitoring IP and the security AXI-Lite bus.

Slices

CONCLUSION AND PERSPECTIVES

In this work, a bus-based MPSoC with adaptive security enhancements (also known as firewalls) is presented. These firewalls protect memories and memory-mapped IPs according to user-defined security policies. By using communication bus properties, it allows developers to get run-time updates with a low area overhead compared to a static solution; furthermore, when security updates are required, there is no invalid data leakage. Mechanisms defined in this work do not need any user event as protection levels are all defined in a dedicated processor. This work corresponds to a trade-off between SECA [START_REF] Coburn | Seca: security-enhanced communication architecture[END_REF] and DPU [START_REF] Fiorin | Implementation of a reconfigurable data protection module for noc-based mpsocs[END_REF] solutions with an implementation on the AXI bus standard implemented in Xilinx tools suite. This work proposes an update feature which is not present in the other bus-based solution. In terms of area overhead over a basic softcore processor, adaptive firewalls have a lower impact (around 11%) than DPU solution (25%, it is based on a shadow memory acting as a buffer-like mechanism); SECA solution has the lowest area overhead (6.20%, mainly due to the centralized security manager) but provides a static solution.

Fig. 2 .

 2 Fig. 2. Architecture of the security/monitoring area

Fig. 3 .

 3 Fig. 3. Security Policy memory layout for Local and Cryptographic Firewalls

Table 1 .

 1 Comparison with existing works Fiorin et al. approach

		SECA	DPU	This work
		[7]	[5]
	Communication	AHB/APB	NoC	AXI-4
	standard
	Adaptive ?	No	Yes	Yes
	Threat	Wide range of	Mainly buffer Wide range of
	model	soft. attacks	overflow	soft. attacks

Table Checking

 Checking Firewall Interface while other information (such as address, format and read/write modes) are sent to the Security Builder. The Correspondence Table indicates the location of the security policy associated with the bus address:

		Proc.	Proc.
		Local	Local
		Firewall	Firewall
	Local	Local	Local
	Firewall	Firewall	Firewall
	IP	IP	IP
					Specific to
	Block		Security Policies		Cryptographic Firewall
	RAM
				Trusted	Untrusted
	AXI4_addr
	AXI4_format
	AXI4_rnw	Module
		AXI4_data	check_out
	AXI
	system bus	AXI-4	AXI-4	DDR Memory IP or

Fig. 1. Firewall structure within a bus-based MPSoC

When a data comes from the AXI system bus, it is stored in the

Table 2 .

 2 Standalone firewalls results

Table 3 .

 3 Table of standalone results

Acknowledgment

The work presented in this paper was realized in the frame of the SecReSoC project number ANR-09-SEGI-013, supported by a grant of the French National Research Agency (ANR).