domain: info.info-au math spi sdv

text:
Madalena Chaves
email: madalena.chaves@inria.fr

Alfonso Carta
email: alfonso.carta@inria.fr

Attractor computation using interconnected Boolean networks: testing growth rate models in E. Coli

Boolean networks provide a useful tool to address questions on the structure of large biochemical interactions since they do not require kinetic details and, in addition, a wide range of computational tools and algorithms is available to exactly compute and study the dynamical properties of these models. A recently developed method has shown that the attractors, or asymptotic behavior, of an asynchronous Boolean network can be computed at a much lower cost if the network is written as an interconnection of two smaller modules. We have applied this methodology to study the interconnection of two Boolean models to explore bacterial growth and its interactions with the cellular gene expression machinery, with a focus on growth dynamics as a function of ribosomes, RNA polymerase and other "bulk" proteins inside the cell. The discrete framework permits easier testing of different combinations of biochemical interactions, leading to hypotheses elimination and model discrimination, and thus providing useful insights for the construction of a more detailed dynamical growth model.

Introduction

Large networks with complex interactions are hard to analyze in detail, but logical and discrete models can facilitate this task. Based essentially on the structure and topology of the network interactions, logical models provide qualitative information on the dynamical properties of system [START_REF] Glass | The logical analysis of continuous, nonlinear biochemical control networks[END_REF][START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF], which can be used for model discrimination and guidance in model improvement. There are many recent examples of applications of discrete models including Drosophila embryo pattern formation [START_REF] Sánchez | A logical analysis of the drosophila gap-gene system[END_REF][START_REF] Albert | The topology of the regulatory interactions predicts the expression pattern of the Drosophila segment polarity genes[END_REF], yeast cell cycle [START_REF] Sevim | Reliability of transcriptional cycles and the yeast cell-cycle oscillator[END_REF], T-cell response [START_REF] Saez-Rodriguez | A logical model provides insights into T cell receptor signaling[END_REF], or an apoptosis network [START_REF] Calzone | Mathematical modelling of cell-fate decision in response to death receptor engagement[END_REF].

Boolean networks are a class of logical models whose variables are described in terms of only two levels (1 or 0; presence or absence; "on" or "off"), which have been useful for biochemical systems [START_REF] Wang | Boolean modeling in systems biology: an overview of methodology and applications[END_REF]. The dynamics of a Boolean model is determined by specifying an updating mode, most commonly synchronous (all nodes updated simultaneously) or asynchronous (only one node updated at any given instant). Since the state space is finite, the dynamics can be represented in terms of a transition graph, which can be studied using some classical algorithms from graph theory [START_REF] Cormen | Introduction to algorithms[END_REF]. Other, more specific tools are available for an exact and rigorous analysis of the transition graph [START_REF] Lorenz | Analysis and characterization of asynchronous state transition graphs using extremal states[END_REF], computation of attractors (or asymptotic behavior) [START_REF] Chaves | Predicting the asymptotic dynamics of large biological networks by interconnections of Boolean modules[END_REF], and other properties [START_REF] Thomas | Biological feedback[END_REF]. In addition, a wide range of computer tools are available for simulation and analysis of discrete models [START_REF] Gonzalez | GINsim: a software suite for the qualitative modelling, simulation and analysis of regulatory networks[END_REF], model reduction [START_REF] Naldi | Dynamically consistent reduction of logical regulatory graphs[END_REF], or model checking [START_REF] Fages | Modelling and querying interaction networks in the biochemical abstract machine BIOCHAM[END_REF].

It is clear that discrete models are not appropriate to finely describe the behavior of a system, since continuous effects are difficult to reproduce with such models (such as whether an oscillation is sustained or damped), but they are useful to verify whether a given network of interactions is feasible and compatible with known properties of the system. This is a first step towards the construction of a more detailed and informed model.

As an application, we will analyze a network of interactions involved in determining bacterial growth of Escherichia Coli, which varies nonlinearly with different factors, such as availability of nutrients or the concentration of the necessary enzymes and proteins needed for cell division [START_REF] Bettenbrock | Analysis of the correlation between growth rates, EIIACrr phosphorylation, and intracellular camp levels in Escherichia coli K-12[END_REF][START_REF] Shachrai | Cost of unneeded proteins in E. coli is reduced after several generations in exponential growth[END_REF]. Mathematical models have been developed to describe and reproduce several regulatory modules and their response to nutrient availability [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF][START_REF] Hardiman | Topology of the global regulatory network of carbon limitation in Escherichia coli[END_REF]. One of the least understood aspects in these studies remains the actual modeling of bacterial growth: while it is clear that growth depends on the general availability of "bulk" proteins, ribosomes, and RNA polymerase, it is difficult to find a reasonable mathematical model that reproduces all these effects [START_REF] Goelzer | Bacterial growth rate reflects a bottleneck in resource allocation[END_REF]. In many cases, growth is considered to be a given constant and the model is designed to reproduce a single phase of bacterial growth.

Here, we propose to test and study a dynamical function for bacterial growth in terms of the major components involved in bacterial cell division, that is, gene transcription (RNA polymerase) and translation (ribosomes). To test the feasibility of mathematical growth functions, we will focus on a qualitative model of the network involved in the carbon starvation response [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF] and its interconnection with a basic model describing the dynamics of ribosomes and RNA polymerase (see Section 3).

We will use two methods for analysis of qualitative systems (see Section 2): first, a method that transforms piecewise affine (PWA) systems into discrete and then Boolean models [START_REF] Chaves | Comparing Boolean and piecewise affine differential models for genetic networks[END_REF][START_REF] Jamshidi | Comparing discrete and piecewise affine differential equation models of gene regulatory networks[END_REF]; and, second, a recently developed method to compute the attractors of an interconnection of two Boolean modules [START_REF] Chaves | Predicting the asymptotic dynamics of large biological networks by interconnections of Boolean modules[END_REF][START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF]. Our analysis generates a general view of the dynamical properties of a model which is a first step towards verifying the feasibility of the model's structure -by comparing to experimental observationsand facilitates hypotheses testing. The results indicate that at least two (positive) qualitative levels for growth rate (such as "high" and "intermediate" rates) are needed in order to reproduce both the stationary and exponential growth phases (see Section 4).

Methodology

In this section, we briefly recall two mathematical methods which are very useful for the analysis of qualitative systems and, in particular, interconnections of Boolean models.

From discrete to Boolean models

Although Boolean variables can only take the values 0 or 1, it is nevertheless possible to construct Boolean models that describe variables with a discrete number of values [START_REF] Van Ham | How to deal with variables with more than two levels[END_REF][START_REF] Chaves | Comparing Boolean and piecewise affine differential models for genetic networks[END_REF]. Consider a discrete model Σ disc = (Ω d , F d ), with variables V = (V 1 , . . . , V n ) ′ , state space Ω d = Π n i=1 {0, 1, 2, . . . , d i }, where d i ∈ N is the number of levels of variable V i (i = 1, . . . , n), and a state transition table

F d : Ω d → Ω d .
The state of the system at the next instant k + 1 is given as a function of the state of the system at the current instant k, according to the rules F d , using the notation:

V + = Fd (V ).
Throughout this paper, the function Fd is obtained from F d by assuming an asynchronous dynamic updating rule, that is, exactly one variable is updated at any given time:

V + ∈ {W ∈ Ω d : ∃k s.t. W k = (F d ) k (V ) = V k and W j = V j , ∀j = k }. (1)
Furthermore, for a more realistic model, we consider that each variable V i can only switch from its current level to an immediately adjacent level [START_REF] Thomas | Biological feedback[END_REF], that is:

V + i ∈ {V i -1, V i , V i + 1}, ∀i. (2)
The idea is to create an extended Boolean model Σ bool = (Ω b , F b ) where each discrete variable V i is represented by d i Boolean variables, for instance, {X i,1 , . . . , X i,di }, so that the state space of the model is

Ω b = {0, 1} d1+•••+dn .
There are several possible ways to convert the discrete into the Boolean variables, but here we chose to use the same criterion as in [START_REF] Chaves | Comparing Boolean and piecewise affine differential models for genetic networks[END_REF] which stipulates that

V i = k ⇔ (X i,1 = • • • = X i,k = 1, X i,k+1 = • • • = X i,di = 0), (3)
meaning that a variable i is at a state k if and only if all the first k Boolean variables are ON. In particular, note that this criterion implies the partition of the state space of the extended Boolean model into permissible and forbidden regions:

Ω p = {X ∈ Ω b : k < l ⇔ X i,k ≥ X i,l }, Ω f = {0, 1} d1+•••+dn \ Ω p .
Thus, to generate the Boolean transition table F b we need to guarantee that no transitions from a permissible to a forbidden state take place. The method described in [START_REF] Chaves | Comparing Boolean and piecewise affine differential models for genetic networks[END_REF] deals with this problem in a natural way, and guarantees that no transitions from permissible to forbidden states take place.

Dynamics of Boolean models

This section contains a brief summary of some useful objects that characterize the dynamics of a Boolean model. There are several possible ways of defining the dynamical updating rules [START_REF] Wang | Boolean modeling in systems biology: an overview of methodology and applications[END_REF] of a Boolean network Σ = (Ω, F b ), but here we will assume asynchronous updates, so the definitions and rules (1) stated for discrete systems also apply, with d i = 1 for all i. Note that (2) is immediately satisfied for Boolean models.

The asynchronous transition graph, G = (Ω, E), of system Σ is a directed graph whose vertices (or nodes) are the elements of Ω, and the edges are given by E. There are thus 2 n nodes in G. Given any two elements a, ã ∈ Ω the edge "a → ã" is in E iff:

ã ∈ {w ∈ Ω : ∃k s.t. w k = (F b ) k (a) = a k and w j = a j , ∀j = k }. A path a 1 ❀ a 2 in G is a sequence of edges linking a 1 to a 2 .
A strongly connected component (SCC) of G is a maximal subset C ⊂ Ω, that contains a path joining any pair of its elements. In general, a SCC may have both incoming and outgoing edges. An SCC with no outgoing edges is called terminal.

An attractor A of G is a terminal strongly connected component, that is, once a trajectory enters A it cannot leave again. Therefore, the attractors can be said to characterize the asymptotic behavior of the network. An asynchronous transition graph always has at least one, but can have multiple, attractors. An attractor can be formed of a single state (we will call it a singleton) or of a subset of Ω.

Interconnection of Boolean models

To study the interconnection of the two systems, we will use a method based on control theory concepts recently developed by one of the authors [START_REF] Chaves | Predicting the asymptotic dynamics of large biological networks by interconnections of Boolean modules[END_REF][START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF]. This method analyzes the asymptotic behavior of the interconnection of two systems directly from the behavior of the two subsystems, without having to construct or analyze the full interconnected system. The advantage is a much reduced computational cost, while still obtaining exact results: indeed, for large (e.g., n ≥ 15) Boolean models, the computation of the asynchronous transition graph and its attractors is unfeasible, as it involves the analysis of a 2 n × 2 n matrix. The idea is to first study each individual system for each set of inputs, obtain the corresponding attractors, and then construct a new object, the asymptotic graph. This new graph is much smaller than the state transition graph of the full model, but it contains all the information on its asymptotic dynamics, namely all the attractors of the full model correspond to attractors in the asymptotic graph. Some notation is next introduced.

Consider two asynchronous Boolean models, Σ A and Σ B , with a set of inputs (U i ) and a set of outputs (H i ):

Σ A = (Ω A , U A , H A , F A ) : Ω A = {0, 1} n A , U A = {0, 1} p A , H A = {0, 1} q A , Σ B = (Ω B , U B , H B , F B ) : Ω B = {0, 1} n B , U B = {0, 1} p B , H B = {0, 1} q B .
The following notation will be used: a ∈ Σ A and b ∈ Σ B denote the states of each system, u ∈ U A and v ∈ U B denote the inputs, and the output corresponding to state a will be represented by h A (a) ∈ H A (resp., h B (b) ∈ H B for state b). The synchronous rules are written:

a + = F A (a; u), and b + = F B (b; v).
For each fixed u ∈ U A , there is a set of attractors of system Σ A , and its elements will be represented by A i u , i ∈ N. Similarly for system Σ B , B j v , j ∈ N. The interconnection of these two systems is formed by letting the input of each system be the output of the other

v = h A (a) ∈ U B u = h B (b) ∈ U A ,
where it is assumed without loss of generality that q A = p B and q B = p A . The new system will be represented by: Σ = (Ω, F bool ) : Ω = {0, 1} n A +n B , F bool : Ω → Ω with the Boolean rules F bool given by the appropriate combination of F A , F B :

F bool (a, b) = (F A (a; h B (b)), F B (b; h A (a))).
Note that F A , F B , and F bool contain the synchronous table of state transitions. Here, we will consider that the dynamics is asynchronous, so that only one variable is updated at a given time. The asynchronous transition graphs of the two modules (one for each fixed input) and that of the full interconnected system will be called, respectively, G A,u , G B,v , and G.

Transition graphs and semi-attractors The first step of the method is to compute all the transition graphs G A,u and G B,v , compute their attractors, and then divide each of these into subsets corresponding to a fixed output. These will be called semi-attractors of the individual system and are defined as follows:

A i uα = the i-th semi-attractor of system Σ A , corresponding to input u, with output α B j vβ = the j-th semi-attractor of system Σ B , corresponding to input v, with output β.

Note that the standard attractor is the union of all corresponding "semi-attractors":

A i u = ∪ all α A i uα .
The asymptotic graph The second step of the method is to construct the asymptotic graph G as whose nodes are the cross-products of semi-attractors:

A i uα × B j vβ .
There is an edge between two of the nodes

A i uα × B j vβ → A i uα × B j α β
if there is a path in the graph G B,α that leads from some state in B j vβ to some state in B j α β . Similarly for an edge A i uα × B j vβ → A ĩ β α × B j vβ . In order to satisfy an asynchronous updating scheme, only one set of variables is allowed to change for each edge. The computational cost can be further reduced by observing that all nodes with u = β and v = α are transient (shown in [START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF]); hence, to compute the attractors of the asymptotic graph we only need to include the edges between nodes satisfying either u = β or v = α.

Attractors of an interconnection

The third step of the method is to compute all the attractors of G as which contain, in fact, a representative of each of the attractors of G. This is theoretically proven in [START_REF] Chaves | Predicting the asymptotic dynamics of large biological networks by interconnections of Boolean modules[END_REF][START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF]:

Theorem 1 [START_REF] Chaves | Predicting the asymptotic dynamics of large biological networks by interconnections of Boolean modules[END_REF] If Q is an attractor of G, then there exists at least one corresponding attractor in

G as , Q as = Q as (Q). Moreover, if Q 1 = Q 2 are two distinct attractors of G, then Q as (Q 1 ) = Q as (Q 2 ).
In other words, we recover all the attractors of the interconnection, without explicitly constructing the interconnected system. In broad terms, Theorem 1 says that any attractor of G generates an attractor in G as , but the converse is not necessarily true and G as may have more attractors than G.

To better illustrate Theorem 1, and show its advantages as well as limitations, a purely theoretical example is next given. For convenience, in the following examples, the attractors are labeled using the decimal representation for the Boolean inputs and outputs, that is:

000 ⇔ u = 1, 001 ⇔ u = 2, . . . , 111 ⇔ u = 8, etc.
2 )], h B (b) = b 2 ,
whose asynchronous transition graphs G A,u and G B,v are shown in Fig. 1, for convenience. Note that the attractors in all graphs are singletons except for B 2 2,2 = {01, 11}. However, since the two states have the same output (h B (01) = h B (11) = 1), in this example the semi-attractors are in fact the actual attractors.

The corresponding asymptotic graph is shown in Fig. 2. To illustrate the computation of an edge, consider the product A i uα × B j vβ = A 1 21 × B 1 11 : since α = 1 = v, the system A does not induce any change in the variables b; in contrast, the fact that β = 1 will induce a trajectory between a state in A 1 21 and an attractor in the graph G A,2 (corresponds to Boolean input u=1). In the graph G A,2 , the state 00 is in the basin of attraction of {01} = A 1 12 . Therefore, there is an edge

A 1 21 × B 1 11 → A 1 12 × B 1 11
. All other edges are similarly computed.

Note that the full interconnected system has four variables and hence its dynamics is given by an asynchronous transition graph G with 2 4 = 16 states. To compute the attractors of G we needed to compute a transition graph with only 2 × 3 = 6 states (2 attractors from system A and 3 from system B). Furthermore, as remarked above, the size of G as can be further reduced by excluding the crossproduct states known to be transient. In this example only A 1 21 × B 1 21 satisfies the condition u = β and v = α, and can be excluded. For higher order systems, such a size reduction can represent very significant savings in computational cost.

The G as of Example I has two attractors:

Q 1 = {A 1 12 × B 1 21 } and Q 2 = {A 1 12 × B 1 11 , A 1 12 × B 2 22 , A 1 21 × B 1 11 , A 1 21 × B 2 22 }.
For this 4-dimensional example, it is easy to check that Q 1 is a true attractor of the full interconnected system (see also Prop. 1), while Q 2 is a "spurious" attractor, that is, not a real attractor of G. To see this, it suffices to note that there is a pathway that leads from a state within Q 2 to Q 1 , and which is not "covered" by G as :

Q 2 ∋ A 1 21 × B 2 22 ∋ (00, 01) G B,1
-→ (00, 00)

G A,1 -→ (01, 00) ∈ Q 1
This Example shows that even very simple (and deterministic) individual asynchronous dynamics can lead to asymptotic graphs that exhibit spurious attractors. However, note that this example was specifically contrived to illustrate the generation of spurious attractors; its Boolean rules are not necessarily biologically plausible. In view of Example I, it would be useful to complement Theorem 1 by conditions permitting to decide whether an attractor of G as is also an attractor of G. An exact result was also proved in [START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF] -i.e., recovering exactly all the attractors in G from the cross-products of semi-attractors, with no spurious generation,-by computing the cross graph which is similar to G as but involves cross-products of all (semi-)SCCs (as opposed to considering only semi-attractors). However, depending on the number of SCCs, the cross graph can often be more costly to compute than the full graph G, hence the usefulness of establishing sufficient conditions for deciding whether an attractor of G as is a "true" attractor.

Some preliminary results were presented in Proposition 1 of [START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF], which are improved below in Prop. 1. To state this, we need to introduce projection functions, for V = A i uα × B j vβ , and R = {V 1 , . . . , V r }:

π(V ) = {(a, b) ∈ Ω : a ∈ A i uα , b ∈ B j vβ }, π(R) = ∪ V ∈R π(V ), π A (R) = {a ∈ Ω A : ∃b such that (a, b) ∈ π(R)}.
The A-output of R is the set:

A-output = {h A (a) : a ∈ π A (R)} ⊂ H A
Similar definitions apply for the projection π B (R) and the B-output of R.

Recall that we are assuming q A = p B and q B = p A , hence H A ≡ U B and H B ≡ U A and the A-output (resp., B-output) of R is also contained in U B (resp., U A ). The new result of Prop. 1 is in parts (ii), (iii), which previously stated "for all u ∈ U A " or "for all v ∈ U B ". The new conditions are much less restrictive, although the proof is similar. If Proposition 1 is not applicable, then one may still verify a posteriori whether R represents an attractor of G by simulating all trajectories starting from all states in π(R) and checking whether any of them leaves R (however, this "direct force" procedure may also involve some computational costs).

Proposition 1 Let R be a terminal SCC of G as . If either one of the following conditions is satisfied: i) R is a singleton (i.e., contains a single product V ); ii) the A-output of R is a singleton and the set π A (R) is an attractor of G A,u for all u in the B-output of R;

iii) the B-output of R is a singleton and the set π B (R)} is an attractor of G B,v for all v in the A-output of R;

then R as = π(R) is an attractor of G.

Proof: We will use the notation (a, b) ❀ G (a ′ , b ′ ) to denote a path connecting the two elements in the transition graph G and (a, b) → G (a, b ′ ) to denote a one-step transition. Part (i) is unchanged from [START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF]. Parts (ii) and (iii) are very similar, so we will only prove part (iii). If the B-output of R is a singleton, say {α}, then any V ∈ R must be of the form

A j α(•) × B (•)
vα , for some v in the A-output of R.

In particular, (see definition of semi-SCC) all A j α(•) belong to the same attractor A j α of G A,α . Suppose now that the set π B (R) is an attractor for all v in the A-output of R. Then, to show that π(R) is an attractor of G as , it suffices to show that: (1) π(R) is a strongly connected set, and (2) π(R) contains all its successors. If (1) and (2) hold, then π(R) is indeed a terminal SCC.

To show (1), let (a, b) and (a ′ , b ′ ) be any two elements of π(R).

Then (a, b) ❀ G (a, b ′ ), since π B (R) is an attractor of G B,h A (a) (v = h A (a) ∈ A-output) (a, b ′ ) ❀ G (a ′ , b ′ ), since a, a ′ belong to the same attractor A j α of G A,α .
To show (2), observe that there are two forms of successors:

either (a, b) → G (a ′ , b) or (a, b) → G (a, b ′ ).
We want to prove that both (a ′ , b) and (a, b ′ ) are in π(R). In the first case, since a, a ′ belong to the same attractor

A j α , it is immediate to see that (a ′ , b) ∈ π(R). In the second case, since b ′ ∈ π B (R) and π B (R) is an attractor of G B,h A (a) , by definition of π B (R) there some exists a ′ such that (a ′ , b ′ ) ∈ π(R). Recall that the B-output is a singleton so h B (b ′ ) = α. This implies (a, b ′ ) ❀ G (a ′ , b ′ ) ❀ G (a, b ′ ), since a, a ′ belong to the same attractor A j α of G A,α .
Therefore, (a, b ′ ) ∈ π(R) as wanted.

Remark. The generalization of points (ii) and (iii) of Proposition 1 to multiple A-outputs and B-outputs is not clear, due to Example I where the spurious attractor Q 2 satisfies A-output=B-output={1,2}. Other examples exist where an attractor of G as of the same form as Q 2 is indeed and attractor of G (see Example 2 in [START_REF] Chaves | Predicting the asymptotic dynamics of large biological networks by interconnections of Boolean modules[END_REF]).

If Proposition 1 cannot be applied, there may be other methods to decide whether an attractor of G as is a true attractor, such as identifying invariant sets of the system that contain the given attractor: examples of this are given below in Propositions 3 and 4.

Example II. To illustrate the relevance of Prop. 1, another theoretical example is now given. The two systems A and B are more conveniently represented by their asynchronous transitions graphs, one for each fixed input (Fig. 3). The dimensions are n a = 2, n B = 3, p A = 1, p B = 2 and their outputs are as follows:

h A (a) = (a 1 , a 2 ) ′ , h B (b) = b 1 .
Note that attractor A 1 2 splits into two semi-attractors, A 1 21 and A 2 23 , and the attractor B 1 2 splits into B 1 21 and B 2 22 . The full interconnected system has five variables and hence its dynamics is given by an asynchronous transition graph G with 2 5 = 32 states. To compute the attractors of G we needed to compute a transition graph with 4 × 7 -8 = 20 states: 4 semi-attractors from system A and 7 from system B, and 8 transient cross-products (see also Fig. 4).

The G as of Example II (Fig. 4) has two attractors:

Q 1 = {A 1 11 × B 1 11 } and Q 2 = {A 1 21 × B 1 32 , A 1 21 × B 2 12 , A 2 23 × B 1 12 , A 2 23 × B 1 32 }. It is easy to check that Q 1 = {00000} is an attractor of G, by Prop. 1(i). Likewise Q 2 = {10111, 10101, 10100, 00111, 00101, 00100}
is also an attractor of G, by Prop. 1(iii): the B-output is a singleton since {h B (b

) : b ∈ π B (Q 2 )} = {1}; the A-output of Q 2 is {h A (a) : a ∈ π A (Q 2 )} = {10, 00}; and, finally, the set π B (Q 2 ) = {111, 101, 100}
is indeed an attractor of both G B,v=10 and G B,v=00 . 3 Application: a model for E. Coli growth mechanism

The bacteria Escherichia Coli are unicellular micro-organisms (present in the human gut, for instance) which grow and divide in the presence of a carbon source, such as glucose or other sugars. In typical experiments, in a carbon rich medium, the bacteria are observed to grow at a constant growth rate, which is referred to as the exponential phase [START_REF] Marr | Growth rate of Escherichia coli[END_REF]. In the absence of carbon, the bacteria enter a stationary phase, with no cellular growth or division. E. Coli use a network of genes and proteins to detect the presence or absence of carbon sources and respond accordingly, by adjusting their gene expression levels.

The major players in this nutritional response network are well characterized (see, for instance, [START_REF] Bettenbrock | Analysis of the correlation between growth rates, EIIACrr phosphorylation, and intracellular camp levels in Escherichia coli K-12[END_REF][START_REF] Hardiman | Topology of the global regulatory network of carbon limitation in Escherichia coli[END_REF][START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF] and references therein) but, in contrast, it has been difficult to find an appropriate dynamical expression for modeling the growth rate of E. Coli [START_REF] Marr | Growth rate of Escherichia coli[END_REF]. In other words, if one wishes to add a model variable to describe growth rate, what should its mathematical rule be? To overcome this problem, models often focus on either the exponential or the stationary phases, thereby considering growth rate to be either constant or zero, respectively [START_REF] Goelzer | Bacterial growth rate reflects a bottleneck in resource allocation[END_REF]. However, such models are not able to describe the transition from one phase to the other, thus failing to provide intuition on a crucial cellular mechanism.

Growth should depend on the capacity of the bacteria to produce all the different proteins necessary to its development and cellular division. In its turn, the synthesis of any protein depends on the transcription and translation steps, which are limited, respectively, by the concentrations of RNA polymerase and ribosomes. To model the many proteins involved in bacterial growth, we will therefore distinguish between three "classes": RNA polymerase, ribosomal proteins, and all others will be collectively denoted as "bulk" proteins (as a reference see also [START_REF] Shachrai | Cost of unneeded proteins in E. coli is reduced after several generations in exponential growth[END_REF], where a distinction is made between ribosomal and nonribosomal proteins). Some models have thus tried to include these effects to obtain a more accurate expression for growth rate. For instance, one may have a dependence on one step: Growth rate ∼ RNA polymerase [START_REF] Sevim | Reliability of transcriptional cycles and the yeast cell-cycle oscillator[END_REF] as tested previously in [START_REF] Carta | A simple model to control growth rate of synthetic E. coli during the exponential phase: model analysis and parameter estimation[END_REF], or in two (or more) steps, each of them separately limiting growth rate, hence the use of the minimum function:

Growth rate ∼ min{ ribosomal proteins, bulk proteins }

as considered in [START_REF] Shachrai | Cost of unneeded proteins in E. coli is reduced after several generations in exponential growth[END_REF], or

Growth rate ∼ min{ ribosomal proteins, RNA polymerase }

as we considered in [START_REF] Carta | A class of switched piecewise quadratic systems for coupling gene expression with growth in bacteria[END_REF]. In this Section, our goal is to test these expressions, by interconnecting a well known nutritional response module with a basic transcription/translation model, using the Boolean interconnection method described in Section 2.3.

E. Coli nutritional stress response module

The nutritional stress response network developed in [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF] involves three groups of variables, each representing a different regulatory effect: DNA supercoiling (determined by the enzymes GyrAB and TopA), carbon response (involving the proteins Crp, Cya), and a global regulator (protein Fis) that sends the carbon availability signal down to the stable RNAs (rrn). The latter are limiting factors in ribosome production, and are thus a measure of the growth of the bacteria.

A very brief description of the main biological steps in response to nutritional stress is as follows (see [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF] and references therein): in answer to carbon depletion, the bacteria increase their cyclic AMP concentration (cAMP); this small molecule will bind to Crp (cAMP receptor protein) to form a complex that controls the expression of different genes, some involved in the synthesis of enzymes that allow the bacteria to make use of other carbon sources, others involved in morphological changes and motility. The complex cAMP-Crp also activates the enzyme Cya (adenylate cyclase), which contributes to produce cAMP from ATP, and represses the global regulator Fis, a protein which is available at high concentration during the exponential phase, and is responsible for the control of many other genes. The protein Fis also represses the complex cAMP-Crp and, among others, it controls two enzymes

= (u 1 , u 2 , u 3 ), v = (v 1 , v 2 , v 3
). The dashed lines represent the interconnection: i.e., the output of one system becomes the input of the other. Bacterial growth rate is internally computed as a function of the external nutrient sources (Glu), ribosomes (here represented by rrn i ), RNA polymerase (pol i ) or "bulk" proteins (which will be basically represented by crp). Growth rate is first translated into two qualitative levels, W 1 and W 2 , which signal downstream. The region under hatching represents the new variables and interactions added to the original model in [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF]. Several different forms for f g will be tested (see text). involved in DNA topology regulation: Gyrase AB (GyrAB) which induces negative supercoils in the DNA and Topoisomerase A which restores supercoiling to "normal" levels. Finally, Fis also stimulates the transcription of stable RNAs, a necessary condition for the production of ribosomes and hence necessary for bacterial growth.

The model developed in [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF] includes a constant external input named "Signal" that represents nutritional stress, that is, the presence ("Signal"= 0) or absence ("Signal"= 1) of carbon sources, while the variable rrn is simply an output, as it does not influence the other variables. Growth rate was summarized into the effect of the complex cAMP-Crp on the other variables, namely Fis, Crp, and Cya. Depending on the value of "Signal", the network reproduced two steady states corresponding to the stationary or exponential phases of E. coli, characterized in Table 1. The two states predicted by this model are consistent with experimental observations: in the exponential phase, Fis is present at high levels, as well as stable RNAs, and the cAMP receptor protein is not strongly present. The opposite happens in the stationary phase.

In our model, the interactions are reorganized in order to include the explicit effect of growth rate. It is known that the complex cAMP-Crp is growth dependent [START_REF] Berthoumieux | Shared control of gene expression in bacteria by transcription factors and global physiology of the cell[END_REF], so we replaced this complex by an equivalent expression that depends on Crp, Cya, and growth, now represented by the arrows u 1 , q r and q y in Fig. 5. The components inside the hatched region in Fig. 5 were not present in model [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF], and the rrn variable did not influence the system. The objective in this paper is thus to refine the effect of growth in the system, as described below in Section 3.2.

The model [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF] consists of a piecewise affine system on six variables, it was further studied in [START_REF] Grognard | Piecewise-linear models of genetic regulatory networks: theory and example[END_REF][START_REF] Tournier | Hierarchical analysis of piecewise affine models of gene regulatory networks[END_REF][START_REF] Corblin | A declarative constraint-based method for analyzing discrete genetic regulatory networks[END_REF] and has been written as an extended Boolean model in [START_REF] Chaves | Comparing Boolean and piecewise affine differential models for genetic networks[END_REF], using the procedure briefly described in Section 2. The first Boolean module is formed by the 8 variables corresponding to genes fis, gyr, and top, since fis is described by 4 Boolean variables and gyr, top by 2 each (see Fig. 5). The rules for this module are given in the Appendix.

Since each variable may have several discrete values, the Boolean models will use var i , i ∈ {1, 2, . . . , d} to denote the corresponding d Boolean variables (see Section 2.1) (similarly for the other variables). The discrete variable can be recovered simply by adding the Boolean variables:

var = d i=1 var i . (8)

The cellular growth module

To test the dependence of growth rate on some of the major model components, we will study a "closedloop system": that is, use the state of the system to construct a mathematical expression for bacterial growth rate and then feed it back to the system, by letting proteins Cya and Crp depend on it. Thus, cellular growth rate (represented by µ) now appears explicitly in the model, as an internal variable that depends dynamically on the state of the system at each instant (see Fig. 5). In agreement with the variables of the system, growth rate will have two positive discrete levels (translated to W 1 and W 2 , see equation ( 9) below). This also implies that the effect of growth on fis, crp and cya has to be updated relative to the original model [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF]. In Fig. 5, there are thus three links (respectively, u 1 , q r and q y ) which are not fixed for now, but for which several possible combinations will be tested, with a view to better understand growth signaling (see Section 4). The motivation for building this closed-loop system is to test the dynamical dependence of bacterial growth rate on the system's variables, a question which is still not well understood. Thus, for this example, an expression for growth rate will be considered valid if the refined system in Fig. 5 is able to reproduce the same results as the (more schematic) model [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF].

As indicated in Fig. 5, the second Boolean module will describe the expression of the genes encoding for crp, cya, rrn, and will further include pol, to represent the expression of RNA polymerase, the enzyme responsible for gene transcription (2 Boolean variables each). The presence of carbon sources will be represented by the external input Glu. Previously [START_REF] Carta | A simple model to control growth rate of synthetic E. coli during the exponential phase: model analysis and parameter estimation[END_REF], we have studied a mathematical expression for bacterial growth rate that is dependent only on RNA polymerase, for a simple 2-dimensional model. However, experimental data [START_REF] Shachrai | Cost of unneeded proteins in E. coli is reduced after several generations in exponential growth[END_REF] suggests that ribosomes play a major role, hence we wish to improve our results by analyzing models that consider different combinations of limiting factors, and checking their compatibility with known results.

The growth variable, µ, and its downstream signals will be given by:

µ = Glu and f g (rrn 1 , rrn 2 , pol 1 , pol 2 , crp 1 , crp 2 ); (9)
W = 2 -µ; W 1 = sign(W ); W 2 = max(0, W -1);
where sign(W ) = 1 if W > 0 and sign(W ) = 0 if W = 0 (by construction, sign(W ) is never negative).

The variables W 1 and W 2 correspond, respectively, to:

W 1 = 1 ⇔ µ ≤ 1, W 2 = 1 ⇔ µ = 0.
and satisfy W 1 ≥ W 2 . Following ( 5)-( 7) and ( 8), different expressions for the function f g will be tested, namely:

f r g = rrn 1 + rrn 2 , f p g = pol 1 + pol 2 , f b g = crp 1 + crp 2 , (10)
f rp g = min(rrn 1 + rrn 2 , pol 1 + pol 2 ) f rb g = min(rrn 1 + rrn 2 , crp 1 + crp 2 ),
where the protein Crp is used as a surrogate for the level of expression of "bulk" proteins. In addition, to describe how the growth rate affects the genetic machinery, two functions need to be chosen: these correspond to the arrows labeled q y and q r (see below), which will also be a function of W 1 and W 2 .

Several possible combinations will be tested and the final results compared to the original model.

System interconnection

The full discrete system will thus have 7 variables, V = (fis, gyr, top, crp, cya, rrn, pol) ′ , with discrete levels d 1 = 4, d j = 2 for j = 2, . . . , 7 and state space:

Ω d = {0, 1, . . . , 4} × {0, 1, 2} × . . . × {0, 1, 2}.
The extended Boolean model will have 16 variables. As described in Section 2.3, the interconnection of two input/output asynchronous Boolean networks such as systems ( 13) and ( 14), is obtained by setting u = h B (b) and v = h A (a). Most of the input/output functions are already fixed by model [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF]. There is a new interaction between the two modules, due to the effect of the growth rate in fis, which is represented by u 1 in Fig. 5:

u 1 ∈ {W 1 , W 2 }, u 2 = crp 1 or crp 2 , u 3 = cya 1 or cya 2 , v 1 = fis 1 , v 2 = fis 2 or fis 4 , v 3 = fis 3 .
The goal in this paper is the discrimination between different variants of the model in Fig. 5, in order to choose the mechanism that better represents bacterial response. The variants cover:

• models for growth rate: f r g , f b g , f rb g , and f rp g ; • interactions between growth signals and the genetic machinery response: q r , q y , and u 1 .

As remarked above (Section 3.1), the interactions q r , q y , and u 1 in some sense replace the effect of the complex cAMP-Crp on the system, by including an explicit dependence on growth rate. To evaluate the new rules we will consider that there are two signaling stages, corresponding to the response of Cya/cAMP (the initial steps in the case of nutritional stress) and of Fis (global regulator). The response of crp will be timed with one or the other:

q r = u 1 or q r = q y .
The following distinct combinations for q y , q r and u 1 will be tested:

(I) q y = W 1 , q r = W 1 , u 1 = W 1 (II) q y = W 1 , q r = W 1 , u 1 = W 2 (III) q y = W 2 , q r = W 2 , u 1 = W 1 (11)
(IV )

q y = W 1 , q r = W 2 , u 1 = W 2 (V ) q y = W 2 , q r = W 2 , u 1 = W 2 (V I) q y = W 2 , q r = W 1 , u 1 = W 1 4 Results
As discussed above (cf. Section 3), the goal is to recover the behavior of the system as described in Ropers et al [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF] (Table 1) but now with growth rate "actually computed" by the bacteria, for the system in closed loop form which uses the state of the system. Various combinations of interactions and growth rate functions were tested, with the results summarized in Table 2 and discussed below.

Table 2: The attractors for each combination of growth rate function and interactions u 1 , q r , q y . Attractors σ i satisfy rrn = pol = 0, while attractors α j , j ∈ {2, 4, 24, 48, 52, 72}, satisfy rrn ≥ 1 and pol ≥ 1 (full characterizations are given in Sections 4.2 and 4.3). The indexes i, j denote the number of distinct states contained in the attractor. All the attractors of G as are also attractors of G: either they satisfy Prop. 1 and/or other methods, as indicated. The highlighted row ( * * * ) represents the model variants which better reproduce Table 1 results (see Section 4.4).

Growth rate function Interactions [START_REF] Chaves | Predicting the asymptotic dynamics of large biological networks by interconnections of Boolean modules[END_REF] Attractors, G as (Glu=1) Attractors, G as (Glu=0) Stationary Exponential Stationary

f r g ,f p g ,f rp g ,f rb g I σ 4 [Prop. 1(i)] α 4 [Prop. 1(i)] σ 4 [Prop. 2] II σ 4 [Prop. 1(i)] α 72 [Prop. 3] σ 4 [Prop. 2] III,VI σ 4 [Prop. 1(i)] α 2 [Prop. 1(i)] σ 4 [Prop. 2] * * * IV,V σ 4 [Prop. 1(i)] α 24 [Prop. 1(iii)] σ 4 [Prop. 2] f b g I,III,VI σ 52 [Prop. 4(i)] α 52 [Prop. 4(ii)] σ 4 [Prop. 2] II σ 48 [Prop. 4(i)] α 48 [Prop. 4(ii)] σ 4 [Prop. 2] IV,V σ 24 [Prop. 1(iii)] α 24 [Prop. 1(iii)] σ 4 [Prop. 2]
As an indication of the computational costs, application of the method presented in Section 2.3 to compute the attractors for model f rp g , case IV, gave the following results: • there are eight constant-input asynchronous transition graphs for each system (G A,u , or G B,v );

• on these graphs there are a total of 22 semi-attractors for system Σ A and 20 for Σ B ;

• the total number of vertices in the asymptotic graph is thus 20 × 22 = 440;

• as remarked in Section 2.3 (and [START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF]), the number of vertices in G as can be further reduced by eliminating those which are known to have no incoming arrow. This leads to only 90 vertices;

• the computational cost of finding the attractors of the interconnected system Σ has therefore been reduced from analysis of a size 2 16 = 65536 to a size 90 matrix;

• one should nevertheless consider the cost of computing this size 90 matrix, which involves reachability calculations in the 2 × 8 individual asynchronous transition graphs (the full process was very fast here, taking between 30-60 seconds for each model variant).

General properties

Some immediate observations from the results are:

• a common point to all model variants is that, in the presence of nutrient (Glu=1), G as always has two attractors which are both attractors of G, by application of Prop. 1(i) or (iii), or other methods (see Prop. [START_REF] Sánchez | A logical analysis of the drosophila gap-gene system[END_REF][START_REF] Albert | The topology of the regulatory interactions predicts the expression pattern of the Drosophila segment polarity genes[END_REF].

• for all model variants, the first attractor (σ i , i ∈ {4, 24, 48, 52}) has rrn = pol = 0 and the second attractor (α j , j ∈ {2, 4, 24, 48, 52, 72}) rrn = pol = 1. The first may be said to represent stationary phase, while the second stands for exponential phase.

• also for all model variants, in the absence of nutrient (Glu=0), there is only one attractor, σ 4 ; this can be verified directly (see Prop. 2 below). The stationary phase attractor σ 4 has four states and is characterized by:

σ 4 : fis = 0, gyr ∈ {1, 2}, top = 0, crp = 2, cya ∈ {1, 2}, rrn = 0, pol = 0. (12)
coinciding with the stationary attractor of model [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF] (see Table 1) with the exception of gyr and cya which oscillate between 1 and 2 (instead of being fixed at 2).

• all models involving the ribosomes or RNA polymerase as a growth rate limiting factor exhibit the same stationary phase and a similar exponential phase attractors, depending only on the choice of feedback interactions.

As an example, for model variant f rp g , case IV, the basins of attraction for σ 4 and α 24 are disconnected. The stationary phase attractor is formed of a single vertex, while the exponential phase attractor is composed of 9 vertices, as shown in Fig. 6. However, all system B semi-attractors coincide (Prop. 1(iii) is satisfied): In practice, the attractor in Fig. 6 can be reduced to (either) one of the horizontal rows, with three components only. All concentrations are fixed, except for fis, gyr, and top which are allowed to oscillate in any given increasing or decreasing order, provided that fis ≥ 1 and top ≤ 1.

In the case where no carbon sources are present, it can be shown that all model variants become the same, and hence exhibit the same stationary phase attractor. This is essentially due to the direct effect of growth rate on the synthesis of RNA polymerase. Proof : In the case Glu=0, we immediately have the steady state values for rrn and pol:

µ = 0 ⇒ pol 1 = pol 2 = 0 ⇒ rrn 1 = rrn 2 = 0
For the interactions W , q y , q r , and u 1 it also follows that:

W 1 = W 2 = 1 ⇒ q y = q r = u 1 = 1.
Together with the rules in Appendix, this leads to:

cya 1 = crp 1 = 1 ⇒ u 2 = u 3 = 1 ⇒ fis i = hf i which simplifies to fis + 1 = f is 2 , fis + 2 = f is 3 , fis + 3 = f is 4 , fis + 4 = fis 3 and not fis 4 .
Thus, at steady state, the values for fis satisfy fis i = 0, for all i, which in turn imply that all the outputs of system A are zero: v i = 0 for all i. The remaining concentrations can now be easily established from the Boolean rules, so it follows that there is only one attractor and that it is σ 4 (12).

Growth Rate limited by ribosomes or RNA polymerase

For the model variants using f r g , f p g ,f rp g , or f rb g , the stationary phase attractor σ 4 is always the same (as described in Section 4.1). The exponential phase attractor, α j , j ∈ {2, 4, 24, 72}, depends on the wiring and has j states characterized by : Note that cases IV and V (α 24 ) are similar to the exponential phase attractor of [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF] (see Table 1) (the only difference is rrn now fixed at 1, which seems reasonable for the exponential phase). Cases I,III,VI (α 2 ,α 4 ) fail to reproduce the levels of fis during exponential phase (here they are fixed at zero). Case II (α 72 ) also exhibits oscillations in crp and cya, which are not observed in Table 1. This attractor does not fit into Proposition 1, but an alternative way to show that it is not a spurious attractor, is to note that the set of states with ribosomes, RNA polymerase and Fis all at discrete level 1 is invariant, so trajectories cannot leave this set; therefore, an attractor with such properties must exist, with the only possible candidate being α 72 .

Proposition 3 Assume that f g ∈ {f r g , f p g , f rp g , f rb g }. The set Q = {x ∈ Ω d : rrn = 1, pol = 1, fis ≥ 1} is invariant.
Proof : From the Boolean rules (see Appendix), it suffices to note that:

rrn = pol = 1 ⇒ µ = 1 ⇒ W 1 = 1, W 2 = 0 ⇒ q y = q r = 1, u 1 = 0 and also rrn = pol = 1, µ ≥ 1 ⇒ rrn = pol = 1.
And then:

u 1 = 0 ⇒ fis 1 = h01 ≡ 1 ⇒ fis ≥ 1.
Therefore, the set Q is invariant.

Growth Rate limited by bulk proteins

For the model variants using f b g , the exponential phase attractors are characterized as follows:

α 24 : fis ∈ {1, 2, 3, 4} gyr ∈ {0, 1, 2}, top = {0, 1}, crp = 1, cya = 2, rrn = pol = 1, α 48 : fis ∈ {1, 2, 3, 4} gyr ∈ {0, 1, 2}, top = {0, 1}, crp ∈ {1, 2}, cya = 2, rrn = pol = 1, α 52 : fis ∈ {0, 1, 2, 3, 4}, gyr ∈ {0, 1, 2}, top = {0, 1}, crp ∈ {1, 2}, cya = 2, rrn = pol = 1.
The stationary phase attractors are similar in all variables except that rrn = pol = 0. Comparison with Table 1 shows many differences with respect to model [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF]. These attractors are also true attractors of G, as shown by application of the following result. Proposition 4 Define the sets P 0 and P 1 :

P 0 = {x ∈ Ω d : rrn = 0, pol = 0}, P 1 = {x : rrn = 1, pol = 1}.
Then:

(i) The set P 0 is invariant independently of the function f g ;

(ii) The set P 1 is invariant if f g = f b g . Proof : Invariance of P 0 follows directly from the Boolean rules for rrn i and pol i . For P 1 , it suffices to note that the Boolean rules imply (see Appendix): cya ≥ 1 and crp ≥ 1 which imply µ ≥ 1. Then, rrn = pol = 1 and µ ≥ 1 imply rrn = pol = 1.

Model discrimination

Based on the observations above and comparison of Tables 1 and2, it seems clear that the growth rate function should depend on the ribosomes and/or RNA polymerase. With this model, at steady state there is no difference between a dependence on ribosomes or RNA polymerase, although the transient dynamics do depend differently on these two species (simulations in Section 4.5). (This may be due to a very simplified model for the transcription/translation steps which is, however, not our aim to study here.) The model variants corresponding to I,III,VI do not satisfy the properties of the exponential phase attractor and can thus be eliminated. The interconnection of type II has most of the correct properties, but it allows the concentration of crp and cya to oscillate, in contrast to the original model. The cases that better fit the original model are IV and V, whose asymptotic behavior is indistinguishable. This is consistent with the observations (summarized in Ropers et al. [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF]) that: immediately upon carbon starvation, or absence of carbon source, transcription of the gene cya is activated, which leads to production of the protein Cya. Phosphorylation of Cya leads to synthesis of cyclic AMP, which in turn will bind to Crp and produce a complex [cAMP-Crp]. This complex will then control a variety of genes which are directly involved in the adaptive response of E. coli to a deprivation of carbon. Among others, it activates crp, and inactivates cya and the global regulator fis.

Furthermore, to establish the transition to exponential phase, and guarantee the presence of global regulator fis, the wiring interactions should be as in cases II, IV, and V, which all satisfy u 1 = W 2 : in other words, since W 2 = 1 corresponds to µ = 0, fis is inhibited only at low growth rate, as also observed in [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF]. The levels of crp and cya are in agreement with those of Table 1 if crp is activated at high growth rate level and the inhibition effect on cya is not later than the activation of crp, i.e., q r = W 2 , q y ≥ q r (recall that W 1 ≥ W 2 ).

In conclusion, to develop a more detailed continuous model, bacterial growth rate should depend on the ribosomes. For simplicity, one may even consider the ribosomes to be the only variable influencing growth rate (besides the external input), because no differences were observed between models with f r g , f rp g , or f rb g .

Dynamical behavior

By virtue of Theorem 1 and Propositions 1 to 4, we know that any trajectory of the interconnected model will eventually reach either the exponential or stationary phase attractors, depending on the initial condition and path in the graph G. To illustrate possible dynamical behaviors, we can generate trajectories in G, by randomly choosing the variable to be updated at the next instant, according to the rules. Note that this simulation does not need the graph G to be constructed; but, on the other hand, such simulations cannot characterize the full behavior of the network. Hence the usefulness of the asymptotic graph G as , which can now be completed with some statistical results on initial conditions and attractors reached.

For the statistical analysis we choose interconnection model IV, and growth rate models f rp g , f r g , and f p g . It must be noted that the asymptotic graph "looses" some trajectories of the full interconnected system as, to construct G as , the system is assumed to evolve in one of the constant-input/constantoutput graphs {a}×G B,α or G A,β ×{b} until reaching an attractor. In simulations, however, the system is allowed to switch before reaching an attractor, meaning that the basins of attraction are not really disconnected as might be suggested by the asymptotic graph. Monte Carlo simulations of the full model (10 4 randomly generated trajectories) assume that all transitions in G are equally probable and show that the two attractors are reached with similar frequencies: for f rp g , a fraction of 0.59 (0.57 for f r g , or 0.58 for f p g ) trajectories converge to the exponential phase attractor.

Note that the invariance results in Propositions 3 and 4 already provide an idea of the basins of attractions, since they imply notably that initial conditions of the form rrn = pol = 0 (resp., rrn = pol = 1) lead immediately to the stationary (resp., exponential) phase attractor. To obtain more information on the distribution of the basins of attraction, we have further analyzed the probability that the system converges to either attractor given an initial condition with variable var i = ℓ (where var i runs over the sixteen Boolean variables of the system and ℓ ∈ {0, 1}). We found that the convergence to either attractor depends essentially on the initial concentrations of RNA polymerase and ribosomes, while all other concentrations play minor roles (in agreement with Propositions 3 and 4). An interesting observation is that, for all variables except the polymerase, and for any initial condition, the probability of converging to exponential phase is higher than to stationary phase. It is also evident that the absence of polymerase immediately prevents convergence to exponential phase. In addition, we observe that all trajectories converging to exponential phase need to start with an intermediate (or higher) level of RNA polymerase (pol ≥ 1). Table 3 summarizes the statistics obtained from the Monte Carlo simulations. Our studies lead to the conclusion that RNA polymerase and ribosomes are both crucial for bacterial Table 3: Initial conditions and attractor reached, for some model variants, with interconnection of type IV.

Initial conditions

Attractor reached f rp g , f r g , f p g pol = 0 stationary pol ≥ 1 and rrn = 0 either pol ≥ 1 and rrn ≥ 1 exponential growth, but exert their roles at different times: initially, the presence of RNA polymerase is necessary to grow and reach the exponential phase (otherwise, if RNA polymerase is absent at time zero, the bacteria enter the stationary phase even in the presence of carbon sources), while ribosomes can be absent; at later times, the presence of ribosomes is essential to guarantee the entry into exponential phase.

Conclusions

Several dynamic model variants for bacterial growth rate that consider limitation by availability of the proteins needed for cell division (RNA polymerase for transcription, ribosomes for translation, or other "bulk" proteins) were tested and compared to a well established model. The main goal was to analyze (qualitative) feasibility of the wiring network, as well as the logical coherence of each model variant. This was accomplished by using a Boolean version of the model for nutritional stress response in [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF], coupled with a basic cellular growth module. We can conclude that Boolean models provide a useful framework for analysis of a system's dynamical behavior, convenient for hypotheses testing and model discrimination. This framework presents several advantages from a computational point of view, as many tools and algorithms are available for the study and rigorous analysis of the networks. In particular, using the interconnection of two Boolean modules, it is possible to compute the attractors of a large network at a much lower cost than with classical graph theoretical tools. However, the drawbacks of this methodology include problems related to identifying the two (or more) Boolean modules as well as the corresponding inputs and outputs, which are not always obvious (see also [START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF]). As the number of modules and inputs increases, also the computational cost will increase and a balance must be found. This is a topic that should be further developed in future work.

A number of interesting points arise from our qualitative analysis. First, it was clear that limitation of growth rate by the ribosomes is needed in order to correctly reproduce the asymptotic modes, as well as transient dynamics, of the original model [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF]. Second, in the presence of nutrient, our closedloop model -where bacteria internally compute their growth rate, rather than responding to an already fixed signal-has the capacity for bistability (i.e., two asymptotic modes, representing exponential and stationary phases). Thus the closed-loop model also recovers the correct response to initial conditions: if both ribosomes and RNA polymerase concentration is very low, then the bacteria cannot grow even in the presence of nutrient. In the absence of nutrient, only the stationary phase attractor remains, as should be expected. Finally, by comparison to [START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF], we were able to discard most of the model variants and retain several properties necessary to reproduce the original model's attractors.

Since our main goal was essentially theoretical, we have not fully explored the directions for model improvement suggested by our analysis. For instance, a more detailed module for transcription/translation including other components besides ribosomes and RNA polymerase, or the modeling of the "bulk" proteins in a more precise way. To conclude, although discrete models are, of course, not appropriate for a detailed description of a system or to answer more specific questions, this analysis constitutes a very useful preliminary study of growth rate models. It provides many indications and clues for future work on constructing a more detailed, continuous model of the system.

The rules for the cellular growth module can be written as follows:

crp + 1 = 1;
crp + 2 = (not q r and crp 1 and not v 1 ) or (q r and crp 1 and not (v 2 or v 3 )); cya + 1 = 1; cya + 2 = (not q y and cya 1 ) or (q y and (hy1 or hy2)); where the auxiliary expressions are hy1 = cya 1 and (not crp 1 or not crp 2 ); hy2 = cya 1 and not cya 2 and crp 1 and crp 2 ;

Figure 1 : 1 = 2 = 1 = [v and not b 2 ) 2 = [v and b 1 and b 2 ]

 1121222 Figure 1: Example I: the asynchronous transition graphs that define the dynamics of the two systems (4). Example I. Consider the following bi-dimensional systems A and B, with n A = n B = 2 and p A = p B = 1: a + 1 = u and (a 1 and not a 2 ),

Figure 2 :

 2 Figure 2: Example I: the asymptotic graph of the interconnection of the two systems (4). The crossproducts inside shaded squares belong to an attractor. The cross-product inside a white circle represents a transient state that can be excluded from the computation.

Figure 3 :

 3 Figure 3: Example II: the asynchronous transition graphs of systems A and B, for each fixed input.

Figure 4 :

 4 Figure 4: Example II: the asymptotic graph of the interconnection of systems A and B defined in Fig. 3. States inside light shaded squares belong to some attractor; there are two attractors in this graph. States inside white circles represent known transient state, which can be discarded from the computation. States inside light shaded circles represent all other states.

Figure 5 :

 5 Figure 5: The interconnection of the fis global regulatory module (left rectangle) and a basic cellular growth model (right rectangle). Each module has three inputs and outputs: u= (u 1 , u 2 , u 3 ), v = (v 1 , v 2 , v 3). The dashed lines represent the interconnection: i.e., the output of one system becomes the input of the other. Bacterial growth rate is internally computed as a function of the external nutrient sources (Glu), ribosomes (here represented by rrn i ), RNA polymerase (pol i ) or "bulk" proteins (which will be basically represented by crp). Growth rate is first translated into two qualitative levels, W 1 and W 2 , which signal downstream. The region under hatching represents the new variables and interactions added to the original model in[START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF]. Several different forms for f g will be tested (see text).

  B12 54 = B 17 74 = B 19 84 = {10111010}, while the A system semi-attractors are characterized by the levels of fis, with either fis = 1,fis = 2 or f is ≥ 3:

Figure 6 :Proposition 2

 62 Figure 6: The exponential growth phase attractor (model f rp g , case IV). Since the B j α4 components are all equal, this attractor can be reduced to (either) one of the horizontal rows, with three components only (see text for more details).Proposition 2 Assume that Glu=0. Then, the asymptotic graph for all model variants exhibits only one attractor, σ 4 .

α 2 :

 2 fis = 0, gyr ∈ {1, 2}, top = 0, crp = 2, cya = 2, rrn = pol = 1, α 4 : fis = 0, gyr ∈ {1, 2}, top = 0, crp = 2, cya ∈ {1, 2}, rrn = pol = 1, α 24 : fis ∈ {1, 2, 3, 4} gyr ∈ {1, 2}, top = 0, crp = 1, cya = 2, rrn = pol = 1, α 72 : fis ∈ {1, 2, 3, 4}, gyr ∈ {0, 1, 2}, top = 0, crp ∈ {1, 2}, cya ∈ {1, 2}, rrn = pol = 1.

1 = pol 1 or rrn 2 ; rrn + 2 = pol 2 and rrn 1 and v 3 ; pol + 1 = 2 =

 1212 (sign(µ) and rrn 1 and pol 1 ) or pol 2 ; pol + sign(µ) and rrn 2 and pol 2 ;

Table 1 :

 1 The two E. coli modes reproduced by the model[START_REF] Ropers | Qualitative simulation of the carbon starvation response in Escherichia coli[END_REF]. If a variable has more than one value, this means that the asymptotic solution is oscillatory among those values.

	"Signal"	fis	gyr top crp cya rrn Phase
	0	1,2,3,4 1,2	0	1	2	0,1 Exponential
	1	0	2	0	2	2	0	Stationary

Acknowledgments

We are especially grateful to Laurent Tournier for many discussions and for providing part of the Matlab codes used here to analyze the asynchronous transitions graphs (specifically, the decomposition into strongly connected components and subsequent hierarchical organization). We also thank Jean-Luc Gouzé and our reviewers for many useful suggestions that helped improve the paper.

This work was supported in part by projects GeMCo (ANR 2010 BLAN0201-01) and ColAge (Inria-INSERM large scale initiative action).

A Boolean rules of the two E. coli modules

The Boolean model for the Fis module is defined by a set of rules which use some auxiliary expressions of the form h-given below: h04 = hf 4n and gyr 1 and not top 2 ; h11 = ((u 2 or u 3 ) and hf 2 ) or ((not u 2 or not u 3 ) and h01); h12 = ((u 2 or u 3 ) and hf 3 ) or ((not u 2 or not u 3 ) and h02); h13 = ((u 2 or u 3 ) and hf 4 ) or ((not u 2 or not u 3 ) and h03); h14 = (not u 2 or not u 3 ) and h04;


LIFE CYCLE ASSESSMENT: GAS EMISSIONS OF GREENHOUSE MATERIAL FLOWS FOR TWO TYPES OF CONSTRUCTION

Keywords: ALC, carbon footprint, systemic approach, materials, material flows. 1

, allows to estimate the environmental impacts of a product, from the extraction of raw materials to the end of life. In the building industry, this method enables to provide quantitative indicators of environmental impacts of the products and materials used.

This study provides a new tool for environmental evaluation on the carbon footprint that applies to buildings of any kind. Five types of material are identified at a first level of description: structural work products (PGO), Second work products (PSO), equipment products (PE), activity products (PA) and current consumption products (PCC). Several descriptive and analytic parameters are associated with these materials.

In this study we intend to assess the environmental impacts of the buildings of a holiday cottage, based on two methods of construction (mineral materials and organic type).

We will illustrate the methodological aspects related to the realization of this tool, based on a systemic and typological approach.

INTRODUCTION

Throughout its life cycle a building has a significant environmental impact (23% of CO 2 emissions and 75 % of the total mass of waste [START_REF] Durable | Observation et statistiques -Construction et bâtiments[END_REF]). To address this concern, a practice known as eco-design has developed. It is the integration of environmental aspects into product design and product development [START_REF] Standardization | [END_REF].

This concept proposes to take into account a number of principles, methods and technical to reduce environmental impacts at each stage of the implementation of the construction.

In the context described above, the aim of this article is directly related to the proposal of an environmental impact assessment methodology, applied to buildings in an analytical and global approach.

COMPARATIVE ANALYSIS OF LIFE CYCLE OF BUILDING MATERIALS

The comparative study is carried out on the materials used for two types of construction (mineral and organic). The construction of the inorganic type is essentially composed of stone. The building is made of concrete block, mortar and steel. As regards the construction of the organic type, the frame and the structure are made of wood.

The different treatments of products such as the steps of manufacturing, transportation, their implementation and the end of service life are examined and analyzed. Their study allows to determine the importance of each type of material, and propose selection criteria based on environmental footprint.

The concept of the tool is a systemic and typological approach, which makes it possible to assess the problems and developments of a system. Several scales are identified and analyzed: The idea is to identify, target and analyze all entities of the largest to the smallest scale. It is from this method is that it establishes the acquisition of data and calculations of carbon footprint.

The data acquisition is done in the form of a table as shown in Figure 2. This is the main draw for data acquisition products by type level. From this table, we input the product, unit, quantity and country of origin. The scope of application covers the manufacturing of products, transportation, packaging, and end of life as show below. However, given the difficulties in collecting data, we do not consider the maintenance and waste treatment, but make the assumption of a scenario of end of life for each material.

During data acquisition, the emission factors must be informed. For manufacturing carbon footprint, it is calculated either from the emission factor of the product (if it is specified) or from emission factors of the materials that compose it. The same applies to the transportation carbon footprint: the calculation is based on the emission factors of different modes of transport of the product (air, sea, river, road and rail). The reason why one should consider the types of products is that they facilitate the interpretation of the results for the building as well as the functional spaces. They make it possible to visualize the environmental footprints by type of product, and to provide assistance on the choice of materials.

ASSESSMENT OF THE ENVIRONMENTAL FOOTPRINT

The carbon footprint simulation was performed by the creation of a new Excel tool and the use databases INIES [START_REF] Inies | Base nationale française de référence sur les impacts environnementaux et sanitaires des produits, équipements et services pour l'évaluation de la performance des ouvrages[END_REF] and Carbon data [START_REF] Ademe | Base Carbone[END_REF]. The environmental impacts were calculated according to the method described in ISO 14040. This method ISO 14 040 is used to connect the inventory data for the environmental damage they create. The results obtained in this study show that some building materials contribute very significantly to the generation of the environmental impact of a construction.

It is interesting to note here (Figure 3) that the carbon emissions are significantly variable depending on the type of product considered and the type of construction.

Figure 3: Comparative analysis of emissions of carbon for two types of construction

It is also interesting to observe the influence of the sub-types of level (Figure 4). The carbon emissions of the organic building are considerably reduced compared to the mineral construction. The largest variances noted are related to the superstructure and bulkheads. Indeed, as for the evaluation of the carbon footprint of the structural work and second work, only the materials that relate to the frame and the structure of the building will have a significantly different value. Emissions of other product types are similar (such as electricity and plumbing), since they are the same regardless of the construction materials considered. In this sense, it is important to consider all types of product level 0 (Structural work and second work) and level 1 (Infrastructure, superstructure...) to be able to target the environmental impact of buildings.

CONCLUSIONS

Based on a selection criterion taking into account the environmental impact, a significant reduction of the carbon pollutants is possible. According to this study, the organic building type construction would have less impact than mineral construction with a gap of 373 t CO 2 eq. (Either -295 t CO 2 eq. for organic construction and 77 t CO 2 eq. for mineral construction). More precisely it is when considering the emissions of materials for the structure and the frame, that one can note a difference between the two types of construction.

In this sense, the study not only highlights the influence of the choice of materials, but mainly showed the point of the LCA approach to assess the environmental impact of a building. Indeed, this methodology allows both to have a synthetic and analytical approach.

SYNTHESES

Le bâtiment tout au long de son cycle de vie, contribue aux émissions de gaz à effet de serre. Face à l'enjeu d'anticiper les impacts environnementaux des activités humaines, une pratique s'est développée : l'éco-conception. Ce concept propose de prendre en compte un certain nombre de principes, de méthodes, et de techniques pour limiter les impacts futurs des produits à chaque étape de la mise en oeuvre de la construction.

La méthode d'analyse du cycle de vie (ACV) décrite par la série des normes ISO 14040 est un outil répandu pour l'évaluation de l'impact environnemental d'un produit. Cette méthodologie peut s'appliquer aux procédés et aux services.

Dans le cadre de cet article, le but est de comparer les divers matériaux entrant dans la composition de deux types de construction, pour utiliser comme critère de choix leur empreinte environnementale.

La présente étude est consacrée à la présentation d'une méthodologie d'ACV selon une approche systémique et typologique, appliquée à deux types de construction : minéral (béton) et organique (bois). Cinq typologies de matériau sont identifiées à un premier niveau de description : Produits Gros OEuvre (PGO), Produits Second OEuvre (PSO), Produits d'Equipement (PE), Produits d'activité (PA) et Produits de Consommation Courante (PCC). Les simulations pour cet article, ont été effectuées pour les émissions de gaz à effet de serre (CO2) des Produits Gros OEuvre (PGO) et des Produits Second OEuvre (PSO).

Pour considérer tous ces paramètres, l'analyse a été réalisée à partir de la création d'un outil sous Excel et en s'appuyant sur les bases de données existantes (INIES et Carbone). Le champ d'application couvre la fabrication des produits, le transport, l'emballage et la fin de vie. Le bilan carbone pour chaque étape de traitement du produit est calculé à partir des facteurs d'émission, de l'unité fonctionnelle et des quantités.

Les résultats obtenus montrent d'une part quel est la construction la plus impactant et d'autre part ciblent les typologies de produit majoritairement contributrices. En effet, d'après les simulations effectuées, la construction de type organique aurait moins d'impact qu'une construction minérale, avec un écart de 373 t CO2 eq. (Soit -295 t CO2 eq. pour la construction organique et 77 t CO2 eq. pour la construction minérale). Plus précisément c'est en considérant les émissions des matériaux pour la structure et l'ossature du bâti, qu'une différence est notable entre les deux types de construction.

En ce sens, l'étude a permis non seulement de mettre en évidence l'influence des choix des matériaux, mais a montré surtout l'intérêt de la démarche de l'ACV pour l'évaluation de l'empreinte environnementale d'une construction. En effet, cette méthodologie permet à la fois d'avoir une approche synthétique et analytique.

Par ailleurs, d'autres éléments d'analyse seront ajoutés à l'outil. Pour considérer l'ensemble des émissions carbones de l'espace bâti, les Produits d'Equipement, les Produits d'activités et les Produits de Consommation Courante, seront traités. Ce même outil permettra à terme d'effectuer d'autres analyses d'impacts environnementaux tels que l'utilisation des ressources ou encore l'acidification atmosphérique.

Figure 1 :

 1 Figure 1: Diagram of systemic and typological approch

Figure 2 :

 2 Figure 2: Architecture of the scope and data acquisition

Figure 4 :

 4 Figure 4: Comparative analysis of carbon emissions (t CO2 eq.) of big work and second work for two types of construction

PGO


LIFE CYCLE ASSESSMENT : GAS EMISSIONS OF GREENHOUSE MATERIAL FLOWS FOR TWO TYPES OF CONSTRUCTION Application of a holiday cottage!

BACKGROUND

PURPOSE

The life cycle assessment (LCA) of the series of ISO 14040 standards [1], allows to estimate the environmental impacts of a product, from the extraction of raw materials to the end of life. In the building industry, this method enables to provide quantitative indicators of environmental impacts of the products and materials used. This concept proposes to take into account a number of principles, methods and technical to reduce environmental impacts at each stage of the implementation of the construction. In the context described above, the aim of this article is directly related to the proposal of an environmental impact assessment methodology, applied to buildings in an analytical and global approach. . This study provides a new tool for environmental evaluation on the carbon footprint that applies to buildings of any kind. Five types of material are identified at a first level of description: structural work products (PGO), second work products (PSO), equipment products (PE), activity products (PA) and current consumption products (PCC). Several descriptive and analytic parameters are associated with these materials. In this study, we intend to assess the environmental impacts of the buildings of a holiday cottage, based on two methods of construction (mineral materials and organic type). We will illustrate the methodological aspects related to the realization of this tool, based on a systemic and typological approach.

The concept of the tool is a systemic and typological approach which makes it possible to assess the problems and developments of a system. Several scales are identified and analyzed :

Your caption can go here. Based on a selection criterion taking into account the environmental impact, a significant reduction of the carbon pollutants is possible. According to this study, the organic building type construction would have less impact than mineral construction with a gap of 373 t CO2 eq. (Either -295 t CO2 eq. for organic construction and 77 t CO2 eq. for mineral construction). More precisely it is when considering the emissions of materials for the structure and the frame, that one can note a difference between the two types of construction.

METHODS

In this sense, the study not only highlights the influence of the choice of materials, but mainly showed the point of the LCA approach to assess the environmental impact of a building. Indeed, this methodology allows both to have a synthetic and analytical approach.

Furthermore, other test elements will be added to the tool. To consider all the carbon emissions of the built space, equipment products (PE), Activity products (PA) and consumer current products (PCC), will be addressed. This same tool will eventually perform other analyzes of environmental impacts such as the resource use or the air acidification.

RESULTS

CONCLUSION REFERENCES

The idea is to identify, target and analyze all entities of the largest to the smallest scale. It is from this method is that establishes the acquisition of data and calculations of carbon footprint, for two types of construction applied to a case study.

The scope of application covers the steps of life cycle assessment. However, given the difficulties in collecting data, we do not consider the maintenance and waste treatment, but make the assumption of a scenario of end of life for each material. It is also interesting to observe the influence of the sub-types of level. The carbon emissions of the organic building are considerably reduced compared to the mineral construction. The largest variances noted are related to the superstructure and bulkheads. Indeed, as for the evaluation of the carbon footprint of the structural work and second work, only the materials that relate to the frame and the structure of the building will have a significantly different value. Emissions of other product types are similar (such as electricity and plumbing), since they are the same regardless of the construction materials considered. In this sense, it is important to consider all types of product level 0 (Structural work and second work) and level 1 (Infrastructure, superstructure...) to be able to target the environmental impact of buildings.

  [1] N. E. I. 14040, Management environnemental. Analyse du cycle de vie, Principes et cadre, 2006. [2] INIES, "Base nationale française de référence sur les impacts environnementaux et sanitaires des produits, équipements et services pour l'évaluation de la performance des ouvrages," 2013. [Online]. Available: www.base-inies.fr/. [Accessed avril 2014]. [3] ADEME, "Base Carbone," [Online]. Available: www.basecarbone.fr. [Accessed avril 2014].

  emissions of carbon for two types of construction Comparative analysis of carbon emissions (t CO2 eq.) of big work and second work for two types of construction The carbon footprint simulation was performed by the creation of a new Excel tool and the use of databases INIES [2] and Carbon data [3]. It is interesting to note here that the carbon emissions are significantly variable depending on the type of product considered and the type of construction.

Acknowledgements

This research has been supported by :