Olivier Michel
email: olivier.michel@univ-paris12.fr

Antoine Spicher
email: antoine.spicher@univ-paris12.fr

Jean-Louis Giavitto
email: giavitto@ibisc.univ-evry.fr

Rule-Based Programming for Integrative Biological Modeling Application to the Modeling of the λ Phage Genetic Switch

Keywords:

HAL is

In all these models, computations are specified by a set of local rules that describes how neighboring data may interact to build some new data. Their specificities lie in how the neighborhood relationship is defined: abstract chemistry with (multi-)sets, P systems with a spatial organization formalized as Venn diagrams without intersection (corresponding to nested compartments), cellular automata with a regular neighborhood, L-systems with tree-like structures, etc. In other words, these models may be seen as specific forms of local application of rewriting rules on different kinds of spatially organized data [START_REF] Giavitto | Data structure as topological spaces[END_REF].

In the following paragraphs, we detail how a rewriting system can be used to model some natural system, starting from their usual point of view in computer science to their use for natural systems.

Rewriting and Computation

A rewriting system (8) is a formal system used as a mean to compute by substituting a part of an entity by another. A rewriting system is defined by a set of rules R. Each rule R i ∈ R: α → β is composed of a left-hand side α (lhs for short), specifying the part to be substituted, and a right-hand side β (rhs for short), specifying the new part replacing α.

Computing by local substitutions is common in computer science: it has been developed in a lot of contexts and applied on various kinds of entities. For instance in multiset rewriting, the lhs is a submultiset which is replaced by the multiset corresponding to the rhs; in term rewriting, the lhs is a subtree which is replaced by the tree corresponding to the rhs; in graph rewriting, the lhs is a subgraph which is replaced by the graph corresponding to the rhs; etc.

Let the variables e, e , . . . denote the considered entities. A rewriting rule defines a binary relationship between these entities. We write e ⇒ Ri e to denote the local transformation of e into e by applying the rule R i : e is defined as e where one of its subpart α is replaced by β. If no α exists in e, the relation is not defined on e.

We write e ⇒ R e to denote that e is transformed into e by applying the rules in R. We need here to be more specific about which rules are applied ( This sequence is called a derivation of e w.r.t. the set of rules R and the strategy St. This transformation of e into e can be seen as the result of the computation specified by R where the derivation corresponds to the sequence of the intermediate results.

Rule-Based Modeling of Natural Systems

The framework proposed by rewriting systems suits well the formal description of natural phenomena. As a matter of fact, many authors (9; 10; 11; 12) have recognized the importance of term rewriting in that field. The state of a biological system is represented by a term t of the form

t 1 + t 2 + • • • + t n
where a subterm t i represents a biological entity composing the system or a message (signal, information, process, etc.) sent to an entity. Following this point of view, the dynamics of the biological system is captured by the rewriting steps of the term. For example, rules of the form e + m → e + m can be used to model the reaction of the entity e when it receives the message m: the lhs corresponds to the message m and to the entity e which m is sent to. The + operator denotes that e and m are somehow neighbor in the system. In the rhs, e specifies the new state of the entity and a message m may (possibly) be sent. Other kinds of local evolutions may be specified in the same way. The direct interaction between entities can be expressed by rules of the form e 1 + e 2 → e 1 + e 2 and the creation of a new entity by a rule of the form e → e + e etc. Intuitively, we can deduce that, given the state of a system as a term t, a derivation of t is understood as the description of the evolution of the system along time; we speak of the trajectory of the system.

Table 1 Dictionary establishing a link between the notions required for the modeling and simulation of natural systems with their counterparts in the field of rewriting systems.

Natural System Rewriting System

Model

Spatial organization A data structure T T ∈ {term, set, graph, . . . } State An element e ∈ T Local evolution laws Set of rules R Interaction =⇒ result A rule α → β ∈ R Simulation Local and atomic evolutions Relation ⇒ R i with R i ∈ R Time steps Relation ⇒ R,St Trajectories Relation ⇒ * R,St
More generally, the modeling and the simulation of a natural system using a rewriting system may be summarized by the analogy given in Table 1. In fact, choosing rewriting systems as a formalism to describe dynamical systems implies to consider the handling of discrete time and local space.

Local Specification of Space. For a particular natural system to be modeled, according to Table 1, it is required to choose a particular type T of data structure that describes its local spatial organization. The state of the system is given by a value of type T . Once the state is defined, a rule R i must be defined, for each interacting subpart, to specify the set R of local evolution laws. Such rules can be understood as following:

-the lhs represents a subsystem whose elements are in interaction, -the rhs corresponds to the computation of the local result of this interaction.

In the previous example, the + operator connects entities and signals as they are spatially and/or functionally organized. It denotes at the same time the interacting subparts of the system as well as it specifies the way that the subsystems are organized. Terms built using arbitrary operators can be used to capture hierarchical organizations. However the + operator is (usually) associative and commutative. This implies that the tree structure of an arbitrary term can be reorganized to order arbitrarily the leaves of the tree structure (the leaves of the tree are the constants of the expression). Using associative-commutative operators we obtain multiset. In a multiset, all elements are considered as being all neighbors which corresponds to the metaphor of a well mixed chemical solution where Brownian motion ensures that each element will encounter all the others. Discrete Time. The handling of time is a key point in the modeling of a natural system. The model of time that is naturally induced by the use of rewriting systems, is a discrete time: the application of a rule corresponds to an atomic modification of the state of the system and to the progression of time.

As shown in Table 1, the trajectories of the systems in its phase space are clearly captured by the definition of the relations ⇒ R,St and ⇒ * R,St . Intuitively, they correspond to event-based simulations of the model as they represent the iterative applications of the local evolution laws specified within R with respect to a rule application strategy St.

The MGS project

The MGS project1 proposes to investigate the previously referred computation models by developing a domain-specific programming language, also called MGS, where they may be expressed in a uniform and generic setting. In MGS, the state of a dynamical system is specified using an original and generic data structure, called topological collection [START_REF] Giavitto | Data structure as topological spaces[END_REF], based on the topological relations between the interacting subparts of the system. Furthermore, the specification of the evolution law is simplified by the definition of a case-based function, called transformation, based on the rewriting of topological collections.

Topological Collection. One of the key features of the MGS language is its ability to describe and manipulate entities structured by an abstract topology. Topological collections provide a unified view of the notion of data structure [START_REF] Giavitto | Data structure as topological spaces[END_REF] seen as an aggregate of elements organized with a neighborhood relationship. The structure of the defined topological space allows to specify the organization of the data structure.

The formalization of topological collections has been thoroughly studied in previous work of the authors (13; 14). Without getting too much into details, we give a few technical details on topological collections. Its formalization relies on the notion of abstract cell complex [START_REF] Munkres | Elements of Algebraic Topology[END_REF] defined in algebraic topology allowing the description of the topology of the collection. A cellular complex makes it possible to have a discrete representation of a topological space through a set of topological cells (an abstraction of elementary spaces characterized by their dimensions) connected to each other through their incidence relationship. This relation is based upon the notion of boundary. This structure is then decorated with values leading to the notion of topological chain [START_REF] Munkres | Elements of Algebraic Topology[END_REF]. This final notion of topological chain allows the association of elements of an abelian group to cells of a complex resulting in a rich algebraic structure.

In the context of the modeling of natural systems, topological collections allow an intuitive representation of the states of the system: the elements of the collection are the atomic components, and the interaction graph between components specializes the topological structure of the collection.

Many topologies are available in MGS. In this article, we focus on multiset collections (corresponding to a complete neighborhood) in sections 3 and 4, GBF collections (corresponding to a regular topology) in section 5.1, and Delaunay collections (whose topology relies on the computation of a Delaunay triangulation) in section 5.2.

Transformation. Topological collections represent an adequate medium to extend rewriting systems as previously presented. Indeed, the neighborhood relationship gives a local point of view of the organization of the elements. The transformations extend the notion of rewriting systems to topological collections, and are used to specify the evolution function of the modeled system as show on Table 1.

In MGS, the specification of a transformation T :

trans T = { σ => f (σ, . . . );
... } corresponds to the definition of a set of rules, where σ is a pattern, matching for a subcollection, and f (σ, . . . ) is an expression that evaluates a new subcollection that will be inserted in place of the matched one. The application of transformation T on a topological collection e using a strategy St, written T[strategy = St](e), consists in computing a collection e such that e ⇒ T,St e .

In the current implementation of MGS, all available strategies are built-in. In the following, we will use two of the many rule application strategies provided by MGS: the Gillespie strategy based on the stochastic simulation algorithm proposed by Gillespie to simulate chemical reactions (in sections 3 and 4), and the maximal-parallel strategy (in section 5) widely used in the context of L-systems and P systems.

Motivations and Organization of the Article

MGS provides features allowing the use of rewriting techniques on a large number of different data structures. This point of view is very effective in the modeling of natural systems as it allows to use a different type of organization in the lhs and in the rhs of a rule. Consequently, rules may specify modifications of the system's structure. Thus, this allows the modeling of a class of dynamical systems particularly complex to describe and simulate, dynamical systems with a dynamical structure (16; 17; 18).

In this article, we claim that the generic approach of MGS allows an easy way to express models of biological systems. Previous publications by the authors were devoted to the formalization of the theoretical foundations of MGS, and to the use of MGS for the modeling and simulation of biological processes at a specific level of description.

Here, we focus on the suitability of MGS for the description of the same biological process, from the molecular level to the level of spatially distributed populations. The article is organized as follows: in section 2, we describe the biological process of the regulation of the genetic switch of the λ phage. This process will be our running example throughout this article.

We propose in section 3 a first model of the switch based on a biochemical description of the regulation. For that purpose, we show how Gillespie's stochastic simulation algorithm can be seen a rewriting strategy in MGS. This first model only consider a single cell.

A first model of a population of cells is proposed in section 4. In this section, the population of cell is considered as being homogeneously distributed in space and therefore represented by a multiset topological collection in MGS.

We still consider a population of cells and complexify our model in section 5 by considering a population in an homogeneous and heterogeneous spatial distribution. For that purpose, we propose two models, one relying on a cellular automaton with a Von Neumann neighborhood and one with an individual-based model on a neighborhood defined by a Delaunay triangulation of the cells in space. Both models are naturally implemented in MGS.

Finally, we conclude in section 6 by emphasizing the importance of rule-based languages and formalisms in systems biology and we stress the importance of spatial relationships in biological processes.

The Paradigmatic Example of the λ Phage Switch

In the forthcoming sections, we propose to illustrate the expressiveness brought by rule-based programming for the description of various models of the same biological process: the regulation of the switch of the λ phage. We first start in this section by briefly describing the behavior of the phage.

The Genetic Switch

The λ phage ( 19) is a virus that infects the cells of the bacterium Escherichia coli. It is a phage that has two possible outcomes:

1. replication and lytic phase where the virus dissolves and destroys the host cell, releasing about 100 virions; 2. integration of its DNA in the DNA of the bacterium, and start of a lysogenic phase.

In the lysogenic phase, the virus will silently replicate at each cell division. Moreover, a lysogeny produces an immunity towards further phage infection, by protecting the bacterium from the destruction during a possible new infection by a phage. Under certain conditions (exposure to U.V. for example) a lytic phase can be induced : the viral DNA is released from the bacterial DNA and starts a normal replication and a lysis.

Based on the local conditions, it is one of the two possible phases that is chosen, the decision being under the control of a small region of the phage genome (a hundred base pairs, comparatively to the 48502 base pairs of the genome of the bacterium) and of two genes (cI and cro) and two promoters This regulatory region is called the genetic switch.

Regulation of the Genetic Switch

The region of the DNA (see Figure 1) of the λ phage is composed of two genes cI and cro coding respectively for the proteins CI and CRO. During the transcription, the RNA polymerase binds to the promoters of these genes (Prm and Pr respectively) to synthesize the mRNA that is then translated into monomer proteins CI and CRO.

These monomers dimerize into CI 2 and CRO 2 which can bind to the operators. Operators OR1, OR2 and OR3 overlap the promoters (see Figure 1). The absence or presence of dimers bound to the promoters, eases or hinders the binding of the RNA polymerase, thus regulating the expression of genes cI and cro. Binding of the dimers to the promoters follows certain rules of affinity:

-CI 2 binds first to OR1, then to OR2 and finally to OR3; -CRO 2 binds first to OR3, then to OR2 and finally to OR1; -If CI 2 is bound to OR1, it facilitates the binding of another dimer CI 2 to OR2 and consequently the binding of the RNA polymerase on the Prm promoter. This property does not hold for CRO 2 ; -The RNA polymerase can be bound to the promoter Pr but CI 2 has to be bound to OR2 so that it can be bound to Prm.

All these rules together lead to a complex behavior and the genetic regulation of the transcription.

Simply stated, each of the two proteins CI and CRO will inhibit each others and increase its own synthesis (they self-inhibit themselves when their expression becomes too high). According to the local conditions, one of the two proteins will take over the other one (in terms of concentration) and the system will enter either a lytic (CRO "wins") or a lysogeny phase (CI "wins").

Biochemical Modeling of a Single Bacterium

Our first modeling of the biological process of the phage genetic switch, relies on a translation of the regulation principles into biochemical reactions. It is a straight translation of the previously described process where the spatial distribution of the biological entities is not taken into account.

Chemical Reactions Based Model

The stochastic simulation of biochemical systems becomes useful when the number of molecules and/or the time interval is very small. It is the case of the regulation of the λ phage where only around ten molecules CI 2 are enough to enter the lysogenic phase [START_REF] Lou | A quantitative study of λ-phage switch and its components[END_REF]. Figure 2 describes the model used for the chemical modeling. It follows the description given above. For the sake of simplicity, only a subset of all the dimers-operator bond combinations is considered. It reduces from 27 possible states (3 3 corresponding to all the combinations of no binding, CI 2 and CRO 2 bound on each operator) to only 7 states of the DNA. Actually, we have focused on the subset of the bonds that are most-likely to happen. CI 2 binding on OR2, with OR1 free, could have been taken into account; but this state is rare. We also consider that cI is expressed when OR2 is occupied by CI 2 and OR3 is free. The cro gene is expressed when both OR1 and OR2 are free. The set of interactions (dimerization, binding of proteins on the DNA and gene expression) are treated as the following biochemical equations:

CI C0 G GGGGGG A . (1) 2 CI C12 GGGGGGB FGGGGGG C21 CI 2 (2)
D + CI 2 C 1 GGGGGGB FGGGGGG C -1 D 1 (3)
D 1 + CI 2 C 2 GGGGGGB FGGGGGG C -2 D 2 (4)
D 2 + CI 2 C 3 GGGGGGB FGGGGGG C -3 D 3 (5) D 2 + P Ct G GGGGGG A D 2 + P + CI (6) CRO C 0 G GGGGGG A . (7) 2 CRO C 12 GGGGGGB FGGGGGG C 21 CRO 2 (8) D + CRO 2 C 3 GGGGGGB FGGGGGG C -3 D 3 (9)
D 3 + CRO 2 C 2 GGGGGGB FGGGGGG C -2 D 2 (10)
D 2 + CRO 2 C 1 GGGGGGB FGGGGGG C -1 D 1 (11) D + P C t G GGGGGG A D 2 + P + CRO (12) D 3 + P C t G GGGGGG A D 2 + P + CRO (13)
Equations ( 1) and ( 7) describe the natural degradation of the CI and CRO monomers. Equations ( 3)-( 5) and ( 9)-( 11) express the bindings of the dimers on the operators; the different states of the DNA are represented: constant D corresponds to the DNA with no bonds, D 1 , D 2 and D 3 to the DNA with 1, 2 or 3 dimers CI 2 bound, and D 3 , D 2 and D 1 to DNA with 1, 2 or 3 dimers CRO 2 bound -see Figure 2. The gene expression is given by reactions ( 6), ( 12) and ( 13), where P stands for the RNA polymerase. Each reaction is parameterized with a stochastic constant, C i for the reactions involving CI, and C i for CRO.

Rule-Based Programming of Chemical Reactions

A usual abstraction in the simulation of biochemical systems consists in considering the system, here the bacterium, as an homogeneous chemical solution where the reactions of the model are taking place. As the solution is considered well mixed, it can be represented by a multiset, that is a topological collection where any element may interact with all the other [START_REF] Fisher | Spatio-logical processes in intracellular signalling[END_REF]. This point of view is the starting of abstract chemical computation models like CHAM and Γ . However, the simulation of "real" chemical reactions requires more than only multiset rewriting to take into account the kinetics, corresponding to each reaction constant C i and C i . D.T. Gillespie has proposed in [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF] an algorithm for producing the trajectories of a chemical system by choosing the next reaction and the elapsed time since last reaction occurred. Let µ be a chemical reaction, the probability that µ takes place during an infinitesimal time step is proportional to:

-c µ , the stochastic reaction constant2 of reaction µ; -hµ, the number of distinct molecular combinations that can activate reaction µ; -dτ , the length of the time interval.

Gillespie proved that the probability P (τ, µ)dτ that the next reaction will be of type µ and will occur in the time interval (t + τ, t + τ + dτ ) is:

P (τ, µ)dτ = a µ e -a0τ dτ
where a µ = c µ h µ is called the propensity of reaction µ, and a 0 = P ν a ν is the combined propensity of all reactions.

This probability leads to the first straightforward Gillespie's algorithm called first reaction method. It consists in choosing an elapsed time τ for each reaction µ according to the probability P (τ, µ). The reaction with the lowest elapsed time is selected and applied on the system making its state evolve. A new probability distribution is then computed for this new state and the process is iterated.

Stochastic Simulation in MGS

We now present how the previous chemical model can be implemented in MGS. We need to represent the state of the dynamical system and its evolution function.

State of the System. As said above, the state of the bacterium is represented by a multiset of values. The considered molecules are abstractly represented using MGS symbols that are back-quoted identifiers. For example, the MGS symbol 'CRO2 corresponds to a dimer CRO 2 .

The initial state consists in a single molecule of the phage's DNA with some copies of RNA polymerase. As it has been noticed in vivo, the probability for CI to gain over CRO is low 3 . In order to have the simulation evolve in favor of CI, we put three copies of the CI protein to the initial state: An application of the transformation using this strategy consists in: 3 An even more realistic behavior is obtained by considering models involving gene cII and cIII [START_REF] Arkin | Stochastic kinetic analysis of developmental pathway bifurcation in phage λ-infected Escherichia coli cells[END_REF].

1. Computing the propensity a µ of each rule r µ by using the user-defined stochastic constant cµ (taken from the C constant of the rule) and by evaluating the number of combinations h µ (that is the number of subcollections matched by the lhs of the rule); 2. Computing for each rule r µ , the value of the elapsed time τ µ w.r.t. the definition of P (τ, µ)dτ ; 3. Applying once the rule with the lowest elapsed time.

Rule r µ with the smallest τ µ value is chosen and fired one time on the collection argument of the transformation, after instantiation of the lhs pattern.

A Simulation Example. During the iterations of the application of the Phage transformation using the Gillespie's strategy, MGS system variable tau is incremented using the elapsed times τµ of the chosen reactions. The value of tau is available anywhere in the transformation in case it is required. In the example given above, iterations stop as soon as tau reaches or goes over T arbitrary units of time. The use of the option iter allows to have a fine control over the number of applications of the transformation.

Figure 3 gives the results of three executions of the MGS program above (the instructions required for the output have not been included). Stochastic constants C i and C i have been determined by biological experiments detailed in [START_REF] Kuttler | Gene regulation in the pi calculus : simulation coperativity at the lambda switch[END_REF]. The first plot shows the system in a state that has not yet evolved into a lytic or lysogenic phase; the next plot shows the system after a switch has occurred in a lytic phase (only molecules of CRO remain); the last plot shows the system in the lysogenic phase (only molecules of CI remain). Each simulation is run until 2000 arbitrary units of time are reached for the Phage transformation on an initial state consisting of three copies of the CI protein, one copy of the DNA and ten copies of the RNA polymerase. On 500 runs of the simulation, the lytic phase dominates in 51% of the cases and the lysogenic phase dominates in 39% of the cases.

Modeling of an Homogeneously Distributed Population of Bacteria

The previous chemical model only considers a single isolated bacterium and the reactions taking place in that bacterium. One of the interesting phenomena of the λ phage infection [START_REF] Arkin | Stochastic kinetic analysis of developmental pathway bifurcation in phage λ-infected Escherichia coli cells[END_REF] is that two subpopulations (following either a lytic or a lysogenic phase) rapidly segregate after the introduction of the virus, that is, within a time interval of the same order than the cell division time which is about 20 minutes. We propose to describe in this section the behavior of a population of bacteria through a quantitative and stochastic model relying on an extension of the Gillespie's first reaction method to nested compartments. As a matter of fact, we consider that all bacteria are homogeneously distributed in a "soup".

A Simple Model of Infection Diffusion

We propose to study the segregation and the growth of the population during an infection process. Here, our main objective is to reuse the genetic switch model described and implemented in the previous section and coupling it with a model of infection diffusion taking place at the level of a population of bacteria. For the sake of simplicity, we keep our model as concise as possible.

This model involves a population of two kinds of entity that freely diffuse: bacteria and viruses. A bacterium can be in one of the following four different states:

1. Safe when it is not infected by a phage. 2. Infected when it has just integrated an external virus by endocytosis. We suppose that only safe bacteria can be infected. 3. Lytic when the infection has led to a lytic phase. 4. Lysogenic when the infection has led to a lysogenic phase.

To represent the growth of the population, we consider the two following behaviors of a bacterium:

-Division: all cells may reproduce. We assume that both daughter cells are sharing the same state of their mother cell. -Death: there are two ways for a cell to die, either when it has reached a certain number of divisions (because of cell aging), or when it is a lytic bacterium that "explodes" releasing its virions in the population.

Viruses are considered as "passive" entities which are degradated.

Representation of a Population of Bacteria

The simulation of the population model requires to deal with chemical reactions that are not uniformly taking place in space: the cell membrane of each bacterium separates inner chemical reactions from their outside environment. A bacterium plays the role of a compartment, and compartments break the hypothesis of homogeneity of the system required by Gillespie's first reaction method. Nevertheless, Gillespie's method can be extended to compartmentalized systems under the following reasonable hypotheses:

-each compartment is an insulator that makes inner and outer reactions being independent; -chemical reactions take place homogeneously within each compartment; -chemical reactions take place homogeneously in the outside of the compartments, that is, at the level of the population.

This last hypothesis requires that the population of compartments is well mixed. In this model, we assume that segregation is space-independent as the infection is uniformly spread at the initial state.

Extension of Gillespie's First Reaction Method. The previous hypotheses ensure that Gillespie's approach can be independently considered at the level of the population and inside each compartment. Following the same reasoning as in section 3.2, let us consider a compartment σ and a reaction µ, and let us respectively denote c σ µ , h σ µ and a σ µ , the stochastic reaction constant of reaction µ in σ, the number of molecular combinations that can activate reaction µ in σ, and the propensity of reaction µ in σ. We can then define the probability P (τ, µ, σ)dτ that µ is the next reaction occuring in compartment σ in the time interval (t + τ, t + τ + dτ ) as:

P (τ, µ, σ)dτ = a σ µ e -b 0 τ dτ
where the quantity b 0 = P θ a θ 0 = P θ,ν a θ ν corresponds to the combined propensity of the whole system. Consequently, a new algorithm can be defined in the same way as the original Gillespie's first reaction method: an elapsed time τ is chosen for each compartment σ and each reaction µ w.r.t. probability P (τ, µ, σ)dτ , and the couple (µ, σ) with the lowest elapsed time is considered.

Integration of the Extended First Reaction Method in MGS.

Whereas nested compartments can be obviously handled in MGS as nested multisets, the extended algorithm cannot be easily implemented as a transformation rule application strategy. By contrast, it can be programmed as a function that uses the previously defined transformation Phage. This point of view has already been investigated in [START_REF] Spicher | Stochastic p systems and the simulation of biochemical processes with dynamic compartments[END_REF] where the authors propose to simulate biochemical processes with dynamic compartments in MGS using the formalism of P systems. Based on these results, we implement our model for the phage infection at the level of the population of bacteria.

We first have to define the transformation Env of reactions taking place at the population level:

trans Env = { 'D ={ C = C degradation }=> . b:bact ={ C = C death }=> . b:bact ={ C = C division }=> division(b) b:safe, 'D ={ C = C infection }=> 'D, #3 'CI, b b:lytic ={ C = C lytic }=> #10 'D }
Viruses are represented by their DNA in the environment. The rules of transformation Env respectively specify the natural degradation of viruses, the local evolutions of bacteria (death and division), the infection of "safe" bacteria by a virus, and the death of bacteria in lytic phase. The predicate bact, safe and lytic are used to distinguish the different states of a bacterium.

The implementation of the extended first reaction method consists in computing the local evolution of each compartment and then only considering the faster evolving compartment. We do not describe the MGS code of the CompartmentsEvolve that is straightforward: it consists in the distribution of the Phage transformation onto each cell and the selection of the fastest reaction.

Given the current state of the system, the final function ExtendedFirstReaction computes the new state of the system by choosing if the evolution has to take place at the level of a compartment or at the level of the population: The whole program (including instructions for producing the output) is about 150 lines of MGS code.

A Simulation Example

Figure 4 presents some simulations of this first model of infection diffusion. Parameters have been chosen in order to observe an exponential growth of a safe population (see left plot on Fig. 4(a)) and the same proportion of lysis and lysogeny in a population of infected bacteria (see on right plot of Fig. 4(a), that mean numbers of CI 2 and CRO 2 per bacteria are somehow the same).

The protocol of our study has consisted in running 500 simulations of the model, during 40 arbitrary units of time, and starting with a population composed of 9 safe bacteria and 1 lytic bacterium. When the lytic cell dies, it releases some viruses and makes the infection start. The simulations have underlined four different outcomes at the end of the simulation:

1. Failure of the infection: no virus remains in the environment and all bacteria are safe. This happens with frequency of about 10%. 2. Death of the population: as shown in Figure 4(b), majority of cells enter the lytic phase (CRO predominates when all the population is infected) and die. Then, remaining lysogenic cells are to few (about 3 on the plot) to maintain the population, and disappear due to the stochastic noise. This appears with frequency of about 12%. 3. Homogeneous lysogenic population: Figure 4(c) presents this case where all bacteria become lysogenic (CRO is not expressed anymore after time 25 and only CI remains). The population keep on growing exponentially; the phage DNA is silently replicated with each cell division. This outcome represents about 40% of the simulations. 4. Heterogeneous lysogenic and safe population: this interesting outcome is the more scarce (about 8% of the cases). The behavior is quite similar to the previous one, but the lysogenic phase appears while the whole population has not been infected yet (see Fig. 4(d)). It results a mixed population with safe and lysogenic bacteria.

Last 30% of the simulations corresponds to states that have not reach one of the previous 4 outcomes within the 40 units of time.

Modeling of an Heterogeneously Distributed Population of Bacteria

The previous description of the phage infection was a first step into the modeling at the level of the population of cells. However, the spatial distribution of cells and viruses was homogeneous. Because of that property, we have been able to extend and use Gillespie's algorithm to nested compartments. Nevertheless, this modeling method has two essential drawbacks. The first drawback is of practical importance: Gillespie's algorithm becomes slower as the number of molecular species increases. Consequently, the complexity of the extended first reaction method increases with the number of bacteria in the population. If the model of a bacterium involves N molecular species, simulating a population of n bacteria is equivalent to handling a solution involving n × N distinct molecular species. In this case, quantitative simulations require computation power that is rapidly out of reach for serious simulations (in size and duration).

The second drawback lies on the strong assumption of uniform and homogeneous organization of the system at all scales (from compartments to populations). The (multi-)set topology together with Gillespie's strategy capture the Brownian motion of well-mixed solutions. Nevertheless, this implicit representation of space cannot be easily extended to heterogeneous systems where space matters. Furthermore, it is not possible, with such models to capture spatial properties of biological systems like pattern formation, spatial segregation, . . .

With the second kind of model presented in this section, we focus on this last issue. We consider a population of spatially distributed bacteria and propose a model allowing the study of the impact of cells organization on the evolution of the population. Such kind of model is of major importance since it may drastically change the interpretation of the considered biological system, as it has been shown in [START_REF] Shnerb | The importance of being discrete -life always wins on the surface[END_REF].

We propose in the next sections two models based on two different kinds of spatial representation: a model of population of cells on a static medium implemented as a cellular automaton and an individual-based model evolving on a dynamic graph topology.

Regular and Static Space Model: Cellular Automaton

Cellular automata provide a very simple way to handle spatial phenomena in biological simulations, as long as the topology of the underlying space do not evolve. We detail first how they can be considered as a rule-based computation model. Then, we present a model of the phage infection based on the use of a cellular automaton.

Cellular Automata in MGS

A cellular automaton (CA) is a regular lattice of cells, where each cell is characterized by a state taken from a finite set. The global evolution of the CA consists in applying synchronously on each cell, a local evolution function that computes the new state of the cell as a function of its current state and of the states of the cells in its neighborhood.

This computation model is naturally translated into a rule-based formalism:

-The regular lattice corresponds to a topological collection with a regular neighborhood. The MGS language provides a specific kind of topological collection called Group Based Fields (GBF for short), to deal with such a regular topology (see references (28; 29) for details on GBFs). -The local evolution function can be encoded in rules of the form: s => s if some condition c holds on the neighborhood Such a rule specifies the transition of a cell from the state s to the state s where the condition c defines the context for the application of the rule. This rule can be straightforwardly translated in an MGS transformation rule: where s and s are constants of the language and P c is a predicate that implements condition c. -As all the cells have to be updated in a synchronous mode, the local evolution rules have to be applied all together and everywhere: this corresponds to the maximalparallel rule application strategy.

s

Symbolic Model of an Infection Spreading

There are many ways to define a CA model of a biological system [START_REF] Ermentrout | Cellular automata approaches to biological modeling[END_REF]. We propose here an abstract and symbolic model of the infection spreading where the CA allows to underline the spatial organization.

To achieve this goal, we consider a CA based on a 2D square lattice. As summarized on Figure 5, the infection is characterized by the state of a cell: a cell is either safe or infected or lytic or lysogenic; the death of a lytic bacterium leads to the release of viruses. The local dynamics consists in updating this infection state using the rules given in Figure 5:

-a safe cell becomes infected with a probability Rv equal to the ratio of virus cells in its neighborhood; -an infected cell becomes either lytic or lysogenic with respect to the chemical model given in section 3; -a lytic cell dies and releases viruses.

The translation in the MGS formalism is straightforward. The following expressions define the different types of value used to represent the state of the CA: These declarations defines the type of the considered lattice as a GBF with values of type States. States are represented by symbols. The GBF topology is the Cayley graph of the finite presentation of a commutative group. The group is generated by the elementary displacements available to "move" from one node of the graph to one of its neighbors. The Grid GBF is defined by a finite presentation: a finite set of generators and equations between the generators. Here, four displacements (generators) are considered following the North, the South, the East, and the West direction. Knowing that the group is commutative and that North and South are opposite, as well as East and West, the Cayley graph of the group is isomorphic to a 2D grid with a von Neumann neighborhood.

The dynamics is implemented by the following transformation describing the three rules of evolution: The whole program (including instructions for producing the output) is about 100 lines of MGS code. In this transformation, we should emphasize that:

trans LocalEvolution = { 'Safe / (random(1.0) < R v ) =>
-The first rule is guarded by a predicate that depends on the value of a random number. This allows to apply the rule with probability R v as required by the model. -In the second rule, the choice between lytic and lysogenic phase depends on the call of the previously defined transformation Phage using Gillespie's strategy. -Finally, if a cell is not updated using any of the 3 rules of the transformation, by the semantics of the transformation application, it remains unchanged.

A Simulation Example

Pictures in Figure 6 show some screen shots of runs of the CA model simulations. As the initial state is designed with a single infected cell (see Fig. 6(a)), the infection propagates as a concentric wave (see Fig. 6(b)) leading to a fixed point.

One may expect a fixed point where all CA cells are either in a lysogenic state or in a viral state (see Fig. 6(c)). Surprisingly, another kind of fixed point also occurs where lysogenic bacteria surround the cluster of viruses and stop the spread of the infection (see Fig. 6(d)). This outcome appears for sets of parameters C i , C i leading to a majority of lysogeny.

In order to study this phase transition, we have simulated the model for different lytic rates. The lytic rate represents the probability for an infected bacterium to switch to a lytic phase. The protocol of these simulations has consisted in 400 runs of the simulation for each lytic rate, varying from 0% to 100%.

Figure 7(a) shows that for lytic rates lower than 65%, there is always a formation of a cluster (the probability is equal to 1); for lytic rates higher that 95%, the infection is never stopped. Between these two thresholds, a phase transition occurs quite quickly (in a range of 30%). Figure 7(b) shows that even if a cluster appears for lytic rates lower than 60%, the activity in the CA (due to the lytic bacteria) is not constant and increases with the lytic rate to reach a maximum (at about 150 iterations). By contrast, when no cluster appears, the system stabilizes faster (at about 75 iterations) while the infection covers the whole lattice. It is interesting to see that the activity and the size of the clusters (given Fig. 7(c)) are not correlated showing that lytic rates of about 60% increase the stabilization time.

Amorphous and Dynamical Space Individual-Based Model

While the previous model exhibits the interesting effect of spatial aggregation, one could wonder whether it is meaningful to have a square lattice representing (and therefore constraining) the neighborhood of the cells. Furthermore, important features like cell replication and cell death (except for lytic cells) are not taken into account. Actually, considering such behaviors in a CA requires to deal with way more complex rules that leave the scope of standard CA description to more complex CA such as Margolus neighborhood or lattice gas automata.

In order to study how the propagation of the infection is related to the population growth, we propose to extend the previous model to a more phenomenological model with a dynamical topology where the spatial coordinates of cells matter. Such a topology is available in MGS using the Delaunay topological collections. Even if this kind of modeling falls in the field of agent-based models (ABM) or multi-agent models (MAS), we consider it as being part of what can be described and implemented in a rule-based language like MGS. We develop this point of view in the conclusion of this article.

Individual-Based Model on a Delaunay Graph

By contrast with the complete neighborhood of multisets, this collection (automatically) computes the neighborhood between the collection's elements according to their spatial coordinates. The neighborhood is computed using a Delaunay triangulation of the coordinates in R n . This kind of neighborhood together with a mechanical model has already been successfully used in biology for the modeling of growth of cells (31; 32). The mechanical model consists in a mass/spring system. It allows to keep a coherent global structure where cells are always close to each other and are able to "push" their neighbors if they lack space, or to fill holes in the structure.

A Phenomenological Model

The model is defined by extending the previous CA model. A cell of the CA lattice corresponds now to a punctual mass that is localized in the R 3 space. A mass is characterized by its spatial position, its velocity, and its biological state (i.e., a value of type States since it represents a biological entity). The neighborhood relation of the CA is replaced by springs between masses, and the presence of a spring between two masses depends on the neighborhood induced by the Delaunay triangulation of the masses. The biological evolution of the masses corresponds to the update of their biological states with respect to the previous transformation LocalEvolution extended with 3 rules to express cellular division, death, and virus degradation. The spatial evolution4 of the masses corresponds to the definition of a new transformation SpatialEvolution that computes the movement of a mass due to its incident springs during a time interval ∆t. Finally, the whole evolution step consists in composing the transformations:

SpatialEvolution[*] • LocalEvolution
Here, the mechanical model is iterated until an equilibrium is reached to consider a quasi-static approximation: the time-scale of the mechanical model is considered faster than the biological one. The whole program (including instructions for producing the output) is about 200 lines of MGS code.

A Simulation Example

Figure 8 shows four graphical interpretations of simulations of this model using the Imoview (33) visualization tool developed in the MGS project. In these simulations, viruses and bacteria do not have the same size: viruses are smaller than cells and thus allow a more efficient diffusion in the population. This prevents from the formation of clusters. Nevertheless, another kind of behavior appears as shown in Fig. 8(d). The population is split into two subgroups: a safe one and a lysogenic one. This outcome is quite similar to the behavior shown on Fig. 4(d) for the nested compartment model. It should be of importance in real infections since lysogenic cells are immunized towards further infection. If a lysis is induced in the lysogenic population, the safe bacteria will prove to be a good source of hosts for the released viruses.

Conclusion: towards Spatial Systems Biology

Rule-based programming has been heavily investigated in the field of natural systems modeling. It was an essential contribution of (34) that has started an important effort in the π-calculus community to focus on problems raised by systems biology (35; 36; 37; 38; 39; 40; 41). This approach has been very efficient and successful, yet only limited to the consideration of biological processes that did not involve space: Brownian motion is supposed to allow each entity to encounter all other entities. Other formalisms have additional structure and spatial relationships: P systems consider nested compartments, L-systems allow through the interpretation of their symbols the modeling of complex spatial relationships, etc. but to our knowledge, besides vertex-vertex systems [START_REF] Smith | Local specification of surface subdivision algorithms[END_REF] which offer rewriting on graph topologies, no language nor formalism in the field of natural systems modeling has considered the explicit representation of space together with an implicit handling in evolution rules.

The MGS project is an attempt to combine rule-based programming with topological collections in a effort to offer a versatile, expressive and efficient mean for the modeling of natural systems. Rules allow to focus on the interactions between the elements of the systems while topological collections provide an explicit description of space. Here, we only have considered a small subset of the topologies available in MGS [START_REF] Michel | There's plenty of room for unconventional programming languages or declarative simulations of dynamical systems (with a dynamical structure)[END_REF].

Nevertheless, we have shown that it is sufficient for the modeling of the biological process at various levels in systems biology: our paradigmatic example of the switch of the λ phage has been described from the transcription network at the level of the biochemical interactions to the level of population of cells embedded in a general three dimensional space. Furthermore, while our approach is clearly individual based, classical ABM/MAS (like NetLogo (43)) focus on the description of agents behavior, rulebased languages and formalisms promote the description of the interactions between the individuals. Among the many consequences of this point of view, it allows us to consider any kind of interactions between individuals (and not only interactions between two agents); since the rule application strategy is explicit, it makes possible to have a fine control over the scheduling of the selection of the individuals in interactions; by changing the type of topological collection, it is easy to change the nature of the underlying space where the individuals are embedded, etc.

Our running example is a very well studied biological process at the level of gene regulation. It is a classical case in systems biology which aims at integrating processes at various time and spatial scales into a single and coherent formal description. However, while classical molecular systems biology focuses on gene or protein interaction networks, cellular and supra-cellular organization levels are rarely considered. Since the fate of biological systems is not only determined by genes, the spatial organization of cells, tissues and organs plays a key role in most physiological processes (see for example [START_REF] Farge | Mechanical induction of Twist in the Drosophila foregut/stomodeal primordium[END_REF]) and must therefore be a part of the models. Recognizing this importance by integrating spatial properties extends systems biology towards spatial systems biology.

We believe that languages and tools have to be further developed to take into account these spatial interactions, and that rule-based languages, like MGS, can play a key role. More specifically, the generalization of topological collections using notions from algebraic topology [START_REF] Giavitto | Topological rewriting and the geometrization of programming[END_REF], allows us already to model complex spatial interactions. The relevance of these concepts is currently under validation through various large scale examples: simulation of morphogenesis processes [START_REF] Barbier De Reuille | Computer simulations reveal novel properties of the cell-cell signaling network at the shoot apex in arabidopsis[END_REF], amorphous and autonomic computing examples and simulation at various levels of the behavior of genetically engineered bacteria in the french team of MIT's iGEM competition [START_REF]iGEM: Modeling a synthetic multicellular bacterium. Modeling page of the Paris team wiki at iGEM'07[END_REF].

Fig. 1

 1 Fig. 1 Description of the regulatory region of the genetic switch.

Fig. 2

 2 Fig. 2 Regulation network of the genetic switch of the λ phage. Monomers CI and CRO are described in circles and diamonds; dimers are composed by monomers. Thin dotted arrows represent transcription and translation; bold dotted arrows represent (de)polimerization; thin arrows represent degradation; bold arrows represent regulation. The drawing follows the description of the DNA given in Figure 1.

  BactInit := 'D :: #10 'P :: #3 'CI :: bag:() This code defines a new variable, BactInit. The empty multiset is represented by bag:(), and operator :: inserts an element in the given multiset. The expression #n v represents n copies of the value v. Chemical Reactions Representation and Evolution. Each chemical reaction is translated into a transformation rule (or two if the reaction is reversible) that is characterized by an option C representing the stochastic constant of the reaction. For example, reaction (2) corresponds to the two following MGS rules #2 'CI ={ C = C 12 }=> 'CI2 and 'CI2 ={ C = C 21 }=> #2 'CI Consequently, the dynamics is captured by the following set of rules in the Phage transformation: trans Phage = { (* Rules for CI *) (* Rules for CRO *) 'CI ={ C = C 0 }=> .; 'CRO ={ C = C 0 }=> .; #2 'CI ={ C = C 12 }=> 'CI2; #2 'CRO ={ C = C 12 }=> 'CRO2; 'CI2 ={ C = C21 }=> #2 'CI; 'CRO2 ={ C = C 21 }=> #2 'CRO; 'D0, 'CI2 ={ C = C 1 }=> 'D1; 'D0, 'CRO2 ={ C = C 1 }=> 'D'3; 'D1 ={ C = C -1 }=> 'D0, 'CI2; 'D'3 ={ C = C -1 }=> 'D0, 'CRO2; 'D1, 'CI2 ={ C = C2 }=> 'D2; 'D'3, 'CRO2 ={ C = C 2 }=> 'D'2; 'D2 ={ C = C -2 }=> 'D1, 'CI2; 'D'2 ={ C = C -2 }=> 'D'3, 'CRO2; 'D2, 'CI2 ={ C = C 3 }=> 'D3; 'D'2, 'CRO2 ={ C = C 3 }=> 'D'1; 'D3 ={ C = C-4 }=> 'D2, 'CI2; 'D'1 ={ C = C -3 }=> 'D'2, 'CRO2; 'D2, 'P ={ C = C t }=> 'D2, 'P, 'CI; 'D0, 'P ={ C = C t }=> 'D0, 'P, 'CRO; 'D'3, 'P ={ C = C t }=> 'D'3, 'P, 'CRO } MGS integrates the first reaction method as a transformation rule application strategy. The transformation Phage is called and iterated using the Gillespie's strategy of MGS by the following expression: Phage[iter = (tau <= T ), strategy = 'gillespie](BactInit)

Fig. 3

 3 Fig. 3 Three different runs of the simulation, plotted with Gnuplot of the switch of the lambda phage. The light (resp. dark) curves show the number of CI (resp. Cro) molecules as a function of time (arbitrary units).

fun

  ExtendedFirstReaction(pop) = ( let pop1, tau1 = CompartmentsEvolve(pop) in let pop2, tau2 = Env[strategy = 'gillespie, postlude = (fun x -> x,tau)](b) in if (tau1 < tau2) then pop1 else pop2 fi )

  (a) Simulations illustrating the set of parameters. (b) Simulation leading to the death of the whole population.

  (c) Simulation leading to a purely lysogenic population.

  (d) Simulation leading to a mixed lysogenic/safe population.

Fig. 4

 4 Fig. 4 Results of simulations of the nested compartments model.

Fig. 5

 5 Fig. 5 States and evolution rules of the cellular automaton. The top of the figure describes the five different states of a cell; the middle and lower figures describe the evolution rules for a safe cell, an infected cell and a lytic cell.

  type States = 'Safe | 'Infected | 'Lytic | 'Lysogenic | 'Virus and gbf Grid = < North, South, East, West | North = -South, East = -West > and collection Lattice = [States]Grid

( a )

 a Initial state: a lytic cell in the center of a population of safe bacteria. (b) State of the system after 25 iterations. (c) Fixed point after 75 iterations: no safe cell remains. (d) Fixed point ater 150 iterations: formation of a cluster.

Fig. 6

 6 Fig.6Results of the simulation of the CA model on a toric 20x20 lattice. The color code of the four pictures is: safe bacteria are plotted in green, infected bacteria in red, lytic bacteria in blue, lysogenic bacteria in yellow, and viruses in dark blue.

  (a) Probability of cluster formation according to the lysis rate. (b) Mean number of iterations for stability according to the lysis rate. (c) Mean size of the clusters according to the lysis rate.

Fig. 7

 7 Fig. 7 Results of the simulation of the CA model on a toric 20x20 lattice. On the three graphs, plots show some properties of the system according to the lytic rate of the model. The size of the error-bars clearly shows the effects of stochasticity on the results.

Fig. 8

 8 Fig. 8 Results of the simulation of the dynamical space model. Safe bacteria are plotted in green, infected bacteria in red, lytic bacteria in blue, lysogenic bacteria in yellow, and viruses in dark gray.

  and in what order). For this purpose, we write e ⇒ R,St e to characterize the rule application strategy where St represents the considered policy. For instance, St may mean that only one rule R i , chosen in a non-deterministic way, is applied; or that rules are applied all together on different subparts of e in a concurrent mode. Considering a set of rules R together with a rule application strategy St, we define the relation ⇒ * R,St as the reflexive and transitive closure of ⇒ R,St . In other words, if e and e are two entities, e ⇒ * R,St e means that there exists a finite sequence of entities e 1 , . . . , e n such that e = e 1 , e 1 ⇒ R,St e 2 , . . . , e n-1 ⇒ R,St e n , e n = e

  as x / P c (x) => s

		with probability Rv		Phage leads a lysogeny
		Phage leads a lysis
	Safe	Infected	Lytic	Lysogenic	Virus

The Web site of the project is http://mgs.ibisc.univ-evry.fr

Evaluating the stochastic constants is one of the hardest task in stochastic simulations of biochemical reactions. The interested reader should refer to(22; 23) that describe two experiences in that field.

Because of the lack of space, we do not include here the MGS code for this model. Nevertheless, the code of all the examples given in this article is available upon request from the authors.

Acknowledgments. The authors gratefully acknowledge the reviewers for their valuable comments on a first version of this article.

This work is partially funded by the BQR of the University of Évry on "Synthetic Biology".


Dieter Fiems
email: dieter.fiems@telin.ugent.be

Koen De Turck
email: koen.deturck@telin.ugent.be

Queues with M/G/∞ input: a survey and some new results

As input traffic at various nodes in packet switched telecommunication networks typically exhibits considerable correlation, there is a continuing interest in queuing systems with correlated arrivals. The class of M/G/∞ arrival models can account for various types of arrival correlation, has a twolevel structure which is appealing in the context of network protocol stacks, and is characterised by only few parameters. Moreover, queues with M/G/∞ input prove to be amenable for a classic queuing analysis. In view of these characteristics, it is not surprising that these models have attracted considerable attention in the past. In this paper, we survey the literature on M/G/∞ type queuing systems, present some new results on these systems, and extend the class of M/G/∞ input models, while retaining the analytic tractability of the corresponding queuing systems.

I. INTRODUCTION

Arrival correlation significantly affects queuing performance and hence there is a continuing interest in tractable queuing models which can accurately capture the correlation in the arrival process. Arrival models of interest include, amongst others, finite state-space Markovian arrival models like the discrete batch Markovian arrival model as well as Markovian arrival models with an infinite but structured state-space. A prime example of the latter class is M/G/∞ input -often referred to as train arrival or session models -which is the subject of this paper.

A continuous-time version of M/G/∞ input was already investigated by Cohen in 1974 [START_REF] Cohen | Superimposed renewal processes and storage with gradual input[END_REF], while the origins of te discrete-time M/G/∞ input studied here can be traced back to the works of Cox [START_REF] Cox | Statistics: an appraisal, chapter "Long-range dependence: a review[END_REF]. Formally, let P = {P n , n ∈ Z} be a sequence of independent and identically distributed Nvalued random variables with common mean λ and ith order moment λ i . Moreover, let S = {S n,k , n ∈ Z, k ∈ N * } be a doubly-indexed sequence of independent and identically distributed N * -valued random variables with mean µ and ith order moment µ i . The M/G/∞-input process is then defined as follows,

A n = ∞ k=0 P n-k =1 1 {S n-k, >k} . (1)
Here 1 {x} is the indicator function which evaluates to 1 if x is true and to 0 otherwise. The sequence A = {A n , n ∈ Z} is stationary ergodic and corresponds to the queue size process of a discrete-time M/G/∞ with arrival process P and service process S. In literature, P is sometimes assumed to be a sequence of Poisson distributed random variables. However, this restriction will not be imposed unless explicitly stated in the remainder. The two-level structure of M/G/∞ input is a natural abstraction of an upper and a lower layer in a network protocol stack. At the upper layer, a connection or session is established which continues for a certain duration. The session lengths correspond to the service times of the M/G/∞ queue. During the session, packets are then transmitted by the lower layer in the protocol stack. This natural abstraction is not the only reason M/G/∞ traffic has attracted much attention. Equally important is the fact that M/G/∞ processes are versatile in capturing correlation in the arrival process while they are characterised by only a few parameters: the mean arrival rate λ (assuming Poisson arrivals) and the session length distribution S(n) = Pr[S 0,1 ≤ n], see [START_REF] Cox | Statistics: an appraisal, chapter "Long-range dependence: a review[END_REF]- [START_REF] Resnick | Data network models of burstiness[END_REF]. Finally, we mention the analytical tractability of single server queues with M/G/∞ input as an additional motivation for studying these arrival processes.

The remainder of this paper is organised as follows. The next section reviews the stochastic characteristics of M/G/∞ input. Basic results on the single-server queue with M/G/∞ input are then presented in Section III. Despite the complex state space of the M/G/∞ process, the expression of the mean queue content of this queuing system is surprisingly simple. This observation leads us to extend the class of M/G/∞ processes while retaining the tractability of the mean queue content. This class of regenerative arrival processes is introduced in Section IV and includes discrete autoregressive processes [START_REF] Kamoun | The discrete-time queue with autoregressive inputs revisited[END_REF] as well as branching arrival processes [START_REF] Fiems | Queues with Galton-Watsontype arrivals[END_REF] in addition to M/G/∞ input.

In terms of applicability, the single server model suffers from the inherent drawback that a single active session always consumes all available bandwidth. To avoid this, either sessions should not continuously produce packets or more servers need to be introduced. In the former case, the arrival process no longer satisfies the regeneration property and no simple expression for the mean queue content can be found while in the latter case, regeneration does not suffice for tractability. Therefore, the focus shifts to approximations for various performance measures. Light-traffic approximations are considered in Section V. Finally, heavy-traffic limits are the subject of Section VI.

II. PROPERTIES OF M/G/∞ INPUT

Taking expectation in (1) immediately yields,

E[A 0 ] = λµ .
Analogously, the variance of A 0 can be calculated directly from [START_REF] Cohen | Superimposed renewal processes and storage with gradual input[END_REF],

var(A 0 ) = λµ + (λ 2 -λ -λ 2 ) ∞ k=0 (1 -S(k)) 2 .
More generally, the covariance of the M/G/∞ input is given by,

cov(A 0 , A n ) = λ E[(S 0,1 -n) + ] + (λ 2 -λ -λ 2 ) ∞ k=0 (1 -S(k))(1 -S(k + n)) .
Assuming P constitutes a sequence of Poisson distributed random variables, these expressions simplify considerably. In this particular case, the distribution of A 0 is known: A 0 is Poisson-distributed with mean λµ. Moreover, the covariance of the M/G/∞ input simplifies to

cov(A 0 , A n ) = λ E[(S 0,1 -n) + ] = λ ∞ k=0 Pr[S 0,1 > n + k] ,
see e.g. [START_REF] Krunz | Modeling video traffic using M |G|∞ input processes: a compromise between Markovian and LRD models[END_REF], [START_REF] Parulekar | Tail probabilities for M/G/∞ input processes (I): Preliminary asymptotics[END_REF].

The covariance function allows for adapting the service time distribution to fit the autocorrelation function of e.g. measurement data. For example, in [START_REF] Krunz | Modeling video traffic using M |G|∞ input processes: a compromise between Markovian and LRD models[END_REF], the service time is chosen to fit the autocorrelation of a number of video traces. M/G/∞ proves to be very versatile in this respect. In fact, if the autocorrelation function ρ(n) = cov(A 0 , A n )/ var(A 0 ) is decreasing and integer convex with 1 > ρ(1) and lim n→∞ ρ(n) = 0, then the probability mass function of the service times of the M/G/∞ input can be expressed in terms of the autocorrelation function as follows [START_REF] Poon | A refined version of M/G/∞ processes for modelling VBR video traffic[END_REF],

Pr[S 0,1 = k] = ρ(k -1) -2ρ(k) + ρ(k + 1) 1 -ρ(1) .
In other words, the service time distribution can be adapted to match the autocorrelation function exactly for a large class of processes.

III. THE SINGLE SERVER QUEUE

We first focus on the discrete-time single-server queue with M/G/∞ input and single-slot service times by means of classic probability generating functions techniques as investigated in [START_REF] Wittevrongel | Correlation effects in ATM queues due to data format conversions[END_REF]. The queue content at consecutive slot boundaries are related by the following Lindley recursion,

U n+1 = (U n -1) + + A n , (2)
whereby A n is defined in [START_REF] Cohen | Superimposed renewal processes and storage with gradual input[END_REF]. Obviously, (U n , A n ) does not constitute a Markov chain as all active session lengths need to be tracked as well. Therefore, let

H (k) n
denote the number of sessions that deliver their kth packet during slot n and let H n denote the infinite-length vector with elements

H (k) n . Then (H n , U n ) constitutes a Markov chain. For ρ = E[A 0 ] = λµ < 1,
the chain admits a stationary solution ( Ĥn , Ûn ). Let P(x, z) denote the joint probability generating function of the marginal distribution of this solution,

P(x, z) = E ∞ k=1 x Ĥ(k) 0 k z Û0 = E x Ĥ0 z Û0 , (3)
with x = [x 1 , x 2 , x 3 , . . .]
. This joint probability generating function then satisfies the following functional equation,

P(x, z) = P (x 1 z) z (P(C(xz), z) + (z -1)(1 -ρ)) , (4)
with

C(x) = [C 1 (x 2 ), C 2 (x 3 ), .
. .] a vector with elements,

C k (x) = 1 -(1 -z) 1 -S(k) 1 -S(k -1)
.

Iterating (4) yields an explicit expression for P, although it contains infinite sums and products. The moments of the queue content can then be obtained by the moment generating property of probability generating functions. In particular, the mean queue content equals,

E[ Û0 ] = ρ + µ 2 (λ 2 -λ) + (µ 2 -µ)λ 1 -ρ . (5)
In addition to the results above, single-server queues with M/G/∞ input have been studied under the assumption of specific session lengths, e.g. deterministic session lengths [START_REF] Xiong | Buffer contents and delay for statistical multiplexers with fixed-length packet-train arrivals[END_REF], [START_REF] Kamoun | Performance analysis of a discrete-time queuing system with a correlated train arrival process[END_REF] or geometrically distributed session lengths [START_REF] Bruneel | Performance study of statistical multiplexing in case of slow message generation[END_REF], [START_REF] Feyaerts | Session delay in file server output buffers with general session lengths[END_REF]. Extending M/G/∞ input, the queuing model in [START_REF] Vuyst | Statistical multiplexing of correlated variable-length packet trains: an analytic performance study[END_REF] and [START_REF] Vuyst | Mean value and tail distribution of the message delay in statistical multiplexers with correlated train arrivals[END_REF] allows for arrival correlation of the sessions. A twostate Markov chain is introduced to capture the correlation of session arrivals (an interrupted batch arrival process IBP); this is a queue with IBP/G/∞ input. In all references above, a first-in-first-out scheduling discipline is assumed. This is not the case in [START_REF] Walraevens | A discrete-time priority queue with train arrivals[END_REF] where priority disciplines are investigated under the (simplifying) assumption of geometrically distributed session lengths.

Note that, despite the size of the state space of the M/G/∞ process, equation ( 5) above is intriguingly simple and only contains the first and second order moments of the arrivals and the session lengths. This observation has sparked the investigation of the regenerative arrival processes introduced next.

IV. REGENERATIVE ARRIVALS

The analysis of the queuing model of the preceding section makes use of the following property. If the queue is empty, then there cannot be active sessions. Indeed, if this were the case, there would be packets present in the queue. In mathematical terms, this means that the arrival process regenerates whenever the queue is empty. Regeneration turns out to be the key property which makes the M/G/∞ single-server queue analytically tractable. In fact, regeneration is the only property required to make the mean queue content analysis tractable as shown in this section.

That is, we now investigate the single-server queue with single-slot service times for the following arrival process. A = {A n , n ∈ Z} is a stationary ergodic sequence of nonnegative integer random variables which adheres the following assumptions.

[A1] The arrival process regenerates when there are no arrivals. That is,

Pr[A k+1 = i k+1 , . . . , A k+ = i k+ |A k = 0, F k -∞ ] = Pr[A k+1 = i k+1 , . . . , A k+ = i k+ |A k = 0] , (6)
for ∈ N * = {1, 2, . . .} and i j ∈ N for all j ∈ Z and where F k -∞ is the natural filtration of the arrival process

{A k }. [A2] Let f n = Pr[A n = 0, A n-1 > 0, . . . , A 1 > 0|A 0 =
0] denote the regeneration distribution. We assume the existence of its second moment,

∞ n=1 n 2 f n < ∞ .
The following is our main result on queuing systems with regenerative arrivals, see [START_REF] Fiems | Mean queue content of discrete-time queues with zero-regenerative arrivals[END_REF] for its proof. For stationary ergodic {A k } satisfying properties [A1] and [A2] and with E[A 0 ] < 1, the mean steady-state queue content E[ Û0 ] is given by,

E[ Û0 ] = E[A 0 ] -E[(A 0 ) 2 ] 2(1 -E[A 0 ]) + ∞ m=0 (E[A m+1 |A 0 = 0]-E[A 0 ]) + ∞ m=0 E[A 0 A m ] -E[A 0 ] 2 1 -E[A 0 ] . (7)
The mean queue content is finite provided the infinite sums in (7) converge. Seemingly somewhat artificial, many arrival processes adhere to the regeneration property. It holds for the above-cited M/G/∞ input, autoregressive and branching-type arrival processes [START_REF] Fiems | Queues with Galton-Watsontype arrivals[END_REF]. However, the class of these zero-regenerative arrival processes is considerably larger as illustrated by the following two examples.

A. Train arrivals with autoregressive session arrivals

A first example is the DAR(1)/G/∞ input model. In contrast to M/G/∞ input, new sessions arrive in accordance with a discrete autoregressive arrival process. This arrival model is characterised by the sequence S as defined in Section I and by (1) a sequence B = {B k , k ∈ Z} of iid Bernoulli distributed random variables with E[B 0 ] = p and (2) an iid sequence {N k } of N-valued random variables. N k denotes the number of new trains in slot k if B k = 0. The number of arrivals in slot k is still given by (1) but now S k is not an iid sequence. Instead, S k is expressed in terms of B k and N k as follows,

S k+1 = B k S k + (1 -B k )N k . (8)
By calculating E[A 0 A n ] and E[A n |A 0 = 0] for this arrival process, the following expression for the mean queue content is found,

E[ Û0 ] = E[N ] 2 µµ 2 -E[N ]µ 2 (1 -E[N ]µ + E[N ]µ(1 -2p) -E[N ] 2 µ 2 (1 -2p) (1 -p)(1 -E[N ]µ) + (E[N 2 ] -E[N ] 2 )µ 2 (1 + p) 2(1 -p)(1 -E[N ]µ) . (9)
For p = 0, the arrival model simplifies to the train arrival model where the number of new trains is a sequence of iid random variables. Plugging p = 0 in (9) yields [START_REF] Kamoun | The discrete-time queue with autoregressive inputs revisited[END_REF]. In addition, assuming single slot train-lengths -this means µ = µ 2 = 1 -we obtain the single-server queuing system with discrete autoregressive arrivals of [START_REF] Kamoun | The discrete-time queue with autoregressive inputs revisited[END_REF]. The expression of the mean queue content then simplifies to,

E[ Û0 ] = E[N 0 ](1 -3p)+E[(N 0 ) 2 ](1 + p)-2(1 -p) E[N 0 ] 2 2(1 -E[N 0 ])(1 -p) .

B. Stationary ergodic arrivals during trains

As a second example we consider a train arrival model in which trains produce packets in accordance with a stationary ergodic process. As before, each train generates at least one packet such that regeneration is ensured. However, trains may produce more packets. Let H k,n, denote the number of packets produced in the lth slot of the nth train arriving in slot k. The processes {H k,n, , = 0, 1, . . .} constitute a doublyindexed (by k and n) sequence of stationary ergodic positive random processes. Let A k denote the number of arrivals in slot k, we then have,

A k = ∞ m=0 P k-m n=1 1 {S k-m,n >m} H k-m,n,m . ( 10
) Let β = E[H 0,1,1 ], β 2 = E[(H 0,1,1
) 2 ] and κ n = E[H 0,1,0 H 0,1,n ] characterise the arrival process during sessions and let ρ = λµβ denote the arrival load. The mean queue content then equals,

E[ Û0 ] = - 1 2 λβµ 2 + ρ + λ ∞ m=0 (1 -S(m)) m n=0 κ n 1 -ρ + (λ 2 -λ)β 2 µ 2 -3λ 2 µ 2 β 2 -λµβ 2 2(1 -ρ) . (11)
Under the additional assumption that the number of packets produced in the consecutive slots is an independent sequence, the former expression further simplifies to,

E[ Û0 ] = - 1 2 λβµ 2 + 2ρ + λµ 2 β 2 -ρβ 2(1 -ρ) + (λ 2 -λ)β 2 µ 2 -3λ 2 µ 2 β 2 -λµβ 2 2(1 -ρ) , (12)
in correspondence with results published in [START_REF] Feyaerts | Session delay in file server output buffers with general session lengths[END_REF].

V. LIGHT TRAFFIC

We now focus on approximations of performance measures in case the regeneration property does not hold. In this section, we consider approximations when the queue is subject to light traffic. Performance of heavily loaded queues is the subject of the following section.

Recalling the construction of M/G/∞ input, there are two natural ways to reduce the arrival load. The reduction is either achieved by reducing the arrival rate of the sessions or of the amount of traffic produced in a session. The latter reduction can be accomplished by so-called π-thinning: every packet arrival in the original system also arrives in the system with reduced load with some probability q and is discarded otherwise.

Tsoukatos and Makowski [START_REF] Tsoukatos | Power-law vs exponential queueing in a network traffic model[END_REF] study approximations when the arrival rate of sessions is reduced. For the multi-server case, approximations for the probability that the queue is empty are obtained. Regarding the moments, a structural property is obtained for both single-server and multi-server systems is obtained, but an explicit expansion is obtained only for the single-server queue. However, the latter results have limited practical value as in this case an explicit expression is available.

Reducing the load by π-thinning is investigated in [START_REF] Turck | A Taylor series expansions approach to queues with train arrivals[END_REF] under the additional assumption that the session lengths are geometrically distributed. Under this assumption and without thinning, this model was independently investigated in [START_REF] Bruneel | Performance study of statistical multiplexing in case of slow message generation[END_REF] and [START_REF] Brandt | A single server model for packetwise transmission of messages[END_REF]. Let q denote the packet arrival probability and p denote the probability that a session continues, p = 1 -µ -1 . The following Taylor series expansions in q for the first two moments of the stationary queue content were obtained in [START_REF] Turck | A Taylor series expansions approach to queues with train arrivals[END_REF],

E[ Û0 ] = λ 1 -p q + λ 2 2(1 -p 2 ) + 2pλ 2 -(1 -p)λ 2(1 + p)(1 -p) 2 q 2 + p 2(1 -p)(p 2 + p + 1) λ 3 + (5p 3 + 3p 2 + 1)λ -p(1 + 2p)(1 -p) 2 2(1 -p) 2 (p 2 + p + 1)(1 + p) λ 2 + 4 p 4 -p 3 + p 2 -1 2(1 + p)(p 2 + p + 1)(1 -p) 2 λ 2 + 2p(1 + 2 p 3 )λ 2 -2 p 3 (1 -p) 2 2(1 + p)(p 2 + p + 1)(1 -p) 3 λ q 3 + O(q 4 ) , (13)
E[ Û 2 0 ] = λ 1 -p q + 3λ 2 2(1 -p 2 ) + 6pλ 2 -3(1 -p)λ 2(1 + p)(1 -p) 2 q 2 + (15p 3 + 13p 2 + 2 p + 3)λ -6p 4 + 6p 2(1 -p) 2 (1 + p)(1 + p + p 2 ) λ 2 - (9p + 2) 3(1 -p)(p 2 + p + 1) λ 3 + 2p(6p 3 + 2p 2 + 3)λ -12p 5 + 15p 4 -2p 3 + p 2 + p -3 2(1 + p)(p 2 + p + 1)(1 -p) 3 λ 2 - 9p 3 -2p -2 3(p 2 + p + 1)(1 -p 2 ) λ q 3 + O(q 4 ) . (14)
The methodology in [START_REF] Turck | A Taylor series expansions approach to queues with train arrivals[END_REF] allows for obtaining higher order expansions and expressions for higher order moments as well. However, the expressions for these expansions grow quickly in size and are therefore only practical for numerical evaluation.

VI. OTHER LIMITING RESULTS

The parsimony and versatility of the M/G/∞ input model is also apparent when considering various limiting regimes. Along with fractional Brownian motion and On/Off sources, M/G/∞ models have been among the most popular input processes, especially in the context of evaluating the impact of long-range dependence on buffer performance.

We summarize results from the fields of heavy-traffic and large-deviations theory [START_REF] Parulekar | Tail probabilities for M/G/∞ input processes (I): Preliminary asymptotics[END_REF], [START_REF] Tsoukatos | Power-law vs exponential queueing in a network traffic model[END_REF], [START_REF] Borst | Fluid Queues with Heavy-Tailed M/G/∞ input[END_REF], [START_REF] Resnick | Steady-state distribution of the buffer content for M/G/∞ input fluid queues[END_REF], [START_REF] Mandjes | Large deviations for Gaussian queues: modelling communication networks[END_REF], [START_REF] Mandjes | A note on queues with M/G/∞ input[END_REF]. As a fully rigourous exposition of the results contained in these papers, would easily exceed the total length of this survey, we sacrifice rigour for ideas and intuition. The queuing model considered in this section is the early-arrival and multi-server variant of (2):

U n+1 = (U n + A n -c) + .
Across the different methodologies, a dichotomy is observed between the cases that the input process is shortrange dependent (srd) resp. long-range dependent (lrd). Recall that a process A is srd (resp. lrd) if the covariance function v(k) . = cov(A 0 , A k ) is summable resp. non-summable. In [START_REF] Parulekar | Tail probabilities for M/G/∞ input processes (I): Preliminary asymptotics[END_REF] it is shown that

∞ k=0 |v(k)| = λ 2 (µ 2 + µ), (15)
hence the input process is srd (lrd) if and only if µ 2 is finite (infinite). Also recall the definition of regularly varying distributions: a mass function f (.) is regularly varying of index h, if and only if for all y > 0:

f (yt) f (t) → y h , t → ∞.
If h = 0 then we call f (.) slowly varying; and if h ∈ [0, 1) then f (.) is subexponentially varying. The heavy-traffic limit considered in [START_REF] Tsoukatos | Power-law vs exponential queueing in a network traffic model[END_REF] is as follows: consider a series of queuing systems indexed by r ∈ N. The rth model is fed by an input process with session arrival rate λ r , with λ r ↑ c/µ as r → ∞. Note that the distribution of the M/G/∞ service times, as well as the release rate c remains the same for every system in the sequence. The heavy traffic limit is then the distribution Q such that q r ∞ /w r ⇒ Q, where q r ∞ is the steady-state buffer content distribution for the rth system, w r is an appropriate scaling function and ⇒ denotes weak convergence.

If the input process is srd, then we find a convergence to Brownian motion, with formulas than go back as far as Kingman's. If on the other hand, S 0,1 has a Pareto distribution with parameter α, 1 < α < 2, then the solution is expressed in terms of a Mittag-Leffler special function, and the tail is found to be Pareto with parameter α -1. We note that the limiting process is Gaussian in the srd case, but non-Gaussian in the Pareto case.

Mandjes considers a different limit in [START_REF] Mandjes | Large deviations for Gaussian queues: modelling communication networks[END_REF]. In particular, he studies the Gaussian process that retains the first two moments of the original input process, that is, the average number of arrivals during stationarity E[A 0 ] and the covariance function. It is then shown that the Pareto case yields a logarithmic many-sources asymptotic the form 1 n log p n (b, c) → O(b α-1 ), provided that the buffer level is also large. This suggests Weibullian instead of Pareto tails for the buffer content, which is another evidence of the fact that the manner in which a limit of a series of stochastic processes is attained is very important.

Direct large-deviations analyses (without first transforming the input process into a Gaussian form) have been undertaken as well. For example, Parulekar and Makowski [START_REF] Parulekar | Tail probabilities for M/G/∞ input processes (I): Preliminary asymptotics[END_REF] obtain large-buffer asymptotics through large-deviation theory, whereas Duffield [START_REF] Duffield | Queueing at large resources driven by long-tailed M/G/∞-modulated processes[END_REF] and Mandjes [START_REF] Mandjes | A note on queues with M/G/∞ input[END_REF] obtain many-sources asymptotics. It is observed that srd inputs invariably lead to exponential tails, whereas with lrd crucial differences are caused by the shape of the session duration distribution. If session durations are subexponentially varying, then the input rate during overflow will significantly exceed the drain rate c, and the time to overflow is proportional to the buffer level. In contrast, with slowly-varying inputs, the buffer will fill very slowly, at a rate only slightly higher than c. Also note that for srd M/G/∞ inputs, the machinery developed for logarithmic asymptotics for general single-server queues with light-tailed arrival processes suffices (as developed for example in [START_REF] Duffield | Large deviations and overflow probabilities for the general single server queue, with applications[END_REF]).

Another strand of research [START_REF] Borst | Fluid Queues with Heavy-Tailed M/G/∞ input[END_REF]- [START_REF] Likhanov | Bounds on the buffer occupancy probability with self-similar input traffic[END_REF] constitutes the analysis of fluid queues fed by M/G/∞ input. Each active session causes the queue level to rise at a fixed rate. A number of observations can be made from these papers; the first being a confirmation of the fact that light-tailed session lengths lead to exponential tails while heavy-tailed service times lead to subexponential tails. In [START_REF] Likhanov | Bounds on the buffer occupancy probability with self-similar input traffic[END_REF], upper and lower bounds are obtained, leading to exact asymptotics for Pareto-distributed session lengths if a certain peak-rate condition is fulfilled.

In [START_REF] Borst | Fluid Queues with Heavy-Tailed M/G/∞ input[END_REF], the more general but interesting situation of several heterogeneous M/G/∞-processes is considered, each with different session rates and session durations, some of which are heavy-tailed. The typical configuration of long sessions that causes overflow (the path to overflow) is identified through an integer program, which paves the way for the exact asymptotics of the workload behavior.