Aurélien Gautreau

Alain Barrat

Marc Barthélemy

Global disease spread: statistics and estimation of arrival times

Keywords: Epidemiology, Complex networks, World Airport Network

We study metapopulation models for the spread of epidemics in which different subpopulations (cities) are connected by fluxes of individuals (travelers). This framework allows to describe the spread of a disease on a large scale and we focus here on the computation of the arrival time of a disease as a function of the properties of the seed of the epidemics and of the characteristics of the network connecting the various subpopulations. Using analytical and numerical arguments, we introduce an easily computable quantity which approximates this average arrival time. We show on the example of a disease spread on the world-wide airport network that this quantity predicts with a good accuracy the order of arrival of the disease in the various subpopulations in each realization of epidemic scenario, and not only for an average over realizations. Finally, this quantity might be useful in the identification of the dominant paths of the disease spread.

Introduction

In our modern world, the existence of various transportation means has strongly affected the way in which infectious diseases spread among humans. In fact, it has become unavoidable to take into account in the study of the geographical spread of epidemics the various long-range heterogeneous connections typical of modern transportation networks, This naturally gives rise to a very complicated evolution of epidemics characterized by heterogeneous outbreaks patterns [START_REF] Cohen | Changing patterns of infectious disease[END_REF][START_REF] Cliff | Time, travel and infection[END_REF]Colizza et al., 2006a,b), as recently documented at the world-wide level in the SARS case (http://www.who.int/csr/sars/en).

Such complex phenomenon can be tackled at different granularity levels, including very detailed agent-based simulations, complicated social and spatial structures, complex contact networks, etc. [START_REF] Anderson | Infectious diseases in humans[END_REF][START_REF] Hethcote | Gonorrhea: transmission and control[END_REF][START_REF] Kretzschmar | Measures of concurrency in networks and the spread of infectious disease[END_REF][START_REF] Keeling | The Effects of Local Spatial Structure on Epidemiological Invasions[END_REF]Pastor-Satorras & Vespignani, 2001a,b;[START_REF] May | Infection dynamics on scale-free networks[END_REF][START_REF] Ferguson | Planning for smallpox outbreaks[END_REF][START_REF] Meyers | Network theory and SARS: predicting outbreak diversity[END_REF][START_REF] Chowell | Scaling laws for the movement of people between locations in a large city[END_REF][START_REF] Eubank | Modelling disease outbreaks in realistic urban social networks[END_REF]. In particular, a very important class of models in modern epidemiology describes the propagation between different interconnected subpopulations with the help of the so-called metapopulation models. In these models, the disease spreads inside each subpopulation (often assumed homogeneously mixed) and is transmitted between different subpopulations by a "coupling" depending on the model. In the case of human infectious diseases, this coupling is caused by the movements of individuals and depends on the transportation means relevant at the spatial scale chosen for the description of the spread.

In the context of large scale spread, air-transportation represents a major channel of epidemic propagation, as pointed out in the modeling approach to global epidemic diffusion of [START_REF] Rvachev | A mathematical model for the global spread of influenza[END_REF]; [START_REF] Baroyan | An attempt at large-scale influenza epidemic modelling by means of a computer[END_REF] and by similar studies on the behavior of specific outbreaks such as pandemic influenza, HIV or SARS [START_REF] Longini | A mathematical model for predicting the geographic spread of new infectious agents[END_REF][START_REF] Grais | Assessing the impact of airline travel on the geographic spread of pandemic influenza[END_REF][START_REF] Brownstein | Empirical Evidence for the Effect of Airline Travel on Inter-Regional Influenza Spread in the United States[END_REF][START_REF] Flahault | A method for assessing the global spread of HIV-1 infection based on air-travel Math[END_REF][START_REF] Hufnagel | Forecast and control of epidemics in a globalized world[END_REF]Colizza et al., 2006aColizza et al., ,b, 2007a)). In this case, the relevant metapopulation model considers as subpopulations the inhabitants of the various world cities and the air travel of infectious individuals results in the spread of the disease from one city to another.

Taking into account the global airline transportation infrastructure in metapopulation models has recently become possible thanks to the availability and analysis of large scale databases [START_REF] Barrat | The architecture of complex weighted networks[END_REF][START_REF] Guimerà | The worldwide air transportation network: Anomalous centrality, community structure, and cities' global roles[END_REF], and to the always increasing computer capacities. The word wide airport network (WAN) is described by a complex weighted graph, in which the airports are the vertices and the weighted links represent the presence of direct flight connections among them (the weights corresponding to the number of available seats on each connection). This network, which consists of more than 3000 nodes and 18000 links, displays small-world properties as well as strongly heterogeneous topology and traffic properties. For many different dynamical phenomena occurring on complex networks [START_REF] Albert | Statistical mechanics of complex networks[END_REF][START_REF] Dorogovtsev | Evolution of networks: From biological nets to the Internet and WWW[END_REF]Pastor-Satorras & Vespignani, 2004;[START_REF] Boccaletti | Complex networks: structure and dynamics[END_REF], the presence of such emerging properties has been shown to imply the breakdown of the standard results. This is particularly true for epidemiology (Pastor-Satorras & Vespignani, 2001a,b;[START_REF] May | Infection dynamics on scale-free networks[END_REF] where classical results about epidemic thresholds do not apply if the heterogeneity of the network is too large. This therefore calls for a systematic investigation of the epidemic spread in the framework of metapopulation models defined on complex networks such as the world airport network. Such investigations can be carried out in various parallel and complementary directions. On the one hand, the availability of large scale datasets has recently allowed for the development of tools for intensive computational epidemiology (Colizza et al., 2006a(Colizza et al., ,b, 2007a)), which can be used for example for scenario evaluations (Colizza et al., 2007a) and risk assessment. On the other hand, these computational tools can also be used for the study of fundamental properties of these metapopulation models, and to understand the main mechanisms of the disease spread and the role of the various complex properties of the transport network. In particular, Colizza et al. (2006a,b) have investigated how the properties of the global air travel network affect the heterogeneity and the predictability of spreading patterns. Very recently moreover, the effect on the epidemic threshold of both the network properties and the endogenous reaction was studied using analytical and numerical tools typical of statistical physics [START_REF] Colizza | Reaction-diffusion processes and meta-population models in heterogeneous networks[END_REF]Colizza & Vespignani, 2007).

In this paper, we tackle the problem of determining the arrival time of an epidemics spread in a metapopulation model consisting in subpopulations (cities) coupled by a transportation network. We stress that this study is different from the ones of [START_REF] Crépey | Epidemic variability in complex networks[END_REF] and [START_REF] Barthélemy | Dynamical patterns of epidemic outbreaks in complex heterogeneous networks[END_REF], in which each node of the network is an individual (and not a full subpopulation), and where links represent the possibility of the contamination between individuals. The main result of the present paper is to exhibit an easily computable quantity depending only on the networks characteristics and which approximates the average arrival time. This quantity is only an approximation of the exact arrival time but we show on the example of the world-wide disease spread that it allows to rank the cities according to the arrival of the disease.

The paper is organized as follows. We present in section 2 the Rvachev-Longini metapopulation model for the epidemics spreading, and in section 3, we study the arrival time of a disease for simple topologies such as a one-dimensional line of cities. We then extend these results to the case of complex networks and we propose a quantity which is a good approximation to the average arrival time in the various nodes or subpopulations connected through an arbitrary network of connections. We finally show in section 4, using numerical simulations of epidemic spreading on the worldwide airport network, that our approximation is not only valid for the average arrival time, but is also able to give the order of arrival of the disease in the various cities with a good precision for each spreading realization. A brief account of some of the results presented here can be found in [START_REF] Gautreau | Arrival Time Statistics in Global Disease Spread[END_REF].

The Rvachev-Longini metapopulation model

The Rvachev-Longini model [START_REF] Rvachev | A mathematical model for the global spread of influenza[END_REF] was initially introduced to describe the spread of the 1968-69 Hong-Kong flu at the worldwide level. This metapopulation model uses two different levels of description: the various subpopulations are given by the inhabitants of the 58 large cities corresponding to the 58 largest airports, and homogeneous mixing is assumed at the individual city level. In this model, the propagation of the disease from one subpopulation to another is due to individuals traveling on the airtransportation network between these airports [START_REF] Rvachev | A mathematical model for the global spread of influenza[END_REF].

The evolution of the number I i of infectious individuals in each city i can thus be written as the sum of two terms

∂ t I i = K({X i }) + Ω({I j }) , (1)
where the first term K of the rhs describes the (epidemic) reaction process inside each subpopulation (city), due to the interaction of individuals in the various possible states X (X = S, L, I, R... depending on the population compartimentalization into Susceptible, Latent, Infected, Recovered... individuals).

The second term of the rhs represents the incoming and outgoing fluxes of infectious individuals to and from other cities j. This model therefore considers a simplified mechanistic approach with a markovian assumption in which individuals are not labeled according to their original subpopulation, and where at each time step the same traveling probability applies to all individuals in the subpopulation, without any memory of their previous locations [START_REF] Rvachev | A mathematical model for the global spread of influenza[END_REF][START_REF] Hufnagel | Forecast and control of epidemics in a globalized world[END_REF]Colizza et al., 2006a,b). The travel term Ω depends therefore on the air-transportation network: if the weight w ij represents the number of passengers traveling from i to j per unit of time (w ij = 0 if there is no direct connection between i and j), and N i the population of the city i, it is reasonable to assume that the probability per unit time that an individual in city i travels to city j is given by w ij /N i . In the case of a simple SI model, the time evolution of the number I i of infectious infividuals in city i is therefore given by:

∂ t I i = λI i (t) N i -I i (t) N i + j w ji N j I j - j w ij N i I i , ( 2
)
where λ is the transmission rate. Similar equations can be written for the other compartments in the various possible models (SIS, SIR...), by modifying the first term K of the rhs in Eq. ( 2) accordingly.

This original formulation has an important drawback. Indeed, it is deterministic, since it considers only expectation values, while both epidemic spreading and travel of individuals are inherently stochastic processes. The number of infectious individuals is then treated as a continuous variable and, even if the initial condition (at t = 0) of the spreading consists in one single infectious individual in a given city i 0 , all cities have a non-zero density of infectious at any time t > 0. In order to tackle the problem of the arrival time of the spread in the various subpopulations, this pathology has to be accounted for. A first possibility, already presented by [START_REF] Rvachev | A mathematical model for the global spread of influenza[END_REF], to avoid this unrealistic situation consists in considering that each travel term w ij I i /N i is present only when w ij I i /N i > 1, i.e. in putting a threshold term θ(w ij I i /N i -1), corresponding to the fact that the expected number of travelers has to be larger than unity. In this approach however the "quantity" of individuals traveling is still a continuous variable. Another possibility, which consists in treating all quantities as integers and going back to the microscopic stochastic processes, can be necessary for detailed modeling purposes (Colizza et al., 2007a). This solution is however computationally quite demanding, so that we will consider in all numerical simulations performed in this paper an intermediate framework: the travel term is treated as in the stochastic generalization described by Colizza et al. (2006a,b) (see also [START_REF] Hufnagel | Forecast and control of epidemics in a globalized world[END_REF], where the number of individuals traveling on each connection is an integer variable randomly extracted at each time step. More precisely, each of the [I i ] (the integer part of I i ) infectious individuals has a probability p ij = w ij ∆t/N i to go from i to j in the time interval ∆t. The numbers ξ ij of infectious individuals traveling on the various connections existing from i to other cities j are given by a set of stochastic variables which follow then a multinomial distribution (Colizza et al., 2006b)

P ({ξ ij } | {I i }) = [I i ]! [I i ] -j ξ ij ! j ξ ij ! × j p ξ ij ij   1 - j p ij   [I i ]-j ξ ij (3)
where the mean and variance of the stochastic variables are ξ ij ([

I i ]) = w ij ∆t[I i ]/N i and V ar (ξ ij ) = w ij N i ∆t 1 - w ij N i ∆t [I i ].
On the other hand, the endogenous growth (inside each city) will be treated for simplicity as a deterministic evolution of densities of infectious and susceptible individuals. We have moreover checked numerically that the inclusion of stochastic effects described by noise terms in the evolution equations (Colizza et al., 2006a,b) do not change our results. Finally, we will present the results essentially for the simple SI model, since arrival times will basically depend on the first stage of the disease development. We have however carried out numerical simulations as well in the SIS and SIR case, with consistent results (see section 4).

It is interesting to note that in realistic cases such as a disease propagation on the WAN, most weights are symmetric (w ij = w ji ) [START_REF] Barrat | The architecture of complex weighted networks[END_REF] but the probabilities of travel from one city to another are not (p ij = p ji ). The travel therefore effectively occurs as a random diffusion on a directed weighted complex network. Moreover, an important difference with the case of random walks on complex networks (Noh & Rieger, 2003) comes from the endogenous evolution inside each subpopulation. Indeed, as soon as an infectious individual reaches a city, a contamination process starts and the number of random infectious walkers is not constant. Such an intricate behavior has also important consequences on the existence of epidemic thresholds, as studied recently by [START_REF] Colizza | Reaction-diffusion processes and meta-population models in heterogeneous networks[END_REF].

As a last general remark, we note that the simple topological distance between nodes should naturally play a role in the context of arrival times, but does not contain all the information needed to characterize such a process, nor does a priori the optimal weighted distance, which takes into account the weights [START_REF] Wu | Optimal paths in complex networks with correlated weights: The worldwide airport network[END_REF] but not the populations nor the endogenous epidemic evolution. Moreover, since most transportation networks are small-world networks (even when taking into account weights), many cities lie at the same topological distance from a given seed, but will be reached at very different times.

Arrival times statistics

In this section, we study the arrival time of the first infectious individual in the various cities, starting with very simple topologies for the transportation network.

Two cities : arrival time distribution

We start with the case of two connected cities (0 and 1). We denote the link weight by w 01 = w, the population of the first city by N 0 = N, the number of infectious individuals in city 0 by I 0 = I. The probability to jump from 0 to 1 during the time interval ∆t is p = w N ∆t. The initial condition is given by I 0 (t = 0) = I 0 , I 1 (t = 0) = 0, i.e., infectious individuals are present only in city 0. We consider that the travel events occur as instantaneous jumps (of probability p for each individual) at discretized times. The probability that the first infectious individual arrives at time t 1 = t = n∆t in city 1 is then given by

P (t 1 = n∆t) = 1 -(1 -p) [I 0 (n∆t)] × n-1 i=1 (1 -p) [I 0 (i∆t)] , (4)
which expresses the fact that at least one successful "jump" from 0 to 1 of an infectious individual occurs at time n∆t, and none at previous times. In Eq. ( 4), the quantity [I 0 ] denotes the integer part of I 0 . In order to obtain the probability density of the arrival time in city 1, we first note that in realworld systems the number of travelers is usually small with respect to the total population of a city: p = w∆t/N ≪ 1 (with any reasonable ∆t). In this limit, Eq. ( 4) becomes

P (t 1 = t) = p[I 0 (t)]e -p 0<i<n [I 0 (i∆t)] .
(5)

We will also assume that the travel probability is sufficiently small so that the number of infectious individuals I 0 (t) in the city 0 grows substantially before the city 1 is contaminated. Writing 1 ≪ I 0 (t) ≪ N (the second inequality is due to the fact that in realistic cases only small fractions of the population are contaminated), we thus consider the following approximation: [I](t) ≈ I 0 (t) ≈ I 0 e λt (as long as t < t 1 ). Using as well the standard approximation ∆t 0<i<n [I 0 (i∆t)] = t 0 I 0 (τ ) dτ , we obtain for I 0 = 1

P (t)dt = w N e λt-w Nλ e λt Θ(t)dt , (6)
where Θ(t) is the Heavyside function which ensures the positivity of the arrival time. The distribution ( 6) is a Gumbel distribution with average

t 1 = 1 λ ln N λ w -γ
, where γ is the Euler constant. The contribution of the negative t in the distribution has to be negligible before 1, which reads 0 -∞ P (t) dt = w N λ ≪ 1. This range of validity corresponds also to the hypothesis of λ t 1 ≫ 1 (the time to increase I 0 by one individual is small with respect to the time scale of the epidemic arrival in city 1).

Figure 1 displays the results of numerical simulations of the propagation between two cities, using discrete stochastic travel events as described in section 2. The whole distribution of arrival times in the second city shows a very good agreement with the form (6), obtained through the continuous approximations detailed above. When w/(Nλ) is not small enough, stronger deviations are obtained; this is expected since the hypothesis 1 ≪ I 0 (t) for t < t 1 is less valid. Nonetheless, the overall shape of the distribution is still in good agreement (Inset of Fig. 1).

One-dimensional line

We now consider a one-dimensional line of subpopulations i ≥ 0, with an initial condition given by one infectious individual in city 0. We denote by w i the travel flux between city i and city i + 1, t i the arrival time of the first infectious in city i, and ∆ i = t i -t i-1 (t 0 = 0). The quantity t n is thus the sum of the random variables ∆ i for i = 1, ..., n, and the average is given by t n = i=1,...,n ∆ i . These random variables ∆ i are a priori correlated and not identically distributed which implies that the central limit theorem can not be used in general.

Homogeneous line

In the case of a homogeneous line with uniform populations and weights (N i = N, w i = w), one expects that the distributions of ∆ n becomes independent of n at large n, with a well-defined average lim n→∞ ∆ n = ∆ , so that t n ≈ n ∆ . The issue then is the computation of ∆ . Two time scales only are present: N/w which represents the average time for an individual to travel from one city to the next, and 1/λ which represents the typical transmission time of the disease. Dimensional analysis thus implies that the adimensional quantity λ ∆ has the form

λ ∆ = F w Nλ . (7)
In order to estimate the unknown function F , a first possibility is to consider that at large n the spread consists in an evolving epidemic front obeying the continuous equation

∂ t I(x, t) = λSI +w∂ 2 x I(x, t).
Looking for a traveling wave solution [START_REF] Murray | Mathematical Biology II[END_REF] leads to a front speed v = 2 λw N , and thus to

∆ f ront ≈ 1 2λ Nλ w . ( 8
)
Another possibility consists in using the results of the previous subsection, and assume that the average of ∆ n remains close to the arrival time in the first city ∆ n ≈ ∆ 1 = t 1 which yields

∆ Gumbel = 1 λ ln Nλ w -γ . ( 9
)
This approximation neglects the fact that I n-1 (t) increases between t n-1 and t n due both to the endogeneous growth in city n -1 and to the arrival of infectious individuals from city n -2, while for the computation of t 1 , only the endogeneous growth of I 0 has to be considered. It can therefore be expected that ∆ Gumbel will overestimate the real ∆ .

Finally, we can also consider a deterministic formulation, in which travel between i and i + 1 occurs only if there are enough infectious individuals in i, i.e. I i (t) has reached a certain threshold θ, and is then treated as continuous.

In this framework, each city is considered infected only if the number of infectious individuals is above θ, and t n is defined as the first time that I n reaches the threshold: I n (t n ) = θ. Between t n-1 and t n , no travel can therefore occur out of n, and we can write

∂ t I n = λI n + w N I n-1 . During [t n-1 , t n ],
under the realistic hypothesis that I n ≪ N n , and that I n-1 is in a first approximation given by I n-1 = e λ(t-t n-1 ) , we obtain the equation 1 = w N ∆ n e λ∆n . The difference between the arrival times in two successive cities is therefore given by

∆ det = 1 λ W Nλ w ( 10
)
where W is known as the Lambert W function.

In order to test these various analytical approaches, we consider numerical simulations of the stochastic model described in section 2, for identical cities of population N = 10 6 located on a one-dimensional line, with uniform travel fluxes w between successive cities. The measured average arrival times are, as expected, proportional to n at large n (not shown). Figure 2 displays the corresponding slope as λ ∆ versus ln(λN/w). Various values of N, λ and w with the same ratio λN/w yield the same λ ∆ , as predicted from the theoretical dimensional analysis (7). As shown in Fig. 2, ∆ f ront is in agreement with the simulations only at small λN/w: the slowest the travel, the less a spatial continuous approximation is valid. On the other hand, both ∆ det and ∆ Gumbel display reasonable agreement, and in particular correctly capture the increase with ln(λN/w) for large λN/w. As expected, ∆ Gumbel slightly overestimates ∆ . We also show in the inset of Fig. 2 the correlations between the ∆ n 's. These correlations vanish rapidly with the distance along the line.

Negative correlations can however be observed between ∆ n-1 and ∆ n . Such phenomenon can be understood as follows. Assume that, in a given spreading realization, ∆ n-1 is small. In this case, I n-2 (t n-1 ) will be unusually small too.

The "reservoir" of infectiousness defined by city n -2 will thus transmit less infectious individuals to n -1 at times t > t n-1 , and therefore the subsequent contamination of city n will be slower, leading to a larger ∆ n .

The numerical simulations allow also to measure the whole distributions of arrival times and of their intervals. In Fig. 3, we show P (∆ n ) for different values of n. The Gumbel shape is only valid for n = 1, and for large n more symmetric (Gaussian-like) distributions are obtained. On the other hand, Fig. 4 shows that the distribution of the arrival times themselves display asymmetric Gumbel-like shapes even at large n.

General populations and weights

We now consider the general case of heterogeneous populations N i and fluxes w i . The application of the same considerations as in the case of a homogeneous line leads to t n = n i=1 ∆ i . Since our main goal is to understand the case of complex networks which usually have the small-world property (i.e. a small diameter varying typically as log N), we will now focus on the properties of t n at small n. Assuming that the average of ∆ i remains close to the formula for the arrival time in the first city ( ∆ i ≈ t 1 (λ, w i , N i ) ), we obtain

λ t n ≈ χ n ≡ ln n-1 i=0 N i λe -γ w i . ( 11
)
This expression links the expected arrival time in city n of a stochastic spreading event to the quantity χ n which depends on the properties of the line on which the spreading propagates. It is symmetric by permutations of weights and populations, and Figure 5 shows in fact that the whole distribution of arrival times, obtained through numerical simulations of the discrete stochastic dynamics, and not only its average, respects this symmetry. In the numerical simulations, heterogeneous populations and travel fluxes are uniformly distributed (w i ∈ [10,2000] and N i ∈ [10 5 , 2.10 7 ]). Figure 5 shows that the distribution of arrival times is invariant when one replaces (i) all the random weights by their geometrical mean w = ( i=0...n w i ) 1/n ; (ii) all the random populations by their geometrical mean N = ( i=0...n N i ) 1/n ; (iii) all weights by w and all populations by N. The ratios of the average arrival times for these different sets t

({w i }, {N i }) / t(w, {N i }) , t({w i }, {N i }) / t({w i }, N)
and t({w i }, {N i }) / t( w, N) stay very close to 1, with deviations at most of the order of 5%.

We test in Fig. 6 the relation ( 11) which predicts the average arrival time measured in numerical simulations of the stochastic model. We use for simplicity uniform populations N = 10 6 , and we measure, for each fixed set of weights, the arrival time in each city, averaged over 10 3 stochastic realizations of the spread. Figure 6 displays λ t i versus χ i for the various cities of a 5 cities line, for 100 different random sets of weights (each point corresponds to one given set of weights). We first note that λ t i is a monotonous function of χ i . The figure also clearly shows that χ systematically overestimates the average arrival time. This is due, as already explained in the discussion of Eq. ( 9), to the fact that the use of ∆ Gumbel neglects correlations, and in particular the travel between cities n -2 and n -1 for the estimate of the arrival time in n. Since this travel term increases I n-1 , it increases the probability that an infectious travels from n -1 to n, and therefore reduces t n . This overestimation effect is not present for the first group of points that lie on the diagonal of the figure and correspond to cities directly connected to the seed. For this first group of cities, the average arrival time is indeed correctly given by χ. On the other hand, the small overestimation can already be observed for cities at distance 2 of the seed. For example, points A and B highlighted in Fig. 6 have the same χ B = χ A = χ. Point A corresponds to a city directly connected to the seed S, with a flux of w A passengers per unit time, such that χ A = ln(Nλ/w A ) -γ, while point B is obtained for a city which lies at topological distance two from the seed, with a city B ′ in between and such that χ B = ln N λ w B e γ + ln N λ w B ′ e γ . The two cases correspond to the same value of χ, but to different average arrival times: the arrival time in B is smaller than in A, even if χ A = χ B , because relation (11) neglects the travel of additional infectious individuals between S and B ′ to compute t B . Despite this systematic effect, Fig. 6 shows that relation (11) allows to define a quantity which depends only of the populations and passenger fluxes, and such that the average arrival time of the spreading is a monotonous function of this quantity. In the following, we will investigate how to generalize this quantity to more complex topologies.

From the one-dimensional line to complex networks

Equation ( 11) gives an estimate of the average arrival time in a city n connected to the seed of the spreading through links with certain fluxes, and corresponds to a sum of the quantity ln(Nλ/w) -γ along the links followed by the disease from the seed to city n. In a network, the sum can be computed along any of the paths linking the seed s to any other city j. Since these different paths correspond a priori to different fluxes and to different arrival times, a natural approach consists in approximating the average arrival time in j, starting from a given seed s, by the minimum of ( 11) over all possible paths:

λ t j ≈ χ (j|s) ≡ min {P s,j } (k,l)∈P s,j ln N k λ w kl -γ ( 12
)
where s is the seed, {P s,j } is the set of all possible paths connecting s to j, and the sum is performed over the links (k, l) along each path. In other terms, we have introduced a new weight on each oriented link (k, l) of the network: ln N k λ w kl -γ. The quantity χ (j|s) is then the weighted distance between the seed s and the node j on this non-symmetric network (such a quantity can easily be computed for any topology using the Dijkstra algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF]). Note that since the weights are real-valued, it is highly improbable that two different paths with the same sum of weights exist, so that equation (12) selects a unique path between s and j.

Differences between propagation on a one-dimensional line and on a network can be expected from the two following important topological differences between these structures: (i) paths go through nodes with degree larger than 2 and (ii) there are usually more than one path between two points. We first focus on the effect of the intermediate nodes with possibly large degrees through which the disease propagates. To this aim, we consider the simple topology consisting of a node 0 connected to k neighbors (as shown in the inset of Fig. 7). The infectious individuals, when they arrive at the central node, have multiple possible travel destinations. The probability for each destination to become infected is thus decreased; if the seed is in 0 for instance, the initial infectious individual may even travel to a different peripherical node before creating an endogenous epidemic growth in the hub. This effect can be quantified by measuring the average arrival time of the disease in a given peripheral city i in various situations. We first consider the case of a spreading phenomenon starting in a randomly chosen peripheral node, i 0 = i, and compare the average arrival time in i with the average arrival time in the one-dimensional situation in which only cities i 0 , 0 and i are present. The ratio of these two times is displayed in Fig. 7 ('seed: not hub') as a function of the degree k of the central node. We also show in the same figure the ratio of arrival times when the seed is 0, i.e. the hub itself ("seed: hub"), and the topology is either the one of the inset, or the simple two-cities configuration in which only 0 and i are present and connected. In all cases, the arrival time in i is increased by the presence of other possible connections from the hub. The effect is stronger for larger degree when the seed is the hub itself. The multiplicity of possible destinations therefore not only decreases the predictability of the spread (Colizza et al., 2006a), but also increases the average disease arrival time in a given city.

Another important point distinguishes a propagation in a network from a one dimensional line: namely, that multiple paths can often be found between two given nodes. In fact, it has already been shown by [START_REF] Crépey | Epidemic variability in complex networks[END_REF] in the case in which nodes represent individuals and not subpopulations, that the average arrival time is decreased by the presence of multiple paths, and that not only the shortest paths contribute to the spread between the seed and other nodes. To quantify this effect in metapopulation models, we consider the simple topology depicted in Fig. 8: a node A (population N A ) is connected to a node B by two different paths of length two. The first path connects A with weight w 1 to an intermediate city C of population N, which is in its turn connected to B with weight w 2 . The second path has weights w ′ 1 and w ′ 2 , and the intermediate city C ′ has population N ′ . We denote the probability that the disease spread reaches the city B at time t by P (t) if the two paths are present, by P ACB (t) if only the first path through C (weights w 1 , w 2 ) is present, and by P AC ′ B (t) if only the second path (through C ′ , weights w ′ 1 , w ′

2 ) is present. Since the epidemics reaches B from A through one path or the other, P (t) can be expressed as

P (t) = 1 - t 0 P ACB (τ ) dτ P AC ′ B (t) + 1 - t 0 P AC ′ B (τ ) dτ P ACB (t) (13)
The data shown in Fig. 4 suggest that we can assume that P ACB and P AC ′ B are Gumbel distributions, of averages χ ACB = ln

N A λ w 1 + ln N λ w 2 -2γ and χ AC ′ B = ln N A λ w ′ 1 + ln N ′ λ w ′

2

-2γ, respectively. We introduce the quantities

w eq = w 1 w 2 e γ N λ and w ′ eq = w ′ 1 w ′ 2 e γ
N ′ λ . The quantity w eq can be seen as the travel flux which, if A and B were directly connected, would yield the same arrival time distribution than the real two links AC and CB with weights w 1 and w 2 : χ ACB = ln N A λ weq -γ. We then have

P ACB (t) = w eq N A e λt- weq N A λ e λt (14)
P AC ′ B (t) = w ′ eq N A e λt- w ′ eq N A λ e λt (15)
and finally from Eq. ( 13) (using the fact that the weights are small with respect to the quantities Nλ, N ′ λ, N A λ)

P (t) = w eq + w ′ eq N A exp λt - w eq + w ′ eq N A λ e λt . ( 16
)
We note that P (t) is also a Gumbel distribution, and the average of the arrival time in B, t mp (mp=multiple paths), can be easily calculated:

λ t mp = χ mp ≡ ln N A λ w eq + w ′ eq -γ , (17)
so that the existence of two paths results in the law

e -χmp = e -χ ACB + e -χ AC ′ B . ( 18
)
We have checked this prediction by numerical simulations of stochastic spreading in the small network formed by nodes A, B, C and C ′ , in the simple case of 8 displays the ratio of the average arrival time in B, for a spread seeded in A, to the average arrival time in B when only one path (the one with weights w) is present. We consider only w ′ ≤ w since the situation with w ′ ≥ w is obtained by exchanging the two paths. As w ′ → 0, this ratio goes to 1, while, in the other extreme case w ′ → w, a decrease of the average arrival time close to 10% is obtained. Moreover, the prediction (18), shown as continuous line, is in excellent agreement with the numerical data.

N A = N = N ′ and w 1 = w 2 ≡ w, w ′ 1 = w ′ 2 ≡ w ′ . Figure
Interestingly, when more than two paths are present, Eq. ( 18) can easily be generalized to a sum over all possible paths (not necessarily of length 2):

λ t mp = χ mp with e -χmp = paths e -χ path , (19)
where χ path is computed on each path by Eq. ( 11). Figure 9 presents a comparison of numerical results with this analytical prediction, for various path multiplicities, showing a very good agreement. Similar results are obtained for paths of larger lengths.

We have also investigated numerically, on simple topologies but with random weights, how well the average arrival time in a city is correlated with χ mp computed as in ( 19). In particular, Fig. 10 displays the average arrival time in city B, connected to the seed A of the spreading by k paths of length l and random weights, as a function of both χ B given by Eq. ( 12) (i.e. computed along the path that minimizes χ) and χ mp which takes into account all paths. While χ mp is slightly more strongly correlated with t than χ B , and closer to its actual value, the improvement is not striking. We have also developed an algorithm to compute χ mp on a complex network, by taking into account the contributions of all paths of topological lengths d and d + 1 (where d is the length of the shortest path, in terms of topology, between the seed and the considered city). While the algorithmic complexity is increased (we used techniques based on the Brandes algorithm [START_REF] Brandes | A Faster Algorithm for Betweenness Centrality[END_REF]) with respect to the computation of ( 12), the obtained improvement was again significative but not striking.

4 The world-wide airport network

We now proceed to numerical simulations of the metapopulation model on the Worldwide airport network * , which displays various levels of complexity and heterogeneity [START_REF] Barrat | The architecture of complex weighted networks[END_REF][START_REF] Guimerà | The worldwide air transportation network: Anomalous centrality, community structure, and cities' global roles[END_REF]Colizza et al., 2006a,b). At the topological level, the degree distribution is broad and can be approximated by a power-law; the links' weights (fluxes) are also broadly distributed and span several orders of magnitude. Finally, the world city populations are also broadly distributed according to Zipf's law [START_REF] Zipf | Human behavior and the principle of least effort[END_REF]. All these heterogeneity levels have been shown to play relevant roles in the spread of epidemics at the worldwide level (Colizza et al., 2006a,b). We perform the numerical simulations of the stochastic spreading on the network containing the 2, 400 largest airports, which takes into account 98% of the total traffic (the WAN comprises 3, 100 airports, but we consider here only the links verifying ln N λ w > 0, which amounts to the removal of the smallest airports; in particular, a certain number of nodes have a larger traffic than inhabitants, which leads to a local breakdown of the relation between p ij and the local population and travel flow).

Average arrival time

Figure 11 displays the average arrival time in the various cities, λ t i , as a function of χ i as defined by ( 12). For each given seed, averages are done over 1, 000 stochastic realizations of the disease spread. The figure clearly shows that the value of χ determines the average arrival time (the two quantities are very strongly correlated) and various cities with the same χ are reached at the same time by the propagation. As in the case of the simple topologies studied above, χ in fact systematically overestimates the correct average arrival time, but can still be considered as a very good approximation.

Figure 11 also highlights the effect of hubs, analyzed in the previous section. Frankfurt is indeed the node with the largest degree in the network, and the arrival times at a given value of χ is larger than for other seeds with a smaller degree. This effect remains however quite small and is essentially limited to the first reached cities. The arrival times seem to be closer to our Ansatz. This stems from the delay introduced by the large degree of the initial seed, as discussed in the previous section, which increases all the arrival times, therefore slightly shifting the data upwards and partially compensating for the fact that χ overestimates these times. We finally note how very similar results are obtained for the SIR model. As expected, the important time scale in the SIR (as well as in the SIS) is λ-µ: this quantity controls the endogenous growth of the epidemic at small times. We have investigated a wide range of values for λ -µ, as well as for λ in the simple SI case. The arrival time is well determined by χ as long as the condition w/Nλ << 1 given in section 3 is satisfied.

Figure 11 conveys the result that the average arrival times are not exactly given by χ, but are at least determined to a large extent, given a seed, by this quantity. An immediate application of this result is given by a possible prediction of the order of arrival of the disease in different cities. In order to compare the list of cities ranked by the average arrival time and the same list ordered according to χ, we compute the indicator known as Kendall's τ . This indicator allows a quantitative analysis of the correlations between two rankings of n objects (Numerical Recipes, 2004) and is given by

τ = n c -n d n(n -1)/2 (20)
where n c is the number of pairs whose order does not change in the two different lists and n d is the number of pairs whose order is inverted. This quantity is normalized between -1 and 1: τ = 1 corresponds to identical rankings while τ = 0 is the average for two uncorrelated rankings, and τ = -1 is a perfect anticorrelation.

We present in the first column of Table 1 the values of τ for the lists of cities ordered according respectively to the average arrival time of the spread and to the values of χ, for various seeds. The second column of the table gives the values of τ when the lists are ordered by average arrival time and by χ mp . In both cases the values are extremely high (two random permutations of a list of size N = 2, 400 would give an indicator normally distributed with mean 0 and variance 4N +10 9N (N -1) = 2.10 -4 (Numerical Recipes, 2004)). In order to better emphasize this strong correlation between the ordered list, we also display in Table 1 the τ of Kendall between the city list ordered by average arrival time and by different distances from the propagation seed. As noted in the introduction, the topological distance d from the seed is not completely irrelevant, but other weighted distances which take into account the diversity of fluxes and populations could be expected to be more strongly correlated with the arrival time. In particular, a first possibility consists in defining the effective "length" of each edge as the inverse of the weight, ℓ ij = 1/w ij : the disease will spread more easily on an edge if many passengers travel across it. The corresponding weighted distance between nodes on the network is noted d 1/w . Since the ratios N/w moreover appear naturally in the travel probabilities, we also consider that a directed length N i /w ij can be defined on each edge, and the distance between each node and the seed can as well be computed (d N/w ). The results indicate significant correlations between the average arrival times and these various distances. The correlations are however much stronger with χ and χ mp , showing that these quantities can be used with a good confidence as an estimate of the average arrival time order. Very similar results are obtained for SIS or SIR models, with respectively τ SIS { t , χ} ≈ 0.875 and τ SIR { t , χ} ≈ 0.866, for a spreading process seeded in Hong-Kong.

Arrival order for a given realization

The previous results are valid for the average arrival times, and it is legitimate to ask about their relevance to the case of a single spreading event. Indeed, in the real-world there is not such a thing as averages over different realizations. The natural extension of our results therefore consists in checking if the quantity χ can predict the order in which the disease will reach the various nodes (cities) of the network in each stochastic realization of the spread. More precisely we can compute, for each pair of nodes (i, j) in the network, and for a given infection seed, the probability that the arrival times in i and j for the same realization of the spread, t i and t j , are correctly ordered by their values of χ, i.e. that they verify (t i -t j )(χ i -χ j ) > 0. We use the notation ∆χ (i, j) =| χ j -χ i | and Fig. 12 displays the probability f c that a couple of nodes with a χ-difference ∆χ (i, j) = ∆χ are reached by the disease in the order predicted by their values of χ. If ∆χ (i, j) = 0, no prediction can be made and we obtain indeed f c = 0.5. On the other hand, if ∆χ(i, j) is large (> 10), the two nodes are reached in the correctly predicted order in every realization of the spread.

Since not all node pairs have very different values of χ (and thus a large value of ∆χ), we also show as in Fig. 12 the cumulative distribution p > (∆χ) of the number of couples of nodes with a given value of ∆χ. For instance, for a spreading process seeded in Hong-Kong (degree 109), 75% of the couples have ∆χ > 2, and these pairs are correctly sorted with a probability larger than 80%, instead of only 50% on average if no information is available. Figure 12 moreover shows that the precision is higher when the seed has a small degree. On the other hand, if the seed is Frankfurt (of degree 240, data not shown), the situation is slightly worsened, as can be expected from the smaller global predictability of the spread starting from a hub (Colizza et al., 2006a), due to the many possible travel destinations available for the infectious individuals. In this case, 80% of the couples have ∆χ > 2, and these pairs are correctly sorted with a probability larger than 70%; the probability f c is larger than 90% for 47% of the pairs. Similar results are obtained when the endogenous growth of the disease is described by SIS or SIR models.

We also present in Fig. 13 the same quantities f c and p > (∆χ) for couples of nodes having the same topological distance from the seed. While this topological distance is indeed too simple a quantity to be really useful in the prediction of arrival times, it can seem reasonable that a node at distance 1 from the seed will anyway be reached before a node at distance 5. If large values of ∆χ were obtained only for nodes at very different topological distances from the seed, the prediction shown in Fig. 12 would not be particularly relevant. We see in Fig. 13 that nodes at the same topological distance from the seed can have very different values of χ and be therefore reached by the disease in an order which can be almost certainly predicted by measuring their values of χ. For example, Fig. 13 shows that more than 70% of the couples are well ranked with a probability larger than 70%, even if the two nodes are at the same topological distance (d = 2, 3, 4 or 5) from the seed. Results are presented for a seed with large degree, and slightly better results are obtained for smaller degree seeds.

Finally, another quantitative indication of the relevance of the quantity χ is obtained, similarly to the previous subsection, by comparing the list of cities ordered either by χ or by the arrival time of the spread in a given realization (and not the average arrival time as in the previous section). Kendall's τ is therefore now a quantity which fluctuates from one realization to the other (for a given seed), and the corresponding histograms are shown, for various seeds, in Fig. 14. The values obtained are systematically larger than 0.5, denoting a strong correlation between the lists ordered according to the arrival time and χ. Moreover, the figure highlights how seeds with larger degrees lead to smaller values of τ . This is in agreement with the slightly better prediction capacities displayed in Fig. 12 when the seed has small degree. Interestingly, this result can also be related to the issue of predictability studied by Colizza et al. (2006a,b): the predictability of a spreading process, as measured by the similarity between two stochastic realizations with the same initial conditions, has indeed been shown to be larger when the seed has a small degree.

Conclusion

In this paper, we have studied a metapopulation model for the spread of epidemics on a large scale in which subpopulations (cities) are linked by fluxes of passengers. Such models are particularly useful for the analysis of epidemics which propagate worldwide along the airline connections [START_REF] Rvachev | A mathematical model for the global spread of influenza[END_REF][START_REF] Hufnagel | Forecast and control of epidemics in a globalized world[END_REF]Colizza et al., 2006aColizza et al., ,b, 2007a)). They couple endogenous evolutions of epidemics inside each subpopulation with diffusion along the network of connections, which opens interesting new perspectives [START_REF] Colizza | Reaction-diffusion processes and meta-population models in heterogeneous networks[END_REF]. In this study, we have focused on the issue of arrival times of epidemics, as a function of the seed and of the network's characteristics. We have proposed a quantity easily computable on any network which depends only on the links weights and nodes populations, and which accounts for the average arrival time of the spread in each city. This quantity allows to sort the various subpopulations according to the order of arrival of each single realization of the spread with a very good accuracy. We note that this quantity is given by a certain weighted distance computed on a directed weighted graph, and in particular that this distance selects a shortest weighted path between the seed and the various nodes. This result could therefore shed some light on the existence and properties of "epidemic pathways" whose relevant role was already suggested by Colizza et al. (2006a,b) and which certainly deserve further work. In particular, the role and relative importance of the multiple paths should be investigated. These predictive tools could play an important role in the set-up of containment measures and policies. It would also be interesting to adapt the proposed quantity to more refined or sophisticated compartmental models [START_REF] Elveback | An influenza simulation model for immunization studies[END_REF][START_REF] Watts | Multiscale, resurgent epidemics in a hierarchical metapopulation model[END_REF], or to generalize it to different scales such as the urban scale, where nodes are places like home, work or malls [START_REF] Eubank | Modelling disease outbreaks in realistic urban social networks[END_REF].

Acknowledgments

It is a pleasure to thank Vittoria Colizza and Alessandro Vespignani for discussions. A.G. and M.B. also thank the School of Informatics, Indiana University where part of this work was performed. We are particularly grateful to Vittoria Colizza for sharing with us the data on the city populations and we also thank IATA (http://www.iata.org) for making their database available to us. A.G. and A.B. are partially supported by the EU under contract 001907 (DELIS). 1) in numerical simulations with w N λ = 10 -2 . Inset : distribution of t (1) with w N λ = 10 -1 . The continuous line is the analytical approximation given by Eq. ( 6). Here and in the following, unless specified, N = 10 6 and λ = 0.1. of the average arrival time in a fixed peripheral city i of the network shown in the inset and the average arrival time in the 1d case, vs the degree of the hub, k. Averages are done over 1000 realizations of the spreading, with parameter w N λ = 10 -2 . Circles: the seed is the hub; squares: the seed is one of the peripheral cities.

Seed { t , χ} { t , χ mp } { t , d} { t , d N/w } { t ,

Fig. 1 .

 1 Fig.1. Two cities model, distribution of arrival times t (1) in numerical simulations with w N λ = 10 -2 . Inset : distribution of t (1) with w N λ = 10 -1 . The continuous line is the analytical approximation given by Eq. (6). Here and in the following, unless specified, N = 10 6 and λ = 0.1.

Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .

 2345 Fig.2. Slope ∆ of t n vs n on a line of length 500. The average t n is taken over 500 realizations of the spreading, and the slope is measured at large n (only the cities with n ≥ 100 are considered). The error bars represent the variance of ∆. λ ∈ [10 -2 , 1] and w ∈ [10, 10 5 ], N = 10 6 . Full line: ∆ det . Dotted line: ∆ Gumbel . Dashed line: ∆ f ront . Inset: Correlations C ij = ∆ i ∆ j -∆ i ∆ j versus | i -j | on a line of length 500, the average is done over 500 realizations.

Fig. 6 .Fig. 7 .

 67 Fig.6. λ t vs χ for 5 cities connected on a line. Each point correspond to the average arrival time in one of the city of the line, for one realization of the weights w, averaged over 10 3 realizations of the epidemic spread. Results are presented for 100 different random sets {w i }.

Fig. 8 .Fig. 9 .Fig. 10 .

 8910 Fig.8. Ratio of the average arrival time in B (the seed is A) and the average arrival time in the 1d case (when only the path ACB is connected), versus w ′ /w (in the simulations w 1 = w 2 = w, w ′ 1 = w ′ 2 = w ′ , and N = N ′ = N A ). Averages are done over 1, 000 realizations of the spreading, with w = 10 3 and w N λ = 10 -2 . All cities have the same population N = 10 6 . Full line: theoretical prediction (17).

Fig. 11 .Fig. 12 .Fig. 13 .

 111213 Fig. 11. λ t versus χ on the WAN for diseases starting in different cities. Each red dot corresponds to a city and averages are done over 1, 000 realizations of the spread. Crosses are an average over cities with the same χ. The starting seeds are Frankfurt (degree 240), Bucharest (degree 56) and Hong-Kong (degree 109). When the initial seed is a hub, the average arrival time is larger, especially in the first reached cities, due to the effect of the degree. We also show λ t versus χ for a SIR model seeded in Hong-Kong.

Fig. 14 .

 14 Fig. 14. Histogram of Kendall's τ for the list of cities sorted according either to the arrival time in a given realization or to χ. Different seeds yield different histograms.

* We have as well considered artificial networks, with similar numerical results. Some arbitrariness in the distribution of weights and city populations is necessary for articial networks so that we prefer to present data obtained with a real-world network, in which all quantities stem from real data.


Lionel Seinturier

Nicolas Pessemier
email: pessemie@lifl.fr

Didier Donsez
email: didier.donsez@imag.fr

Towards a Reference Model for Implementing the Fractal Specifications for Java and the .NET Platform

Introduction

So far, several implementations of the Fractal specifications have been proposed. These implementations propose frameworks for programming with Fractal in a target language (Java, C, Smalltalk, C++). The general principles of implementing the specifications are common to all these frameworks. However, as far as we know, no concrete piece of software or no common set of internal interfaces have ever been shared between them.

In this paper, we report on a reference model which has been set up to support several implementations on the Fractal Specifications. This model has been derived to build a Java personality called AOKell, and FractNet, a personality for the languages of the .NET framework. The originality of this model is to be based on the concepts of Aspect-Oriented Programming (AOP) [START_REF] Kiczales | Aspect-oriented programming[END_REF].

Section 2 presents the general architecture of the reference architecture. Section 3 reports on the two personalities, AOKell and FractNet. Section 4 shows some performance measurements. Section 5 concludes this paper.

Architecture

Background

Two dimensions can generally be found in a Fractal implementation: the business dimension and the control dimension. The former is responsible for implementing the core functionalities of the application, while the latter provides a level of supervision and management on the application.

Fractal defines some artefacts to implement the business dimension with components: the notion of a provided interface, a required interface, a component type or an interface type. The control dimension, so called controllers, deals with functionalities such as starting/stopping the components, managing their bindings, their attributes, or the containment hierarchy. The scope of the control dimension is fully open: there is no restriction on what the control dimension should or shouldn't do. It is up to the framework developer to design and implement the control functions s/he needs.

The role of a Fractal framework implementor is then: (1) to support the API defined in the specifications, (2) to provide some structures for implementing controllers, (3) to provide some mechanisms for integrating the control and the business dimensions, i.e. to apply the control functions to the components. For example, Julia [START_REF] Bruneton | An open component model and its support in Java[END_REF], the Fractal reference implementation, uses mixins [START_REF] Bracha | Mixin-based inheritance[END_REF] to perform this integration. The originality of the approach presented in this paper is to use aspect-oriented programming (AOP) [START_REF] Kiczales | Aspect-oriented programming[END_REF] for integrating the control dimension with the business one.

Aspects [START_REF] Kiczales | Aspect-oriented programming[END_REF] are software entities modularizing concerns which are said to be crosscutting. The code of such concerns is not properly modularized in one unique entity (class, procedure, function, etc.), but is spread around many different locations. This phenomenon, known as code scattering, hinders the development and the maintenance of applications. AOP promotes a way of remodularizing these concerns, with the notion of an advice code, and a way of integrating them in the rest of the application, with the notion of a pointcut. By this way, AOP complements existing programming styles, object-oriented, procedural, functional, etc., by removing scattering and producing software which is more modular. AOP has been the subject of many studies in the past years and as been applied to many domains, including middleware (just to name a few works in this area, see for example [START_REF] Colyer | Large-scale AOSD for middleware[END_REF], [START_REF] Zhang | Quantifying aspects in middleware platforms[END_REF] or [START_REF] Kulesza | Reengineering of the JaWS web server design using aspect-oriented programming[END_REF]).

Aspect-Oriented Architecture

Our reference model defines three layers: the two above-mentioned dimensions, application and control, and a middle layer which is aspect-based (see Figure 1). This aspect layer is responsible for modularizing and integrating the control functions into the application. This layer supervises and controls the applications and delegates the concrete realizations of these functions to the control layer. The aspect layer defines one aspect per control functions. For example, the seven basic controllers of Fractal (attribute, binding, factory, component, content, naming, super) are each associated to an aspect. Any new controller which needs to be added must comes with its aspect. The aspects perform two basic tasks: (1) code advising, (2) code injection. Advising is an AOP related term which refers to the ability of intercepting some events of the application (such as method calls, method executions, field writes, etc.) and performing treatments before and/or after these events. Code injection refers to the ability of extending the implementation of a component with new methods, interfaces or fields.

Code advising is used in the aspect layer to control and supervise the application layer. For example, method executions can be intercepted, blocked or released depending on a particular life cycle policy. Code injection is used to enhance the functionalities of a component with new interfaces, such as an interface for managing binding, names, or any other control interface.

In our model, while aspects supervise and control, they do not realize themselves the control functions. They delegate it to the control layer. For example, the aspect in charge of the management of bindings, injects a stub for the BindingController interface which delegates the implementation of this interface to an object of the control layer. This pattern is a frequently used practise of AOP: it allows decoupling the integration logic (aspect layer) from the concrete realization of this logic (control layer).

This architecture are been declined into two implementations, one for Java (AOKell) and one for the .NET platform (FractNet). They are presented in the next sections. Furthermore, some parts of AOKell have been reused in FractNet. The performances of these implementations are summarized in Section 4.

Implementations

This section reports on the two implementations, AOKell and FractNet, which have been derived for the previous reference model. They are briefly presented in sections 3.1 and 3.2. Section 3.3 provides a short description of the code elements which have been shared between the two implementations.

AOKell

AOKell1 is the Java implementation of our reference model. AOKell performs a compile-time (CT) integration of control functionalities into components. As a matter of comparison, Julia [START_REF] Bruneton | An open component model and its support in Java[END_REF], the Fractal reference implementation, performs a load-time (LT) integration based on bytecode engineering performed with the ASM library [START_REF] Bruneton | ASM: A code manipulation tool to implement adaptable systems[END_REF]. While LT integration is somewhat more flexible than CT integration, CT integration leads to a solution which is more type safe and easier to debug.

Two versions of the aspect layers are available with AOKell: one based on AspectJ and one based on Spoon. AspectJ [START_REF] Kiczales | Getting started with AspectJ[END_REF] is the leading tool for AOP. A compiler is available and many tools, including plugins for IDEs, have been developed. While the most common usage of AspectJ is for weaving aspects at compile-time, AspectJ can also be used as a load-time weaver. Spoon [START_REF] Pawlak | Spoon: Annotation-driven program transformation the AOP case. 1st Workshop on Aspect-Orientation for Middleware Development @ Middleware'05[END_REF] is a general purpose tool for transforming Java programs. A strongly-type template mechanism is available for writing any kind of transformations including the ones (code advising and code injection) which are similar to the ones performed by aspect weavers. Spoon is architectured as a back-end of the javac compiler. The functionalities of the two versions are strictly equivalent. However, as we will show it in Section 4, the Spoon-based version of AOKell performs better than the AspectJ-based one. Indeed, being a general purpose tool for code transformation, Spoon allows generating code which is more aggressively optimized.

As an additional feature, AOKell allows developing the control layers as assemblies of component. For that, we introduce the notion of a control component. A control component implements a given control function (e.g. managing bindings, lifecycle, attributes, etc.) and is integrated into the application with an aspect. The idea is that controllers are seldom autonomous, but that, most of the time, they require the collaboration of other control components. When programming these controllers as plain old Java objects (POJO), the interaction schemes between controllers end up being dug into the code as references to other objects. This hinders modularity and leads to code which is poorly structured. The idea is then to apply to the control layer the same principles which have been applied to applications: capture the architecture of the control layers with an ADL (in our case Fractal-ADL), and implement control membranes as assemblies of control components. AOKell provides a set of 13 componentized control membranes, corresponding to the ones which can be found in Julia. The development of a new control membrane is then a matter of writing a new assembly of control components. Although the componentization of membranes is available with AOKell, this is not a mandatory feature: developers can still write controllers and membranes as POJOs.

FractNet

FractNet2 is the implementation of our reference model for the .NET platform. While the aspect layer is specific to FractNet, the control layer of AOKell is reused as this with FractNet. The Java implementation of this layer is compiled with the Visual J# compiler to produce a regular .NET assembly which can be linked with the rest of the application.

The aspect layer of FractNet is implemented with the AspectDNG [START_REF] Gil | AspectDNG. DotNetGuru[END_REF] aspect weaver. As-pectDNG performs a compile-time weaving of aspects on compiled .NET assemblies. This ap-proach has two main advantages: an application can be woven even if its source code is not available, and AspectDNG can weave applications written in any language supported by .NET (C#, J#, VB.NET, Managed C++, etc.). By this way, Fractal applications with FractNet can be written in any languages that can be compiled into a .NET assembly.

Shared code elements

The source code of AOKell and FractNet is organized in three main packages (see Figure 2): glue implements the aspect layer defined in the reference architecture, lib implements the control layer, component is a library of predefined components (this library contains the bootstrap component). A fourth package, tools, provides some utility tools which are not mandatory for building and running a Fractal application with AOKell or FractNet. As summarized in Table 1, the glue package is specific to AOKell and FractNet, whereas the component and lib packages are common to both implementations. These packages are written in Java. They are compiled with the standard javac compiler in the case of AOKell, and with the Visual J# compiler in the case of FractNet.

AOKell

FractNet The figures given in the last column of Table 1 correspond to the percentages in terms of lines of source code for each package. AOKell with Spoon contains about 10,000 lines of code. The glue package, which contains 12% of the source code in this case, is the only one to be specific to either AOKell or FractNet. The remaining 88% are shared between the two implementations.

Performance Measurements

This section reports on a microbenchmark which has been conducted with a simple application containing two components: a client component and a server component. The server component provides an interface with eight methods. Each method owns a different signature, either without parameters, or with primitive parameters, or with object references parameters, and/or with return types.

The measures are taken on 2Ghz Pentium 4 PC running Windows XP Pro, Sun JDK 1.5.0 and the Microsoft .NET Framework 2.0. A warm-up phase is performed before taking measures to avoid bootstrapping and class loading costs. The test consists of series of calls emitted from the client component to the server component. In table 2, the figures correspond to the times taken by 8,000,000 calls (1,000,000 per method defined in the interface provided by the server component). The given figures correspond to the average value of 4 runs. This microbenchmark provides two series of measures depending on whether the lifecycle controller is activated or not. As a matter of comparison, we also provide some reference costs in pure Java or pure C#, and with three interception techniques: AspectJ, JBoss AOP and statically generated proxy classes in Java.

The conclusion which can be drawn from this microbenchmark is that, in the case of AOKell, Spoon delivers better performance than AspectJ. The difference between AOKell and Julia mainly comes from some difference in the semantics of the lifecycle controller implemented by these two frameworks. The performances of the Java and .NET implementations of Fractal are roughly similar, the Java ones being a little more performant.

Conclusion

This paper presents two implementations of the Fractal specifications: one, AOKell, for programming Fractal applications with Java, and one, FractNet, for programming with the languages of the .NET framework. Compared to other existing implementations of Fractal, the novelty of our approach is to provide a reference model which defines some general design principles for deriving frameworks implementing the specifications. The originality of this model is to be based on concepts of aspect-oriented programming.

Three layers can be found in the model: the application layer, the aspect layer and the control layer. The implementation of the control layer is shared between the two implementations, leading to a high level of code reused: about 88% of the FractNet code comes from AOKell (written in Java, this code is compiled for .NET with the Visual J# compiler). The aspect layer is specific to each implementation. Two versions of the layer exist for AOKell, one written with AspectJ [START_REF] Kiczales | Getting started with AspectJ[END_REF], and another one written with Spoon [START_REF] Pawlak | Spoon: Annotation-driven program transformation the AOP case. 1st Workshop on Aspect-Orientation for Middleware Development @ Middleware'05[END_REF] (the latter performs better than the former). The aspect layer of FractNet is written with AspectDNG [START_REF] Gil | AspectDNG. DotNetGuru[END_REF]. The benefit of using AspectDNG is that any .NET assembly can be woven with the control layer. Thus, programming Fractal applications with FractNet can be conducted, not just in one language, but in any language which compiles code into .NET assemblies.

As a matter of future work, we plan to use our frameworks, both AOKell and FractNet, to build dynamic service-oriented platforms [START_REF] Escoffier | Developing an osgi-like service platform for[END_REF]. Another direction would consist in defining a SPI (Service Provider Interface) associated to the reference model. This SPI would provide a more concrete framework for reusing finer grained piece of code between different Fractal implementations. Concerning FractNet, we plan to migrate to another .NET aspect weaver since the AspectDNG project seems to be discontinued.

Figure 1 :

 1 Figure 1: High-level view of the reference model.

Figure 2 :

 2 Figure 2: AOKell and FractNet packages organization.

Table 1 :

 1 AOKell and FractNet code sharing.

	% (lines of source code)

Table 2 :

 2 Cost of invoking and executing an operation (x 8,000,000).

	Operation execution time
	no lifecycle	lifecycle

http://fractal.objectweb.org

http://www-adele.imag.fr/fractnet

Acknowledgments

The work on AOKell is partially funded by France Telecom under the external research contract #46131097.


Nicolas Pessemier
email: pessemie@lifl.fr

Lionel Seinturier
email: seinturi@lifl.fr

Thierry Coupaye
email: thierry.coupaye@rd.francetelecom.com

Laurence Duchien
email: duchien@lifl.fr

A Safe Aspect-Oriented Programming Support for Component-Oriented Programming

In this paper we show that Aspect-Oriented Programming (AOP) can be safely supported by Component-Oriented Programming (COP) by providing a way to control the openness of a component with regards to AOP techniques. Our proposal reconciles the intrusive nature of AOP with the "black box property" of components in COP. We build a compromise between modularity and openness applying the open modules approach to components. The experiment has been achieved on FAC, our model that unifies the notions of component and aspect. We show that most of open modules principles are directly available within our approach, we then study requirements for others. Once all these principles integrated, we are able to tune the accessibility of the content of a component to AOP during system runtime. Thus, components become grey boxes with dynamic variation points accessible to AOP techniques.

I. INTRODUCTION

Component-Oriented Programming (COP) proposes to enhance object-oriented programming by separating concerns into clearly defined entities, called components. Reusable components with contractually specified interfaces are defined and composed together [START_REF] Szyperski | Component Software: Beyond Object-Oriented Programming[END_REF]. Nevertheless, whatever the decomposition adopted to represent a system, it has been shown that some concerns are mixed within a same component (code tangling), and that some concerns are scattered across several components [START_REF] Duclos | Describing and using non functional aspects in component based applications[END_REF], [START_REF] Lieberherr | Programming with Aspectual Components[END_REF]. These concerns which are called crosscutting concerns hinder the reusability, the maintainability, and the evolvability of applications.

To tackle these issues, some approaches have proposed a support for Aspect-Oriented Programming (AOP) in component-based systems [START_REF] Lagaisse | Component-based open middleware supporting aspect-oriented software composition[END_REF], [START_REF] Mezini | Conquering Aspects with Caesar[END_REF], [START_REF] Suve | JAsCo: an aspect-oriented approach tailored for component based software development[END_REF]. AOP is a wellknown paradigm to overcome this issue by modularizing crosscutting concerns using aspects [START_REF] Kiczales | Getting started with AspectJ[END_REF]. The main issue of these approaches is that AOP is applied regardless of the components themselves, the aspects are woven on the objects which implement the components. This intrusiveness breaks the encapsulation property of components and consequently their implicit contracts. This appears to be a major issue in COP where contracts and encapsulation are fundamentals.

At the object level, solutions have been proposed to overcome the issue of intrusiveness of AOP. For example, Aldrich introduced the notion of open modules, a new module system to open a program to AOP while keeping modularity by hiding implementation details of the module [START_REF] Aldrich | Open modules: Modular reasoning about advice[END_REF]. A module is defined as a set of entities which share a set of access points for the join points (points in the program execution flow where aspects will apply) exported by the module. Using this module system, the content of a module can be preserved by designating only the variation points where aspects can act.

In this paper we propose to push the open modules approach a step further by applying it to COP. The objective is to provide a safe way to support AOP in COP by controlling the openness of a component with regards to AOP. This control over variation points of a component can be seen as a compromise between modularity and openness. Our study focus on an extended component model for components and aspects which is presented in [START_REF] Pessemier | A model for developing component-based and aspect-oriented systems[END_REF]. Our model, named FAC, unifies COP and AOP notions together by representing AOP notions as component ones. We show that when components and aspects are unified, some open modules properties are directly handled as first-class entities in our model. We have then extended our model to handle all open modules properties. A component can declare its variation points. Since our model is dynamic, this declaration can evolve at runtime. The remainder of this paper is organized as follows. Section II provides some background on our unified model for aspects and components and on the open modules approach. Section III shows how we have applied open modules to FAC. Finally, Section IV concludes.

II. BACKGROUND

This section provides some background on our previous work [START_REF] Pessemier | A model for developing component-based and aspect-oriented systems[END_REF] on the unification of AOP and COP, and introduces the principles of the open modules approach.

A. FAC: An unification of AOP and COP towards COP

Since COP fails in supporting crosscutting concerns [START_REF] Duclos | Describing and using non functional aspects in component based applications[END_REF], [START_REF] Lieberherr | Programming with Aspectual Components[END_REF], our motivation was to give a support to AOP in COP but also to take advantage from the strong encapsulation property of COP in AOP. Therefore, our proposal is built as a twofold integration of AOP and COP which has the benefit of representing AOP notions using COP ones. Thus, we introduce three main concepts which are related to the general notions of component, binding and composite-component that generally appear in COP [START_REF] Szyperski | Component Software: Beyond Object-Oriented Programming[END_REF].

• An Aspect component is the representation of an aspect as a component. It offers as a provided interface a piece of advice code (the additional behavior to weave on other components). Basically, an aspect component applies around incoming and outgoing calls on component interfaces. Because we represent an aspect as a component we call our approach symmetric. Aspects and components are components, and can interact together using bindings (Figure 1 represents an aspect component connected to other components using various types of bindings). Traditionally in AOP two dimensions are considered: the base and the aspect dimension. Aspects are woven on the base and the base if oblivious of the aspect dimension. In our approach these two dimensions are unified to facilitate the interactions between components and aspects and their evolution.

• An Aspect domain is the representation of the domain of action of an aspect. It is represented as a composite component which contains the aspect component and the components on which it is woven. The notion of aspect domain can be seen as a reification of the notion of pointcut in AOP. Usually a pointcut is a description of a set of join points on which an advice code is woven. In our model we reify this notion as a first class entity (a composite-component) which contains the aspect component woven, and all the components affected by the aspect component. This clarifies the domain of impact of an aspect component in a system. This explicit relationship between advised code and aspects is a notion currently missing in AOP. Figure 1 gives an example of the aspect domain of the transaction Aspect Component which is woven on components C, D and E.

• An Aspect binding is the representation of the implicit link which exists between an aspect and a component on which the aspect is woven. The aspect binding notion can be seen as a more fine-grained notion than the aspect domain to capture the interaction between a component and a particular aspect component. Our philosophy is to consider only one dimension (aspects are components) and two types of interactions (regular bindings and aspect bindings). Aspect bindings are used to connect an aspect component with a component. By this way, each component can locally manage the aspects applied on its incoming and outgoing interfaces.

We have successfully mapped this general model with its three notions to the Fractal component model. Fractal [START_REF] Bruneton | An open component model and its support in Java[END_REF] is a reflective and extensible component model, where bindings can be set and unset dynamically (at runtime); reflection is available through the use of special kinds of meta-interfaces called control interfaces.

We have extended this model by introducing our three main notions. This extension is called Fractal Aspect Component (FAC for short). We have used the provided notions of component and composite-component to represent our notions of aspect component and aspect domain. We have introduced a new control interface called the weaving interface. This interface, which appears on each component of the system, is in charge of setting/unsetting aspect bindings and of weaving aspect components. It has the benefit of locally managing the ordering of aspects for a component.

In addition to the advantages mentioned in the description of our three notions, the mapping onto the Fractal component model allows setting/unsetting aspect bindings dynamically. This makes our weaver dynamic and this opens the way to dynamic adaptation [START_REF] Pessemier | A three level framework for adapting component-based systems[END_REF].

The complete description of our general model for component and aspect and its mapping to the Fractal component model [START_REF] Bruneton | An open component model and its support in Java[END_REF], named FAC is beyond the scope of this paper and can be found in [START_REF] Pessemier | A model for developing component-based and aspect-oriented systems[END_REF]. The next sub-section introduces the open modules approach which allows controlling the degree of openness of a module to AOP.

B. Open modules approach

The concept of open modules has been introduced by Aldrich to limit the access to join points of a system, which are accessed intrusively in AOP. With this approach any exposed join points has to be declared within a module, a special set of classes, to be accessed by aspects. The open modules approach is defined by Aldrich as follows: "Open Modules describes a module system that:

• Rule 1 allows external advice to interactions between a module and the outside world (including external calls to functions in the interface of a module)

• Rule 2 allows external advice to pointcuts in the interface of a module

• Rule 3 does not allow external modules to directly advise internal events within the module, such as calls from within a module to other functions within the module (including calls to exported functions)." We intentionally add rule numbers on items to facilitate the discussion in the following sections. The complete description of this module system can be found in [START_REF] Aldrich | Open modules: Modular reasoning about advice[END_REF].

More recently, the concept of open modules has been applied to AspectJ [START_REF] Ongkingco | Adding Open Modules to Aspectj[END_REF]. In this study, authors have extended the concept with new interesting features. Most of them are specifically related to AspectJ in order to support AspectJ pointcuts. Nevertheless it seems to us that some can be generalized out of the context of AspectJ. Thus, the most important feature with regards to our application of open modules to COP is the ability to open and not only to reduce the visibility on join points of a module. This can become extremely useful when using for instance debug aspects, or when considering dynamic adaptation using AOP. Among new features provided by this study, an interesting one is the ability to designate to which a pointcut is exposed to using a pattern language based on package hierarchy.

To help the discussion of the next section we will call the ability to designate which aspect has access to a module Rule 4, and the ability to open a module by exposing join points Rule 5.

A. Similarities

The first obvious similarity is related to the notion of a module and a component. A module in the open modules approach is a collection of classes which share a set of access points to AOP. A component in COP is a contractually specified entity which provides and requires services by means of interfaces. A component is a black box which naturally hides its implementation details as required by Rule 3 of open modules. Given that in FAC, join points are incoming and outgoing calls on component interfaces, the definition remains correct with regards to Rule 1, 2 and 3. Rule 1 and 2 are understood in FAC by the fact that an aspect component only applies to client and server interfaces (Rule 2). Following the definition of Rule 1, an aspect component interacts between a module (component) and the outside world (other components). Rule 3 is preserved as soon as we do not want to break encapsulation in FAC. Join points are not points inside a component. However, when applying aspect component behavior on component external interfaces, we may consider that the original behavior of the component is altered by the aspect component. Thus, the behavior expected from a given component may be different. It seems important that the weaving of aspect component on provided and required interfaces of a component should be better controlled in order to provide a safer integration of crosscutting concerns. We elaborate more on that particular point in Section III-B.

The second similarity is related to Rule 4 which has been defined to the particular use of AspectJ but can be also used within our system. This rule allows to clearly designate which aspect can apply on a given exported join point of a module. A regular expression is given as a parameter of the expose to declaration which designates a set of packages that are authorized to access the module. In FAC we have a very similar notion: the aspect binding. An aspect binding is set between a component and an aspect component using the weaving interface of the component. Because the weaving interface is a kind of meta interface, we can consider that the access policies defined on components are of the same type of meta-informations than the ones corresponding to the expose to definitions used in the extension of AspectJ supporting open modules. This means that we can consider Rule 4 as naturally handled by each individual component which are able to choose the aspect to be impacted by.

At this point we have seen that Rule 1, 2, 3, and 4 are naturally handled by our model and its mapping to FAC. In the next section we study the requirements to manage remaining rule, Rule 5.

B. Dissimilarities

We have seen that in our model and in its mapping to Fractal, FAC, the considered join points are incoming and outgoing calls on component interfaces. Thus, Rule 3 is implicitely preserved when considering components as modules. However, this also means that the join points inside the component are not exposed. The original idea of open modules is to define some pointcuts and join points and make them available by means of interfaces of a module. In our case, the content of a component is implicitly protected, but we still need a support to give an access to other join points, i.e., join points inside a component. Nevertheless, associating the level of implementation of components and the level of interaction between components (more architectural) takes part in our global vision of what really means applying AOP to COP that we have exposed in [START_REF] Pessemier | A three level framework for adapting component-based systems[END_REF]. The unification of these two levels will allow us to look inside components and to externalize some join points. We believe that these internal join points must also be controlled by the weaving interface. Thus, we would be able to fulfill needs for accessing internal join points, while preserving Rule 3 by means of the weaving interface to control what is accessed or not. This approach has limitations with regards to legacy components that would not have be able to be instrumented to support AOP techniques. In our approach, we only consider a full-fledged component and aspect approach, where the design of the system follows the same formalism, the same design model. The remaining rule (Rule 5) states that the open modules should not be limited to reducing the access to aspects but also to opening it. In Section II-A we have defined the role of the weaving interface as an interface to manage the setting/unsetting of aspect bindings, the weaving and the ordering/re-ordering of aspect components. We have extended the role of this interface to also manage the openness of a component to AOP. A conceptual view of the role of the weaving interface is presented in Figure 2. The weaving interface is able to prevent the weaving of aspects on a particular operation of an interface. Because FAC is a fully dynamic framework, these policies can be changed during runtime. A component can then be adapted to open or reduce the access to its join points. Moreover, since we have discussed it in Section III-A, weaving an aspect component on component external interfaces may change the expected behavior from other components as soon as it intercepts communications between interfaces. The idea of controlling the external join points which are accessible or not by other components seems interesting even if it was not originally considered by Aldrich in the first version of open modules.

E A V I N G I N T E R F A C E C O M P O N E N T = M O D U L E I N T E R F A C E A I N T E R F A C E B I N T E R F A C E C I N T E R F A C E D i n c o m i

IV. CONCLUDING REMARKS

We have seen that AOP and COP can be reconciled on the particular issue of the intrusive property of AOP versus the strong encapsulation property of COP. To do so, we have applied the open modules approach to FAC, a unified model for components and aspects. Following the open modules philosophy our approach is able to open a component to AOP while keeping its content hidden from the outside. This compromise opens the way to a safe integration of AOP in COP. It is safe in the sense that the intrusiveness of AOP is finely managed on each individual component. Moreover in the case of FAC, we have seen that the openness of a component to AOP can be managed at runtime, allowing a component to adapt to unanticipated requirements. Black box components become grey box components providing variability points where AOP can access.

This integration of open modules approach to FAC takes part of our overall vision of applying AOP to COP presented in [START_REF] Pessemier | A three level framework for adapting component-based systems[END_REF]. This vision is based on three levels: (1) An architectural level which is achieved with FAC where aspects notions are mapped on component ones; (2) A control level where aspects can be used to inject the control level of components into the remaining level, (3) the level of implementation of components. In [START_REF] Seinturier | A component model engineered with components and aspects[END_REF] we show how AOP can be used for [START_REF] Bruneton | An open component model and its support in Java[END_REF]. In this paper we present a link between the architectural and the implementation level: Opening a component to expose internal join points and then, weaving architectural aspects to this internal join points (link between (1) and ( 3)). Our next step is to make all this three levels work together, allowing a coherent weaving of aspects whatever the chosen level.

Fig. 1 .

 1 Fig. 1. Aspect binding best practice

Fig. 2 .

 2 Fig. 2. A conceptual view of Open Modules applied to FAC. A component is a module which controls the access to its provided and required interfaces by means of the weaving interface. The content of the component is not exposed to aspects.

ACKNOWLEDGMENTS

This work was partially funded by France Telecom under the external research contract number 46 131 097.