C Briat

H Hjalmarsson

K H Johansson

U T Jönsson

G Karlsson

H Sandberg

Nonlinear state-dependent delay modeling and stability analysis of internet congestion control

It is shown that the queuing delay involved in the congestion control algorithm is state-dependent and does not depend on the current time. Then, using an accurate formulation for buffers, networks with arbitrary topologies can be built. At equilibrium, our model reduces to the widely used setup by Paganini et al. Using this model, the delay-derivative is analyzed and it is proved that the delay time-derivative does not exceed one for the considered topologies. It is then shown that the considered congestion control algorithm globally stabilizes a delay-free single buffer network. Finally, using a specific linearization result for systems with state-dependent delays from Cooke and Huang, we show the local stability of the single bottleneck network.

I. INTRODUCTION

Network congestion is an important problem inducing packet loss and large communication delays, which limit the performance of the overall network. The Transfer Control Protocol (TCP) used in most network adapters/servers includes a congestion control mechanism but, despite of this, it is generally recognized that it will be one of the bottlenecks as network capacity and heterogeneity in traffic behavior are scaled up. This has spurred a wide variety of efforts to design improved congestion control algorithms, e.g. [START_REF] Brakmo | TCP Vegas: end-to-end congestion avoidance on a global Internet[END_REF], [START_REF] Kelly | Rate control in communication networks: shadow prices, proportional fairness and stability[END_REF], [START_REF] Jin | FAST TCP: motivation, architecture, algorithms, performance[END_REF], [START_REF] Paganini | Congestion control for high performance, stability, and fairness in general networks[END_REF].

Design and analysis of such algorithms require accurate models. Here there has been a transition from packet-level models [START_REF] Jacobson | Congestion avoidance and control[END_REF] towards aggregated fluid flow models [START_REF] Kelly | Rate control in communication networks: shadow prices, proportional fairness and stability[END_REF] which fit better into the control theoretic and optimization frameworks. Recently, there has been significant progress in integrating these two model types. In [START_REF] Jacobsson | An improved link model for window flow control and its application to FAST TCP[END_REF] it is shown that it is important to include the self-clocking property of current window based TCP, i.e. the fact that the rate of transmission is controlled, or 'clocked', by the received acknowledgements, and a new link model is proposed. This is further developed in [START_REF] Tang | Window flow control: Macroscopic properties from microscopic factors[END_REF], [START_REF] Jacobsson | ACK-clock dynamics: Modeling the interaction between ACK-clock and network[END_REF], [START_REF] Andrew | Queue dynamics with window flow control[END_REF] leading to a model that is able to capture, e.g. sub-RTT burstiness and link interactions with great accuracy. The resulting model is shown to be able to predict previously unknown network behavior, e.g. that inelastic cross-traffic can severely deteriorate the network response and that FAST TCP may become unstable in some situations when there are heterogeneous round-trip-times. See [START_REF] Jacobsson | Dynamic modeling of internet congestion control[END_REF], [START_REF] Möller | Window-based congestion control[END_REF] for details.

The first of our contributions is to continue the development of more accurate network models. In Section II we develop a new model for general networks using a delay-based protocol. It is shown that the delay operator acting on the measurement (the queuing delay value) is not the delay at time t but the delay at time before t. This points out a recurrent misconception in the literature and is one of the main contributions of the paper. Then, a more accurate model for FIFO buffers is proved allowing for an explicit expression of all the output flows separately. Such a model can already be found in the literature (see e.g. [START_REF] Ros | Loss characterization in high-speed networks through simulation of fluid models[END_REF]) for simulating purpose only, in general without any theoretical justification. Finally, it is shown that using the general model it is, in principle, possible to describe any network topologies. When the network is at equilibrium, the model reduces to a standard model identical to the one in [START_REF] Pagaganini | A control theoretic look at internet congestion control[END_REF], suggesting that the proposed model is better disposed to describe the networks transient behavior.

The queuing delays imply that networked systems possess state dependent time-delays. This fact, together with nonlinear properties makes networked systems hard to analyze, restricting available results. We make several contributions in this respect. Firstly, Section III details an analysis of the particular time-delay involved in our new model. Indeed, thanks to the accurate buffer modeling, it is possible to provide an explicit expression for the queuing delays in the multiple bottleneck problem, allowing then for an analysis of the queuing delays in a more general setting than the single bottleneck problem. We show that, under general conditions, the delay derivative is lower than 1, permitting the use of many tools for time-delay systems available in the literature [START_REF] Niculescu | Delay effects on stability. A robust control approach[END_REF], [START_REF] Gu | Stability of Time-Delay Systems[END_REF].

Further, Section IV details some stability analysis results under various assumptions. The state dependent character makes stability analysis difficult. For example, existing stability results are local, where linearization is performed around the equilibrium time-delay. The latter has not been justified in the literature on stability of networked systems but by using a result in [START_REF] Cooke | On the problem of linearization for statedependent delay differential equations[END_REF], we show that this is correct. Furthermore, we show that the FAST-TCP protocol globally asymptotically stabilizes the nonlinear delay-free single bottleneck problem. We also show that the single bottleneck problem with one user can be made locally asymptotically stable for any network configuration provided that the protocol parameters are chosen accordingly.

II. NETWORK MODEL

This section is devoted to the definition of a network model. The first sections are devoted to the presentation of the models used to described the different elements of the network, namely the buffers and the users protocol. A particular accent is put on the modeling of the queuing delay. Finally, a single-buffer network model with multiple users is given with a discussion on the extension to multiple buffers networks.

A. Buffer Model

A buffer with infinite capacity can be modeled as

τ (t) = Q(φ(t), τ (t)) (1)
where τ (•) and φ(•) are the queuing delay and the input flow respectively. The mapping Q is given by

Q(φ(t), τ (t)) := 1 c (φ(t) -r(φ(t), τ (t))) (2a) r(φ(t), τ (t)) := c if C(t) φ(t) otherwise (2b)
with condition C(t) = [φ(t) ≥ c or τ (t) > 0] and where r(•) and c denote the output flow and the maximal output capacity respectively. Note that the buffer model ( 1)-( 2) is a differential equation with discontinuous right-hand side [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF].

In the sequel, we will assume that the latter model is wellposed and admits solutions.

B. Protocol Model

We next derive a protocol model, which defines how the congestion control mechanism interacts with the buffer model above. Consider a user sending data over a single buffer. The data flow at time t experiences a queuing delay τ (t). At time t + τ (t), the data has been served by the buffer. We suppose that the data are then instantaneously acknowledged to the sender. Hence, we assume that there is no propagation delay and that the round-trip time (RTT) of the congestion control loop, as seen by the sender at time t, is equal to τ (t).

τ (t s ) t s = t r -τ (g(t r )) t r = f (t s ) = t s + τ (t s ) Buffer r(φ, τ ) φ τ (g(t r ))
Fig. 1. Mappings illustration and reference times Fig. 1 describes the relations between the sending time t s and the reception time t r through the map

f : t → t + τ (t) (3)
allowing to compute the reception time t r from the sending time t s . This function is a forward operator, predicting the future from the present. Conversely, it will be shown that there also exists a mapping g : t → g(t) allowing for the computation of the sending time t s from the reception time t r . This operator is a backward one, allowing to retrieve the past from the present. Moreover, these operators are related to each other through the equalities

t r = t s + τ (t s ) (4a) = g(t r ) + τ (g(t r )) (4b) t s = g(t r ) (4c)
as illustrated in Fig. 1. Using these relations, it can be deduced that the RTT satisfies

RT T = t r -t s = τ (t s ) = t r -g(t r ) = τ (g(t r )) (5)
The interest of the operator g(•) lies in its ability to retrieve the time instant where the data entered the buffer from the knowledge of the acknowledgment reception time only.

The next lemma provides further results on f and g: Lemma 2.1: The function g exists provided that the input flow φ(t) > 0 for almost all t ∈ R + . It coincides with the inverse function of the mapping f : t → t + τ (t) where τ evolves according to (1)- [START_REF] Kelly | Rate control in communication networks: shadow prices, proportional fairness and stability[END_REF]. Moreover, we have the following relations:

D + [g] = 1 1 + D + [τ ](g) (6a) = c/φ(g) if C(g) 1 otherwise (6b) D + [f ] = φ/c if C 1 otherwise (6c)
where

D + [f ](t) = lim sup h↓0 f (t + h) -f (t) h
is the upper right Dini derivative.

Proof: We will prove here Lemma 2.1 as well as relations (4). Since t r = f (t s ), if we seek a function g satisfying t s = g(t r ), this means that g corresponds to the inverse function of f . Since by definition t and τ (t) are continuous, the invertibility of f is equivalent to its monotonicity on R + . The function τ is differentiable almost everywhere, so the Dini derivative is considered instead of the usual derivative:

D + [f ] = φ/c if C 1 otherwise (7)
Since c > 0 and φ ≥ 0 then D + [f ] ≥ 0. Hence the only condition for f to be monotonic is φ = 0 except, possibly, at a countable number of points. Under the weak assumption φ > 0 almost everywhere (there exists a communication using the buffer) we can conclude on the global existence of the bijective inverse function g = f -1 almost everywhere. Equations (4a)-(4b) and (4c) are direct applications of the relations f • g = g • f = Id where Id is the identity operator.

To prove (6b), let us consider equation (4b) evaluated at t. Applying the Dini derivative yields

1 = D + [g] + D + [g]D + [τ ](g) (8)
or equivalently

D + [g] = 1 1 + D + [τ ](g) (9)
Finally substitute (1) inside the latter equation to get (6b).

Remark 2.1:

The possibility of retrieving t s from t r may seem strange at first sight. Actually, this is an immediate consequence of the behavior of the passive FIFO buffer.

F 1 (s) F 2 (s) F 12 (s) B 1 (s) B 2 (s) B 12 (s) e -T s q 12 (t) y 1 (t) y 2 (t) y 12 (t) ỹ12 (t) Fig. 2

. Example of a multiple bottleneck network model

Each data is stored in order of arrival and this order is preserved when data are served. The preservation of the order is responsible of the monotonicity of f . Lemma 2.2: The delay measured by the protocol at time t is given by τ (g(t)). Moreover, the delay acting on the state τ is also equal to τ (g(t)).

Proof: When an acknowledgment is received at time t r the information on the RTT contained in the acknowledgment is τ (t s ) and is thus a delayed value of the queuing delay since t s ≤ t r . Using (4c) and (4b) we can conclude that the measured queuing delay at time t r is equal to τ (g(t r )) = τ (t r -τ (g(t r ))), thus normalizing the time as t = t r concludes the proof. Using the above results, the following generic model for a delay-based protocol can be given:

żi (t) = P i (z i (t), τ (g(t))) x i (t) = h i (z i (t)) (10)
where z i , x i are the state of the protocol and the number of outstanding packets for user i respectively. They are related through the output function h i . The measured delay τ (g(t)) has been defined in Lemma 2.2. Remark 2.2: The model ( 10) is a backward model of the protocol where the reference of time is the reception time of ACK's. It is also possible to define a forward expression using the sending time as a reference. However, even if the latter is more intuitive, the former allows to analyze properties of the delay more easily.

Remark 2.3: The model ( 10) belongs to the class of nonlinear systems with state-dependent delays with few results on the stability in the control community [START_REF] Cooke | On the problem of linearization for statedependent delay differential equations[END_REF], [START_REF] Verriest | Stability of systems with state-dependent and random delays[END_REF], [START_REF] Walther | A periodic solution of a differential equation with state-dependent delay[END_REF], [START_REF] Györi | Exponential stability of a state-dependent delay system[END_REF]. The analysis of such systems is reputed to be difficult due to the strong nonlinear nature of the state-dependent delay.

C. Single Bottleneck Network Model

Gathering the buffer model ( 1)-( 2) and the protocol model [START_REF] Jacobsson | Dynamic modeling of internet congestion control[END_REF], the closed-loop system with multiple users is governed by:

τ (t) = Q [ i φ i (t), τ (t)] φ i (t) = Φ x i t -T f i , τ (t) + T i żi (t) = P i z i (t), τ t -T b i -τ g t -T b i x i (t) = h i (x i (t)) (11)
where T f i , T b i are respectively the forward and backward delay for user i and

T i = T f i + T b i .
Since the protocol output is a number of outstanding packets to maintain in the network and the buffer input consists of flows, it is necessary to define an operator Φ(•) converting the number of outstanding packets into flows. It has been shown in [START_REF] Jacobsson | An improved link model for window flow control and its application to FAST TCP[END_REF] that the following operator:

x i (t) = t g(t-T b i )-T f i φ i (s)ds (12)
relates the number of outstanding x i packets from a given flow φ i . However, the integral operator [START_REF] Ros | Loss characterization in high-speed networks through simulation of fluid models[END_REF] defines the flow for a given number of outstanding packets. Thus Φ(•) coincides with the inverse operator of [START_REF] Ros | Loss characterization in high-speed networks through simulation of fluid models[END_REF]. It is not possible to describe it explicitly, thus approximations are generally used in the literature. Examples can be found in [START_REF] Jacobsson | An improved link model for window flow control and its application to FAST TCP[END_REF], [START_REF] Jacobsson | Dynamic modeling of internet congestion control[END_REF] where it has been shown that the use of ( 12) provides accurate results in simulations and experiments.

D. Extension to Multiple Bottleneck Network Model

Let us consider the topology considered in Fig. 2 at the top of the page, where three users share two buffers. Users 1 and 2 use buffers 1 and 2 respectively while user 12 uses both of them. The blocks F i and B i , i ∈ {1, 2, 12} denote respectively forward and backward constant propagation delays.

Suppose now we would like to derive a network model to describe this topology. It turns out that it is not possible to use the model ( 1)-( 2) since the modeling does not allow for an explicit expression of the output flow q 12 corresponding to the input flow y 12 of buffer 1. Indeed, (1)-(2) only defines an aggregate output flow r(•, •) for all input flows.

The problem consisting in finding an explicit formulation for the output flow r(φ, τ ) is solved in the following proposition:

Proposition 2.1: The output flow corresponding to input flow φ i (t) in model ( 1) is given by

r i (φ, τ ) =    φ i (g(t))c j φ j (g(t)) if C(g(t)) φ i (t) otherwise (13)
Moreover, the aggregate output flow is given by r = i r i and satisfies (2b).

Proof: The starting point is equation [START_REF] Ros | Loss characterization in high-speed networks through simulation of fluid models[END_REF]. First, assume w.l.o.g. that propagation delays are neglected and define the initial time to be g(t) with initial condition τ (g(t)) > 0. Since the time needed for the buffer to serve the data initially contained at g(t) is equal to the initial condition τ (g(t)), this implies that the data sent by the user during the interval [g(t), t] are exactly those stored in the buffer at time t. To analyze the output flow corresponding to a given input flow, we will virtually separate the contribution of each input flow to the value of the delay. This is possible using the superposition principle since the flows enter linearly in the buffer model (2a). So, defining the contribution of the flow i by τ i we have:

τ (t) = i τ i (t) with τ i (t) = 1 c t g(t) φ i (s)ds (14)
and, assuming a continuous flow, differentiation yields

τi (t) = 1 c [φ i (t) -φ i (g(t))D + [g](t)] =    1 c φ i (t) - φ i (g(t)) j φ j (g(t)) if C(g(t)) 0 otherwise (15)
We can see that the total output flow r(•) can be decomposed into the sum of atomic r i (•) defined by [START_REF] Pagaganini | A control theoretic look at internet congestion control[END_REF].

Remark 2.4: Not surprisingly, the output flow (13) consists of a delayed and scaled version of the corresponding input flow. Indeed, the proportion of each input flow is preserved and the buffer hence plays the role of a delay operator on the proportion of flows, not on the flows themselves.

Remark 2.5: Using this model, it is in principle possible to build any network topology, such as the one in Fig. 2. Moreover, when the model is at equilibrium, its behavior can be well represented by an LTI system. Propagation delays can then commute with buffers and buffers can be put in parallel using appropriate forward and backward routing matrices. Such a transformation results in the analysis setup studied in [START_REF] Pagaganini | A control theoretic look at internet congestion control[END_REF].

III. DELAY CHARACTERIZATION

This section aims at studying important intrinsic properties of the delay. It will be shown that in the single bottleneck problem, the delay derivative cannot exceed one. However, this nice property may not be preserved when more general topologies are considered.

A. Queuing Delay Analysis

This section is devoted to the analysis of the characteristics of the queuing delay τ (g(t)) in the single bottleneck problem, i.e. model [START_REF] Möller | Window-based congestion control[END_REF]. Note that the presence of the term g(t)

does not change the delay values range. However, its derivative will be strongly affected by this time-scale. It is well known that fast varying delays (whose derivative exceeds 1) are responsible of sudden loss of stability (and also violate causality [START_REF] Verriest | Causal behavior of switched delay systems as multi-mode multi-dimensional systems[END_REF]) compared to their time-invariant or slowlyvarying counterparts [START_REF] Gu | Stability of Time-Delay Systems[END_REF], [START_REF] Kharitonov | On the stability of systems with uncertain delay[END_REF], [START_REF] Papachristodoulou | Stability analysis of linear systems with time-varying delays: Delay uncertainty and quenching[END_REF]. A fast decreasing delay is not penalizing at all, so attention should be paid on the upper bound of the delay derivative only.

Since the queuing delay involved in the protocol ( 10) is equal to τ (g(t)), the analysis of the delay derivative dτ (g(t))/dt is less straightforward. Such a derivative analyzes variations of the delay at sending times with respect to variations of the reception time.

Theorem 3.1: The delay derivative obeys

D + [τ (g)](t) = 1 - c φ(g(t)) if C(g(t)) 0 otherwise (16)
and is bounded from above by one.

Proof: The Dini derivative of τ (g(t)) is given by

D + [τ (g)](t) = D + [τ ](g(t)) • D + [g](t) (17)
Using (6b), we get

D + [τ (g)](t) = D + [τ ](g(t)) 1 + D + [τ ](g)(t) (18)
Substituting the explicit expression of D + [τ ](g(t)) given by (2) evaluated at g(t) yields ( 16). The delay derivatives D + [τ ](t s ) and D + [τ (g)](t r ) are depicted in Fig. 3 where we can see the large difference between their behavior. Indeed, the delay derivative from the sending viewpoint dτ (t)/dt exceeds one when the input flow exceeds the double of the buffer maximal output capacity and hence can be made arbitrarily large. Moreover, τ (t) cannot decrease faster than a unitary slope.

On the other hand, the delay derivative viewed from the reception viewpoint dτ (g(t))/dt is bounded from above by one and is allowed to decrease arbitrarily fast. This illustrates the large qualitative difference between the delay behavior and shows that considering dτ (t)/dt instead of dτ (g(t))/dt is, in fact, very misleading.

B. Delay-derivative for Interconnected Buffers

In this section, we emphasize an important difficulty related to delay derivative in the multiple bottleneck problem. Indeed, in this case, there exist configurations for which the delay derivative can theoretically exceed one. To illustrate this, let us consider the Fig. 2 and focus on the price received by user 12. Users 1 and 2 are considered as unregulated cross-traffic (disturbances) and we ignore propagation delays. Using similar arguments as in Section III-A we get: where the g i 's are the inverse functions of f i : t → t + τ i (t), i = 1, 2. Thus the RTT can be expressed as:

t r = t s + τ 1 (t s ) + τ 2 (t s + τ 1 (t s )) (20a) t s = g 1 (g 2 (t r )) (20b) D + [τ 1 (g 1 (g 2 ))] =          (φ 1 (g 1 (g 2 )) -c 1 )c 2 φ 1 (g 1 (g 2 ))y 2 (g 2 ) + y 12 (g 1 (g 2 ))c 1 if C 1 (g 1 (g 2 )) and C 2 (g 2 ) 1 - c 1 φ 1 (g 1 (g 2 )) if C 1 (g 1 (g 2 )) and not C 2 (g 2 ) 0 otherwise (19) replacements User Packet φ(s) 2c 3c/2 c ր +∞ ↓ -∞ → 1 c/2 1/2 -1/2 1 -1 0 Delay Derivatives D + [τ ](t s ) D + [τ (g)](t r )
RT T = τ 1 (t s ) + τ 2 (t s + τ 1 (t s )) = τ 1 (t r -τ 1 (g 1 (g 2 (t r )))) +τ 2 (t r -τ 2 (g 2 (t r ))). (21)
The obtained closed-loop system will contain the two delayed terms. It has been shown in Section III-A that the derivative of τ 2 (g 2 (t)) does not exceed 1. So, we focus on the delay τ 1 (g 1 (g 2 )) for which we have the following theorem: Theorem 3.2: Consider Fig.

2.

The derivative of τ 1 (g 1 (g 2 )) is governed by model [START_REF] Walther | A periodic solution of a differential equation with state-dependent delay[END_REF] where

C i (α) = [φ i (α) > c i or τ i (α) > 0].
The aggregate input flows of buffers 1 and 2 are given by φ 1 = y 1 + y 12 and φ 2 = ỹ12 + y 2 respectively with

ỹ12 (g 2 ) =    y 12 (g 1 (g 2 ))c 1 φ 1 (g 1 (g 2 )) if C 1 (g 1 (g 2 ))
y 12 (g 2 ) otherwise [START_REF] Kharitonov | On the stability of systems with uncertain delay[END_REF] Proof: The proof is given in Appendix A. Proposition 3.1: If c 1 ≥ c 2 , causality is always preserved. However when c 2 > c 1 , with zero cross-traffic (i.e. y 1 = y 2 = 0), the causality is preserved if and only if the input flow y 12 satisfies

y 12 < c 1 + c 2 1 c 2 -c 1 . ( 23
)
Proof: This is an immediate conclusion drawn from the model [START_REF] Walther | A periodic solution of a differential equation with state-dependent delay[END_REF]. The latter result states that under some mild assumptions (the input flow y 12 may exceed the output capacity c 1 by c 2 1 /(c 2 -c 1 )), the causality is preserved. However, when the difference c 2 -c 1 is large, the assumption becomes stronger and can be more easily violated. In such a case, according to some recent results, the notion of time-delay system collapses and the model might be no longer valid, but the network can still have a correct behavior. Due to space limitation, it is not possible to provide a detailed and general analysis of the delay behavior in presence of cross-traffic. Indeed, a general and relevant analysis is more involved and is left to a future paper.

IV. STABILITY ANALYSIS

This section is devoted to the stability analysis of the congestion control. It will be shown that the congestion control mechanism leads to a globally asymptotically stable network in the delay-free single bottleneck problem (Theorem 4.2). Then, using a linearization result for systems with statedependent delays [START_REF] Cooke | On the problem of linearization for statedependent delay differential equations[END_REF], local stability of the network is proved and stated through Theorem 4.3, Lemma 4.1 and Theorem 4.4. These latter results are the main results of this section on stability analysis.

A. Linearization of the nonlinear model

The nonlinear model [START_REF] Möller | Window-based congestion control[END_REF] can be difficult to analyze depending on the choice of the flow operator Φ. To the authors' best knowledge no general global stability results exist for delay-based protocol (such as FAST-TCP [START_REF] Wei | FAST TCP: Motivation, architecture, algorithms, performance[END_REF]) taking into account all the involved phenomena. Some results can be found in [START_REF] Pagaganini | A control theoretic look at internet congestion control[END_REF], [START_REF] Choi | Global exponential stability of FAST TCP[END_REF], [START_REF] Lestas | Scalable robust stability for nonsymmetric heterogeneous networks[END_REF], [START_REF] Jacobsson | An improved link model for window flow control and its application to FAST TCP[END_REF]. The main reason for these difficulties comes from the presence of nonlinearities, statedependent delays, large number of heterogeneous constant delays, hybrid behavior and an arbitrarily large dimension; making the overall system analysis very intricate.

A usual simplification relies on the linearization of nonlinear model but since the system involves a state-dependent delay, the linearization seems to be a difficult task. Fortunately, the following fundamental result, first provided in [START_REF] Cooke | On the problem of linearization for statedependent delay differential equations[END_REF], restated below, explains how the linearization can be performed.

Theorem 4.1 ([16]): Let us consider a continuous function f (x, y) mapping R n × R n to R n . The nonlinear system

ẋ(t) = f (x(t), x(t -h(x(t))) (24)
with continuous state-dependent delay h(x) ≥ 0 and equilibrium point x * is locally asymptotically stable if the following LTI time-delay system

ż(t) = A 0 z(t) + A 1 z(t -h(x * )) (25)
is asymptotically stable with

A 0 = ∂f ∂x (x * , x * ) and A 1 = ∂f ∂y (x * , x * ).
Surprisingly, no derivative of h(•) is involved. This means that the shape of h(•) around the equilibrium point has no influence on the local stability, only the delay value at equilibrium is important. So, this theorem turns an apparently difficult local stability analysis problem into a very simple one. Indeed, many results are available to study exactly the stability of LTI systems with delays [START_REF] Niculescu | Delay effects on stability. A robust control approach[END_REF], [START_REF] Gu | Stability of Time-Delay Systems[END_REF], [START_REF] Michiels | Stability and stabilization of time-delay systems. An eigenvalue based approach[END_REF]. The latter theorem justifies the commonly considered stability analysis problem [START_REF] Pagaganini | A control theoretic look at internet congestion control[END_REF] where the network is at equilibrium and the queuing delays are set constant.

B. Global Stability Analysis of the Delay-Free Model

In this subsection, we will consider the problem of stability analysis of the model [START_REF] Möller | Window-based congestion control[END_REF] where Φ : x → x/RT T [START_REF] Jacobsson | An improved link model for window flow control and its application to FAST TCP[END_REF]:

τ (t) = 1 c i x i (t -T f i ) τ (t) + T i + δ(t) -1 (26)
where δ(t) ∈ [0, 1) is the scaled cross-traffic. The second mode of the buffer is not recalled here but will be considered in the analysis. We will also consider the following model for the FAST-TCP protocol:

ẋi (t) = γ - τ (g(t -T b i )) T i + τ (g(t -T b i )) x i (t) + α i (27)
where γ, α i > 0 are a tuning constant term and the reference of the protocol respectively. The total constant propagation delay for user i is denoted by

T i = T f i + T b i . Theorem 4.2: The unique equilibrium τ * = α s c(1 -δ * ) , x * i = α i 1 + T i τ * (28)
of the nonlinear delay-free system

τ (t) = 1 c i x i (t) τ (t) + T i + δ(t) -1 ẋi (t) = γ - τ (t) T i + τ (t) x i (t) + α i (29)
is globally asymptotically stable for any α i ≥ 0, α s = i α i , c, γ > 0 and equilibrium cross-traffic δ * ∈ [0, 1).

Proof: The proof is given in Appendix B.

C. Local stability analysis in the single-user problem

1) Derivation of the Linear Model: Using Theorem 4.1 we can state that the linearized model of ( 26)-( 27) around the equilibrium point (28) is given by

ż(t) = µ -ξ 1 0 ξ 2 -ξ 3 A z(t) + µ 2 -ψ 1 0 A h ν(t -τ * -T ) (30) where z(t) = col(y(t), ν(t)), y = x -x * , ν = τ -τ * , µ = (τ * + T ) -1 , ξ 1 = γτ * , ξ 2 = c -1 , ξ 3 = α(cτ * ) -1 and ψ 1 = γT x * .
2) Delay Independent Stability of the Linear Model: This section is devoted to the analysis of the delay-independent stability of the system

ż(t) = Az(t) + A h ν(t -h) (31)
where h ≥ 0 is an unknown constant delay, A and A h are defined in (30).

The following result states a result on the delay independent stability of (31):

Theorem 4.3 (Delay-Independent Stability): The system (31) is delay independent stable if and only if the condition T /τ * = c(1 -δ * )T /α < 1 holds.

Proof: Since the delay acts only on ν, it is thus possible to rewrite the system as a loop interconnection of the SISO transfer function H(s) = C(sI -A) -1 A h with C = 0 1 and the delay operator e -sh . From [START_REF] Gu | Stability of Time-Delay Systems[END_REF], it is possible to conclude that the system is delay independent stable if and only if both A and A + A h C are Hurwitz and ||H|| H∞ < 1 1 . It is straightforward to show using Routh-Hurwitz criterion that both A and A + A h C are Hurwitz. Moreover, it is also easy to show that ||H|| H∞ = H(0) = T /τ * . This concludes the proof. From this result, we can conclude that as long as the queuing delay at equilibrium is larger than the propagation delay, the system ( 26)-( 27) will be asymptotically stable. However, if this condition is not satisfied, then the system (31) becomes delay-dependent stable.

3) Delay Dependent Stability of the Linear Model: As discussed in the latter section, the system can be delaydependent stable in some specific scenarios, i.e. there exists a positive scalar h max > 0 such that the system (31) is unstable for all h ≥ h max . This is stated in the following lemma where a conservative (lower bound) on the critical delay value is provided: Proof: Indeed, since we know that the system is stable for h = 0 and h = ∞, this means that there exists a h max (x * , τ * ) > 0 (depending on the linearization point) for which the system is stable for all constant h ∈ [0, h max ). To prove this, let us consider the polynomial a(s, w) = det(sI -A -A h Cw) and look for zero-crossings of the quasipolynomial a(jω, e -jωh ) = 0 which are clues for sign changes of zeros of a(•, •). This yields

-ω 2 + µ 2 ξ 1 ξ 3 + µ 3 ξ 2 ψ 1 cos(ωh) = 0 (32a) µ(ξ 1 + ξ 3 )ω -µ 3 ξ 2 ψ 1 sin(ωh) = 0 (32b)
after separation of the real and imaginary parts. Let us focus on (32b), which rewrites

ξ 1 + ξ 3 hµ 2 ξ 2 ψ 1 = sin(ωh) ωh , ω > 0.
Since the right-hand side is bounded from above by one, then by choosing a sufficiently small h > 0, the left-hand side can be made greater than one. In such a case, equality (32b) cannot be satisfied and hence there is no ω > 0 such that (32) holds. Thus choosing a sufficiently small h prevents any zero-crossing and stability of the system is guaranteed. It is sufficient to choose h > 0 which satisfies

ξ 1 + ξ 3 µ 2 hξ 2 ψ 1 > 1
or equivalently:

h < τ * (1 -δ * )T + 1
γT .

Thus the system will be stable for all h ∈ [0, h max ) with

h max = τ * (1 -δ * )T + 1 T γ .

4) Local Stability of the Network:

We have shown that the system (31) is delay-dependent stable for any values for the linearization point. Moreover, we have derived a conservative lower bound of the exact delay-margin. We can provide then the main result on stability analysis: Theorem 4.4 (Network Local Stability): Any equilibrium point (x * , τ * ) of the single-user nonlinear delay-system ( 26)-( 27) is:

1) locally delay-independent stable if and only if τ * > T .

2) locally delay-dependent stable if τ * < T and τ * (T -1) + T 2 ≤ 0. Proof: Since we know the exact delay value τ * at each linearization point, it is sufficient to show that τ * + T is contained in [0, h max ). We already know that if τ * > T the system is stable, so the case τ * ≤ T remains to be studied and two cases may occur. The first one is when τ * (T -1) + T 2 ≤ 0 and in this case the h max is always greater than τ * + T independently of the value of γ. Finally, when τ * (T -1) + T 2 > 1, then in order to have h max > τ * + T , we need that γ < 1 τ * (T -1) + T 2 .

V. CONCLUSION

In this paper, a precise modeling for networks with multiple buffers has been provided with a particular accent on the queuing delays. It has been shown that the exact measured queuing delay is delayed by a function depending on the buffer state, and thus the network could be represented as a system with state-dependent delays. Moreover, this modeling has allowed for a precise analysis of the derivative of the queuing delay showing that, in the single bottleneck problem, this delay is always bounded from above by one. However, this property might not be satisfied for more complex topologies. A second part has been devoted to the stability analysis of a network using a FAST-TCP-like algorithm. It has been shown, using a fundamental linearization result for systems with state-dependent delays from Cooke and Huang [START_REF] Cooke | On the problem of linearization for statedependent delay differential equations[END_REF], that a single bottle-neck single-user network is locally asymptotically stable at any equilibrium point provided tuning parameters are chosen accordingly.

Future works will be devoted to a more detailed analysis of the queuing delays and on nonlinear stability analysis of the network in order to obtain global stability results or estimate the stability region.

Fig. 3 .

 3 Fig. 3. Delay derivative viewed at sending (plain) and receiving (dashed).

Lemma 4 . 1 (

 41 Delay-Dependent Stability): The system (31) is delay-dependent stable for any constant delay h < h max with h max < τ * (1 -δ * )T + 1 T γ and for any values c, γ > 0, α ≥ 0 and δ * ∈ [0, 1).

3 )

 3 locally delay-dependent stable if τ * < T , τ * (T -1) + T 2 > 0 and γ < 1 τ * (T -1) + T 2 .

This is a necessary condition since SISO systems are involved. In general, these conditions are only sufficient.

A. Proof of Theorem 3.2 The RTT equation [START_REF] Verriest | Causal behavior of switched delay systems as multi-mode multi-dimensional systems[END_REF] shows that two different delays are involved in the expression of the RTT, namely τ 2 (g 2 ) and τ 1 (g 1 (g 2 )). The former has the same structure as the one analyzed in Section III-A, and we know that its derivative is bounded from above by 1, so let us focus on the second one. Its Dini derivative is given by

Similarly as in Section III-A, the derivative of D + [τ 1 (g 1 )](g 2 ) can be computed and is equal to

is the total input flow of buffer 1. On the other hand, we have

where φ 2 (•) = y 2 (•) + ỹ12 (•) is the total input flow of buffer 2. However, in order to take explicitly into account the flow circulating between the buffers, we need to provide an exact expression for ỹ12 . Such an exact formulation has been provided in the Proposition 2.1 in Section II-D. Using this result, we get [START_REF] Kharitonov | On the stability of systems with uncertain delay[END_REF]. Finally evaluating (33) using (34), ( 35) and ( 22) yields [START_REF] Walther | A periodic solution of a differential equation with state-dependent delay[END_REF]. This concludes the proof.

B. Proof of Theorem 4.2

For given α i , c > 0 and δ * ∈ [0, 1), the unique equilibrium point can be determined by solving for the stationary solutions of (29). Shifting the model (29) around the equilibrium with the change of variable

where we have set ε = 0. Note that this model is valid only when the buffer is non empty or when the input flow exceeds the effective maximal output capacity of the buffer. In the other case, the equation ν(t) = 0 must be used instead. We will study the stability of the origin of the latter system. To this aim define the Lyapunov function

where p ν > 0 is scalar. Computing the derivative of V along trajectories solution of (36) we get the following expression

Thus when ν > -τ * the Lyapunov function is exponentially decreasing. However, this does not say anything when ν = -τ * . Indeed, in such a case, we have V (t) = -p ν i α i τ * /T i < 0 which only says that V is strictly decreasing, but we have no exact information on how the states evolve. Moreover, when the input flow i y i is lower than (1 -δ * )c, the system is in the noncongested mode which is not taken into account in the inequality (37). In order to characterize the stability of the hybrid system over the whole state-space (y, ν) ∈ [-x * , +∞) × [-τ * , +∞), we will differentiate three cases:

i and ν = -τ * (empty buffer with input flow lower than the effective maximal output capacity

i and ν = -τ * (empty buffer with input flow greater than the effective maximal output capacity (1 -δ * )c) 3) ν > -τ * (non empty buffer); the exponential stability has been proved above. Case 1: This case describes the case when the aggregate flow is small (lower than the maximal output capacity of the buffer) and the buffer is empty. This configuration typically could arise when a network is 'starting'. Considering the network equations with initial conditions ν(t 0 ) = -τ * and

for t ≥ t 0 . Thus after some time t 1 > t 0 , we have i [y i (t 1 ) + x * i ] = (1 -δ * )c and we switch to case 2. Case 2: The model is given by (36) with initial conditions (y(t 1 ), ν(t 1 )) = (y(t 1 ), -τ * ) for t ≥ t 1 . Note that ν(t 1 ) = -τ * , so it cannot decrease more and note also that y i (t) > y i (t 1 ) for all t > t 1 by virtue of the first equation ẏi (t) = γx * i τ * . Since the contribution of the y i to the ν dynamics is positive (i.e. ν increases when y i increases), then it will exist a time t 2 ≥ t 1 for which ν(t 2 ) > 0 and we switch to case 3 for which stability has been proved. This concludes the proof.


Z Mladenova

N Korsun

P Petrov

K Peteva

T Geonova

M Iturriza-Gómara

Molecular epidemiology of rotaviruses in Bulgaria: annual shift of the predominant genotype

Rotavirus molecular epidemiology investigations provide important information about the incidence of rotavirus diseases and rotavirus strains in circulation in the prevaccine era. The purpose of this investigation was to study the burden of rotavirus disease, rotavirus strain diversity, and epidemiology specificities of rotavirus infections in Bulgaria. A total of 3,130 stools collected between 2005 and 2008 were tested by immune enzyme tests. G-P genotype identification of rotavirus strains were performed by reverse transcriptase polymerase chain reaction (RT-PCR). Rotavirus etiology was confirmed in 32.4% of the samples tested. Rotaviruses affected predominantly children under 5 years of age (95.5%), with a peak prevalence between the ages of 7 and 36 months. Four of the five globally distributed rotavirus strains (G1P[8], G2P[4], G4P[8], and G9P[8]) constituted 97.7% of all rotavirus strains in circulation. However, annual shifts of predominant rotavirus G-P genotypes were observed from season to season-G4P[8] was predominant in rotavirus season 2004/2005 (56.8%), but was replaced by G9P[8] in 2005/2006 (77.7%), and G2P[4] (41.6%) and G1P[8] (39.5%) in the following two consecutive rotavirus seasons. Yearround circulation of rotaviruses in the country with increased incidence in the winter-spring season and unexpected peaks preceding the rotavirus seasons were observed. Molecular epidemiology data are needed in Bulgaria for health policy makers in order to introduce routine rotavirus vaccination. The monitoring of rotavirus genetic diversity in Bulgaria in the postvaccination period will contribute to a successful rotavirus vaccination program.

Introduction

Despite the success of the current public health services and hygiene control in water supply and sanitation, acute gastroenteritis remains the second most frequent infectious disease after respiratory infections all over the world [1,[START_REF] Bryce | WHO estimates of the causes of death in children[END_REF]. Rotaviruses are the most common etiologic agent of severe infantile diarrhea. Annually, rotaviruses are responsible for 140 million diarrheal episodes and 2 million hospitalizations and approximately 611,000 deaths worldwide, thus, it remains one of the most important challenges in public health [START_REF] Parashar | Rotavirus and severe childhood diarrhea[END_REF].

Rotaviruses are nonenveloped, triple-layered icosahedral particles, of ∼75 nm in diameter, comprising the genus Rotavirus in the Reoviridae family. The viral genome is composed of 11 segments of double-stranded RNA. The VP7 glycoprotein and the VP4 protease-sensitive protein form the outer capsid; these elicit neutralizing antibodies and form the basis for the dual classification system for rotaviruses into G and P types [START_REF] Estes | Rotaviruses[END_REF]. To date, 22 G and 30 P types have been reported in humans, mammals, and birds [START_REF] Estes | Rotaviruses[END_REF][START_REF] Trojnar | The first complete genome sequence of a chicken group A rotavirus indicates independent evolution of mammalian and avian strains[END_REF][START_REF] Schumann | Evidence of interspecies transmission and reassortment among avian group A rotaviruses[END_REF]. Among these, only 11 G serotypes/genotypes and 12 P serotypes and 15 P genotypes have been found in humans [START_REF] Kapikian | Rotaviruses[END_REF][START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF][START_REF] Martella | Identification of a novel VP4 genotype carried by a serotype G5 porcine rotavirus strain[END_REF]. Although the ability of rotavirus genes to segregate independently upon dual infection could, in theory, lead to 132 G-P combinations, five G and P combinations, G1P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF], G2P [START_REF] Estes | Rotaviruses[END_REF], G3P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF], G4P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF], and G9P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF], are epidemiologically important in human infections [START_REF] Gentsch | Serotype diversity and reassortment between human and animal rotavirus strains: implications for rotavirus vaccine programs[END_REF][START_REF] Santos | Global distribution of rotavirus serotypes/genotypes and its implication for the development and implementation of an effective rotavirus vaccine[END_REF]. However, geographical and temporal differences exist in the distribution of G and P types and their combinations: P [START_REF] Schumann | Evidence of interspecies transmission and reassortment among avian group A rotaviruses[END_REF] in association with G1, G2, G3, and G9 and G8 in combination with P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF] and P [START_REF] Schumann | Evidence of interspecies transmission and reassortment among avian group A rotaviruses[END_REF] have been shown to be the second most detected P and G type in India and some parts of the African continent, respectively [START_REF] Santos | Global distribution of rotavirus serotypes/genotypes and its implication for the development and implementation of an effective rotavirus vaccine[END_REF][START_REF] Ramachandran | Unusual diversity of human rotavirus G and P genotypes in India[END_REF][START_REF] Cunliffe | Rotavirus G and P types in children with acute diarrhea in Blantyre, Malawi, from 1997 to 1998: predominance of novel P[6]G8 strains[END_REF][START_REF] Chouikha | Group A rotavirus strains circulating in the eastern center of Tunisia during a tenyear period (1995-2004)[END_REF][START_REF] Benhafid | Epidemiology of rotavirus gastroenteritis among children <5 years of age in Morocco during 1 year of sentinel hospital surveillance, June 2006-May 2007[END_REF][START_REF] Esona | Genomic characterization of human rotavirus G8 strains from the African rotavirus network: relationship to animal rotaviruses[END_REF]. Recent studies have reported the emergence of uncommon G or P types (G5, G6, G10, G11, G12, P [START_REF] Parashar | Rotavirus and severe childhood diarrhea[END_REF], P [START_REF] Kapikian | Rotaviruses[END_REF], P [START_REF] Martella | Identification of a novel VP4 genotype carried by a serotype G5 porcine rotavirus strain[END_REF], P [START_REF] Santos | Global distribution of rotavirus serotypes/genotypes and its implication for the development and implementation of an effective rotavirus vaccine[END_REF], and many more) in different countries [START_REF] Gentsch | Serotype diversity and reassortment between human and animal rotavirus strains: implications for rotavirus vaccine programs[END_REF][START_REF] Santos | Global distribution of rotavirus serotypes/genotypes and its implication for the development and implementation of an effective rotavirus vaccine[END_REF][START_REF] Bányai | Eight-year survey of human rotavirus strains demonstrates circulation of unusual G and P types in Hungary[END_REF][18][19][START_REF] Esona | Novel human rotavirus genotype G5P[7] from child with diarrhea, Cameroon[END_REF][START_REF] Bányai | Genetic diversity and zoonotic potential of human rotavirus strains, 2003-2006, Hungary[END_REF][START_REF] Ramani | Whole genome characterization of reassortant G10P[11] strain (N155) from a neonate with symptomatic rotavirus infection: identification of genes of human and animal rotavirus origin[END_REF]. These uncommon rotavirus G/P genotypes emerge in humans through zoonotic transmission, and reassortment with common human rotavirus strains can potentially lead to the emergence of a pandemic strain [START_REF] Iturriza-Gómara | Reassortment in vivo: driving force for diversity of human rotavirus strains isolated in the United Kingdom between 1995 and 1999[END_REF].

In 2006, two rotavirus vaccines, a monovalent live attenuated G1P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF] vaccine (Rotarix ® , GlaxoSmithKline Biologicals) and a pentavalent human-bovine G1-G4,P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF] reassortant vaccine (Rotateq ® , Merck and Sanofi Pasteur MSD), were licensed. Two of the largest clinical trials in history showed that these vaccines are safe and highly effective for preventing life-threatening dehydrating rotavirus diarrhea, providing either serotype-specific or crossreactive protective immune responses [START_REF] Ruiz-Palacios | Safety and efficacy of an attenuated vaccine against severe rotavirus gastroenteritis[END_REF][START_REF] Vesikari | Safety and efficacy of a pentavalent human-bovine (WC3) reassortant rotavirus vaccine[END_REF].

The introduction of rotavirus vaccine into the national immunization programs of many countries have been accompanied by an estimation of the rotavirus disease burden and rotavirus strain diversity and an evaluation of changes over the rotavirus incidence and rotavirus genotypes in circulation in the postlicensure period.

In Bulgaria, early investigations of rotavirus infections were performed by using electron microscopy and polyacrylamide gel electrophoresis, and established a high incidence of rotavirus gastroenteritis in children <7 years of age: 7-11% among sporadic cases, up to 67% during two outbreaks in the Shumen region [START_REF] Shindarov | Five year study of rotavirus gastroenteritis in Bulgaria[END_REF][START_REF] Dimitrov | Distribution of rotaviruses in Bulgaria. In: Viral gastroenteritis[END_REF]. The studies also revealed that children <1 year of age were the most affected (70%), and confirmed the winter seasonality of rotaviruses, with a peak in November-January and a 'shift' of predominant rotavirus electropherotype from one rotavirus season to another. In a previous study of the molecular epidemiology of rotaviruses in Bulgaria, which included a limited number of rotavirus-positive stool samples (n = 71), a change of the predominant rotavirus genotype during two consecutive rotavirus seasons was observed [START_REF] Tcheremenskaia | Molecular epidemiology of rotavirus in Central and SouthEastern Europe[END_REF].

The aim of the present study was to study the incidence and the epidemiological features of rotavirus gastroenteritis among children less than 5 years of age in the country over an extended period of time (2005)(2006)(2007)(2008) and to describe the rotavirus strain diversity.

Materials and methods

Patients and stool samples

A total of 3,130 patients hospitalized because of acute gastroenteritis in nine hospitals around the country were enrolled in the 4-year survey (January 2005 to August 2008). Stool samples were collected during the first few hours following hospitalization in sterile feces containers and were stored at +4°C before sending to the laboratory. The investigation presented herein covers three full rotavirus seasons (2005/2006, 2006/2007, and 2007/2008) and an 8-month period (January-August) of the rotavirus season 2004/2005. A rotavirus season was defined as a 12-month period starting from September of a calendar year until the end of August in the following year. Of a total of 3,082 samples, those from whom age was known were grouped into two categories, those from children under 5 years of age (90.5%, n = 2,789; age range 25 days to 60 months) and those form school-aged children, youths and adults (9.5%, n = 293; age range 61 months to 70 years). The age of the remaining 48 patients was unknown. An acute gastroenteritis episode was defined as the appearance of three or more loose/watery stools or forceful vomiting over a 24-h period. A 10% (wt/vol or vol/vol) stool suspension was prepared in RNase-free water. Stool suspensions were vortexed and then clarified by centrifugation for 10 min at 4,000g and the supernatant was stored at +4°C until further processing.

Detection of rotavirus group A antigen

The commercial enzyme immunoassay RIDASCREEN ® Rotavirus (R-Biopharm, Darmstadt, Germany) for the qualitative determination of rotaviruses group A was used for stool specimens screening following the manufacturer's instructions.

Rotavirus G and P determination

RNA extraction

After vortexing and centrifugation at 4,000g, the supernatant was used for viral RNA extraction by TRIzol ® LS reagent (Invitrogen, USA) or the QIAamp Viral RNA Mini Kit (QIAGEN, Hilden, Germany), as described by manu-facturers. Extracted viral RNA was resuspended in 50 μl RNase-free water and was stored frozen at -70°C.

G and P genotyping

The rotavirus G and P genotypes were identified by using two different protocols. Positive and negative controls were also included. One-step semi-nested multiplex reverse transcriptase polymerase chain reaction (RT-PCR) using SuperScript One-Step RT-PCR with Platinum ® Taq (Invitrogen, USA) and one set of primers described by Das et [START_REF] Cunliffe | Rotavirus G and P types in children with acute diarrhea in Blantyre, Malawi, from 1997 to 1998: predominance of novel P[6]G8 strains[END_REF][START_REF] Das | Characterization of rotavirus strains from newborns in New Delhi, India[END_REF][START_REF] Gentsch | Identification of group A rotavirus gene 4 types by polymerase chain reaction[END_REF] was used for the period 2005-2007. Since 2007, a random primed reverse transcription step followed by a two-step PCR procedure was implemented where the consensus primers VP7-F/VP7-R and VP4-F/VP4-R were used for the first PCR reactions for G and P typing, respectively, while the second PCR assays were performed by using mixes of one consensus and type-specific primers as follows: for G genotype determination (VP7-R and aBT1/ G1, aCT2/G2, mG3/G3, aDT4/G4, aAT8/G8, mG9/G9, G10, G12) [START_REF] Iturriza-Gómara | Reassortment in vivo: driving force for diversity of human rotavirus strains isolated in the United Kingdom between 1995 and 1999[END_REF][START_REF] Gouvea | Polymerase chain reaction amplification and typing of rotavirus nucleic acid from stool specimens[END_REF][START_REF] Iturriza-Gómara | Rotavirus genotyping: keeping up with an evolving population of human rotaviruses[END_REF], and for P genotype detection (VP4-F and [START_REF] Santos | Global distribution of rotavirus serotypes/genotypes and its implication for the development and implementation of an effective rotavirus vaccine[END_REF]) [START_REF] Gentsch | Identification of group A rotavirus gene 4 types by polymerase chain reaction[END_REF][START_REF] Iturriza-Gómara | Diversity within the VP4 gene of rotavirus P[8] strains: implications for reverse transcription-PCR genotyping[END_REF][START_REF] Simmonds | New oligonucleotide primers for P-typing of rotavirus strains: strategies for typing previously untypeable strains[END_REF].

1 T-1D, 2 T-1, 3 T-1, 4 T-1, 5 T-1, P

Agarose gel electrophoresis

Detection of PCR products were achieved by staining with ethidium bromide following electrophoretic separation on 2% agarose gel in 1× TBE buffer for 1.5 h at 100 V and visualization under UV illumination. As a control for the PCR product lengths, a 50-bp molecular size marker (Invitrogen, USA) was used.

Results

Incidence of rotavirus gastroenteritis cases

During the four rotavirus seasons included in the study, a total of 3,130 fecal specimens were tested for the presence of rotavirus antigen: 164

Distribution of rotavirus G and P genotypes

Of 1,015 rotavirus-positive samples, 912 were further characterized through G and P typing. G and P types were determined from 92.4% (843/912) of the samples, whilst the remaining 7.6% (69/912) of specimens were partly (59 samples) or fully (10 samples) untypeable. A total of 781 of 912 (85.6%) specimens had a single strain and 62 (6.8%) samples contained mixed strains. A total of ten different G combinations and two P-type combinations were identified among the mixed infections, representing the commonly circulating genotypes. The G and P genotype distributions among the 912 rotavirus-positive patients is shown in Table 2.

Among single rotavirus infections, the majority of the cases (97.7%; 763/781) were caused by four of the most common rotavirus strains which have a global distribution: G1P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF], G2P [START_REF] Estes | Rotaviruses[END_REF], G4P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF], and G9P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF]. The most frequently detected strains were G9P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF], followed by G2P [START_REF] Estes | Rotaviruses[END_REF], G1P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF], and G4P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF]; however, differences in the distribution and relative incidence of these genotypes were seen year-onyear (Table 2).

Other strains comprised common G and P types but in unusual combinations, which are a result of reassortment events between common human rotavirus strains, were found in 1.8% of samples (14/781). Unusual rotavirus strains, including G8P [START_REF] Estes | Rotaviruses[END_REF], G9P [START_REF] Martella | Identification of a novel VP4 genotype carried by a serotype G5 porcine rotavirus strain[END_REF], and G12P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF], were detected in 0.5% of samples (4/781).

Age and monthly distribution of rotavirus diarrhea cases

The age distributional analysis revealed that the highest incidence of rotavirus infection was in infants between the ages of 7 and 24 months (Table 3). A significant incidence of rotavirus diarrhea was also seen in children 2-3 years old (29.8%; 126/423) and in infants under 6 months of age (27.1%; 135/498).

The seasonal distribution of rotavirus cases showed yearround circulation of rotaviruses in Bulgaria, which was unusual for countries with temperate climate (see Fig. 1). The incidence of rotavirus-positive episodes by month showed that the lowest incidence was in the spring-early summer months (May-July; 10-19%) and the highest was in the winter (January-March; 50-68%). A smaller summer peak was also observed in August: 23.5% (12/51) and 38.4% (76/198) in 2006 and 2007, respectively. Furthermore, in 2007, the summer peak was of greater magnitude than the winter peak.

Discussion

This study on the molecular epidemiology of rotaviruses circulating in Bulgaria is the first long-term systematic study of rotavirus strain diversity. The objectives of the study were to evaluate the incidence of rotavirus infections in pediatric gastroenteritis and to provide information on the rotavirus strain diversity circulating in Bulgaria over a 4-year period.

The burden of rotavirus disease

Rotavirus was detected in 32.4% of the samples tested, an incidence which is comparable with the detection rates reported worldwide (29-45%) and in Europe (27-51%) [START_REF] Parashar | Rotavirus and severe childhood diarrhea[END_REF][START_REF] Soriano-Gabarró | Burden of rotavirus disease in European Union countries[END_REF].

In Bulgaria, all hospitalizations due to acute gastroenteritis are recorded. According to the hospital discharge data of the Bulgarian Ministry of Health, for the period investigated, a total of 78,426 acute diarrhea cases were registered [36]. Of them, a total of 9,077 cases (incidence rate between 9.5 and 12.7%, depending on the year investigated) were of bacterial origin, while less than 0.2% (n = 122 for the entire period) were with suspected parasitic etiology. The remainder, 87.1-90.3%, accounted for approximately 69,000 diarrheal cases left undiagnosed. Of these, 60-80% of the cases were children under the age of 9 years, which means that, annually, 10,000- age. We established that the incidence rate of rotavirus gastroenteritis varied between 27.5 and 43.1% (average 32.4%), depending on the rotavirus season. From these data, it can be extrapolated that rotavirus infections account for 3,000-4,000 hospitalizations and more than 100,000 mild episodes of gastroenteritis each year in Bulgaria. Nonetheless, a significant diagnostic gap still remains. Investigation of other virological agents of gastroenteritis not currently investigated for in Bulgaria is likely to narrow this gap considerably and allow for a more accurate estimation of the true burden of disease attributable to the different intestinal pathogens. Although a reduction in the episodes of gastroenteritis attributable to rotaviruses was seen between 2005 and 2008, it is not clear whether the recent introduction of rotavirus testing and sampling for the current study may be a confounding factor, and only continued surveillance will allow the monitoring of this trend in the long term.

Rotavirus strain diversity

An interesting observation was the shifts of predominant rotavirus G-P genotypes from season to season. The genotype data analysis has shown that G4P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF] rotaviruses were responsible for more than 65% (25/38) of all hospitalizations during the 2004/2005 rotavirus season, but in the following years, their incidence significantly diminished and remained low until the end of the study period. In the same season, G9P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF] emerged, displacing G4 strains to become the dominant genotype (84.5%). During the following rotavirus seasons, G9 strains remained but with decreasing incidence. G2P [START_REF] Estes | Rotaviruses[END_REF] strains emerged in the middle of the 2006/2007 season, reaching an incidence of 50%. The incidence of infections with G1P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF] strains was unusually low in Bulgaria at the beginning of the study period, but increased gradually year-on-year and, by 2007/2008, became the predominant genotype. Similar seasonal fluctuations in the distribution of the predominant rotavirus strains are common events and are reported all over the world-in the United Kingdom, Hungary, Slovenia, Ireland, Belgium, Paraguay, Bangladesh, and other countries from southeastern Europe [START_REF] Bányai | Eight-year survey of human rotavirus strains demonstrates circulation of unusual G and P types in Hungary[END_REF]18,[START_REF] Iturriza-Gómara | Reassortment in vivo: driving force for diversity of human rotavirus strains isolated in the United Kingdom between 1995 and 1999[END_REF][START_REF] Parra | Diversity of group A rotavirus strains circulating in Paraguay from 2002 to 2005: detection of an atypical G1 in South America[END_REF][START_REF] Van Damme | Distribution of rotavirus genotypes in Europe, 2004-2005: the REVEAL study[END_REF][START_REF] Lennon | Changing profile of rotavirus in Ireland: predominance of P[8] and emergence of P[6] and P[9] in mixed infections[END_REF][START_REF] Rahman | Predominance of rotavirus G9 genotype in children hospitalized for rotavirus gastroenteritis in Belgium during 1999-2003[END_REF][START_REF] Rahman | Prevalence of G2P[4] and G12P[6] rotavirus, Bangladesh[END_REF]. These fluctuations may be driven by herd immunity and the emergence of a particular genotype may reflect the buildup of a sufficiently large population of susceptible children.

Epidemiology of rotavirus infection

The analysis of the epidemiological features of rotavirus gastroenteritis revealed two important observations. First, in the present study conducted in Bulgaria, rotaviruses affected predominantly children under 5 years of age (95.5%), with a peak incidence between the ages of 7 and 36 months (76.3%; 771/1011), in concordance with other countries, where the highest incidence of rotavirus diarrhea was reported in the age group 1-60 months. However, cases of rotavirus diarrhea were identified among older children and adults with an incidence of 15.4%. This is similar to recent reports from the United Kingdom, where rotavirus infections, both symptomatic and asymptomatic, were detected across the age groups [START_REF] Amar | Detection by PCR of eight groups of enteric pathogens in 4,627 faecal samples: re-examination of the English case-control Infectious Intestinal Disease Study (1993-1996)[END_REF]. Secondly, we observed year-round circulation of rotaviruses in the country, with increased incidence in the winter season and unexpected peaks preceding the rotavirus seasons. An investigation of the presence of anti-rotavirus group A IgM antibodies in adults conducted in the United Kingdom [START_REF] Cox | Serological survey of anti-group A rotavirus IgM in UK adults[END_REF] has shown a high titer of antibodies all year round, despite their seasonal distribution in the child population. Thus, this might suggest the year-round circulation of rotaviruses and an asymptomatic or milder course of rotavirus disease during the summer months, which may be a consequence of herd immunity acquired during the winter months. Furthermore, atypical late-summer peaks of rotavi-rus gastroenteritis were observed in the investigation, which preceded the newly coming rotavirus season. For instance, in August 2006, the incidence of rotavirus hospitalizations slightly increased up to 23.5%, while in the same month in 2007, it reached 38.4%. This may suggest the introduction of a new rotavirus strain/variant into the immune-naïve population, as previously reported during the emergence of G9 strains in the United Kingdom [START_REF] Iturriza-Gómara | Reassortment in vivo: driving force for diversity of human rotavirus strains isolated in the United Kingdom between 1995 and 1999[END_REF].

In summary, based on the genotype characterization of VP7 and VP4 genes, we established that the infections caused by a single rotavirus strain accounted for 92% of cases. Four of the five common rotavirus G-P combinations, G1P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF], G2P [START_REF] Estes | Rotaviruses[END_REF], G4P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF], and G9P [START_REF] Rahman | Characterization of a novel P[25],G11 human group A rotavirus[END_REF], constituted 97.7% of all rotavirus strains in circulation. Therefore, each of the two licensed rotavirus vaccines should be effective in the prevention of rotavirus disease in Bulgaria. The monovalent Rotarix ® (GlaxoSmithKline Biologicals) and the pentavalent human-bovine reassortant vaccine Rotateq ® (Merck and Sanofi Pasteur MSD) will provide adequate (heterotypic or homotypic) protection toward the five globally distributed G-P combinations. Moreover, the introduction of rotavirus vaccine into the Bulgarian National Immunization Program, either as mass vaccination or as a recommended vaccine, will significantly decrease rotavirus diarrhea cases among children or associated medical and non-medical expenses. In conclusion, enhanced and long-term investigation of the incidence of rotavirus infections, the changes in rotavirus strain diversity, and of the epidemiological characteristics (variation in regional and temporal distribution) of rotavirus infections in Bulgaria are needed in order to be able to monitor the success of any future vaccination program in Bulgaria.

Fig. 1

 1 Fig. 1 Monthly distribution of acute gastroenteritis cases tested (gray bars) and cases with confirmed rotavirus etiology (black bars). The ovals delineate the atypical summer peaks of rotavirus gastroenteritis

  al. and Cunliffe et al. (for G genotypes) and Gentsch et al. (for P genotypes)

Table 1

 1 [START_REF] Chouikha | Group A rotavirus strains circulating in the eastern center of Tunisia during a tenyear period (1995-2004)[END_REF],000 acute gastroenteritis cases of unknown origin are reported in Bulgaria among hospitalized children less than 9 years of Distribution of stool specimens tested and rotavirus-positive samples during the four rotavirus seasons included in the study

	Table 2 Distribution of rotavi-
	rus strains according to G-P	Genotype	No. of rotavirus strains			Total
	genotype combination and sea-
	son of circulation		Rotavirus season
			2004/2005	2005/2006	2006/2007	2007/2008	No.	%
		Single infections, total	38	161	288	294	781	85.6
		Common human rotaviruses					763	97.7
		G1P[8]	7	16	35	139	197
		G2P[4]	3	1	142	88	234
		G4P[8]	25	7	1	7	4 0
		G9P[8]	2	136	109	45	292
		Uncommon human reassortants				14	1.8
		G1P[4]	0	0	0	5	5
		G2P[8]	1	0	0	5	6
		G4P[4]	0	0	0	1	1
		G9P[4]	0	0	0	2	2
		Unusual rotavirus strains					4	0.5
		G8P[4]	0	1	0	0	1
		G9P[9]	0	0	1	1	2
		G12P[8]	0	0	0	1	1
		Mixed infections, total	3	6	17	36	62	6.8
		G type	1	6	5	20	32
		Р type	2	0	7	11	20
		dual G + dual Р type	0	0	5	5	1 0
		Untyped strains, total	3	8	36	22	69	7.6
		G untypeable	0	5	7	10	22
		Р untypeable	3	2	27	5	37
		G + Р untypeable	0	1	2	7	1 0
		Total	44	175	341	352	912	100

  The ovals delineate the atypical summer peaks of rotavirus gastroenteritis

	Age (months)	Patients with acute gastroenteritis tested	Rotavirus-positive patients
			No.	%
	0-3	194	48	24.7
	4-6	304	87	28.6
	7-12	704	294	41.8
	13-24	858	351	40.9
	25-36	423	126	29.8
	37-48	185	31	16.8
	49-60	121	29	23.9
	> 61	293	45	15.4
	Without information	48	4
	Total	3,130	1,015

Table 3

 3 Age distribution of the patients with acute gastroenteritis tested and patients with confirmed rotavirus infection during the four rotavirus seasons

Eur J Clin Microbiol Infect Dis (2010) 29:555-562

200 Jan uar y Fe bru ary Ma rch Ap ril Ma y Jun e Jul y Au gus t Se pte mb er Oc tob er No vem ber De cem ber Jan uar y Feb rua ry Ma rch Ap ril Ma y Jun e Jul y Au gus t Se pte mb er Oc tob er No vem ber De cem ber Jan uar y Feb rua ry Ma rch Ap ril Ma y Jun e Jul y Au gus t De cem ber Jan uar y Feb rua ry Ma rch Ap ril Ma y Jun e Jul y Au gus t Se pte mb er Oc tob er No vem ber Acknowledgment This work was partially supported by the European Rotavirus Surveillance Network. The study was also supported with grants from the World Health Organization (WHO) and the European Society of Clinical Microbiology and Infectious Diseases (ESCMID) for the provision of fellowships to ZM in the MRC/MEDUNSA Diarrhoeal Pathogens Research Unit, Medical University of Southern Africa, Pretoria, Republic of South Africa, and the Enteric Virus Unit, Virus Reference Department, Centre for Infections, Health Protection Agency, London, United Kingdom.

The authors are thankful to D. Steele (PATH, USA), J. Gray (HPA, United Kingdom), L. Fiore (Instituto Superiore di Sanita, Italy), F.M. Ruggeri (Instituto Superiore di Sanita, Italy), K. Bányai (Veterinary Medical Research Institute, Hungary), A. Sanchez Fauquier (Instituto de Salud Carlos III, Spain), and Mathew Esona (Division of Viral Diseases, NCIRD, CDC, USA) for their encouragement and help during the entire work.