key: cord-0004759-up0lh5n6
authors: Arino, Julien; Portet, Stéphanie
title: Epidemiological implications of mobility between a large urban centre and smaller satellite cities
date: 2015-01-14
journal: J Math Biol
DOI: 10.1007/s00285-014-0854-z
sha: 824394fc7732bb9838da13929fa0135921841a5b
doc_id: 4759
cord_uid: up0lh5n6

An SIR infectious disease propagation model is considered that incorporates mobility of individuals between a large urban centre and smaller satellite cities. Because of the difference in population sizes, the urban centre has standard incidence and satellite cities have mass action incidence. It is shown that the general basic reproduction number [Formula: see text] acts as a threshold between global asymptotic stability of the disease free equilibrium and disease persistence. The case of Winnipeg (MB, Canada) and some neighbouring satellite communities is then considered numerically to complement the mathematical analysis, highlighting the importance of taking into account not only [Formula: see text] but also other measures of disease severity. It is found that the large urban centre governs most of the behaviour of the general system and control of the spread is better achieved by targeting it rather than reducing movement between the units. Also, the capacity of a satellite city to affect the general system depends on its population size and its connectivity to the main urban centre.

more prevalent and easy. Among the consequences of this acceleration and generalization of mobility is a great speeding up and easing of the spatio-temporal spread of infectious diseases, as evidenced notably by the 2003 SARS epidemic or the 2009 H1N1 pandemic (Khan et al. 2009 ).

Studying the effect of mobility on the spread of an infectious disease can be accomplished in many ways, illustrating that the term mobility encompasses a wide variety of processes taking place on vastly different timescales. One of these methods are socalled metapopulation models, which consider distinct spatial locations called patches. In each patch, an epidemic model describes the spread of an infectious disease among members of the local population. The locations are then coupled in a graph, with the patches as vertices and arcs representing the possibility for individuals in the various epidemiological states to travel between locations.

Many aspects have been investigated in metapopulation models; see, e.g., the review by Arino (2009) . One aspect that has not been the object of much attention is the effect of the interconnection of patches with very different population sizes on the dynamics of disease propagation in the coupled system. Wang and Zhao (2008) consider a model for two patches with periodic coefficients of infection and movement rates. Fromont et al. (2003) consider Feline Leukemia Virus in a metapopulation of cats living either in farms or in villages. In these papers, the authors consider coupled patches with some having mass action incidence and others with standard incidence.

In the same spirit, a metapopulation model for the spread of infection within cities is used here to study the movement of people between an urban centre and neighbouring satellite cities and the effect this has on the spatial and temporal spread of the disease. Large cities are indeed often associated to smaller cities of two main types. Suburban areas are most of the time very close to the major city. They can have distinct administrative structures, but a large fraction of their inhabitants commute to the large city (or to other suburban areas) during the work week, if not every day. They also often offer some of the amenities of the large city, such as shopping centres or recreational areas, so inhabitants of the latter also visit them regularly. They are difficult to distinguish from the large city and from a modelling perspective, the large city and its suburbs can generally be considered as a single unit. The other type of structures peripheral to a large urban area are so-called satellite cities. These are generally cities in their own right, located further away from the main city than are suburban areas. However, they are connected to the main city in many ways. For instance, while they generally have hospitals, these facilities are often able to deal only with general cases, with severe cases having to be treated in the large city's facilities. Although not as closely integrated as a major city and suburban cities, individuals make regular movements between the locations.

An SIR-type metapopulation model is formulated, with standard incidence in the large urban centre and mass action incidence in the satellite cities, to tackle the following questions:

1. From Arino (2009) and other work on the subject, it appears that most of the time, the linear autonomous coupling of a collection of patches that are identical, save perhaps for parameter values, does not lead to behaviours that are more complex than those observed for the models in isolation. On the other hand, a change of incidence function within a single population can lead to complicated dynamics; see, e.g., Arino and McCluskey (2010) . Does the coupling of patches with different incidence functions affect the dynamics of a metapopulation disease transmission model? 2. From a more practical point of view, urban centre and satellite cities systems abound. In such configurations, what drives disease propagation in the entire system? Can satellite cities "infect" the urban centre, can they protect themselves from infection by the urban centre? Does the size of a city play a role in its vulnerability to invasion by an infectious disease or its role in facilitating the spread of this disease?

Some mathematical analysis is conducted. Then the model is specialized numerically to several cities in the province of Manitoba (Canada): Winnipeg, the capital and largest urban centre with about 700,000 inhabitants, and satellite cities located some distance from Winnipeg.

In each city, the population is divided between susceptible (S), infectious (I) and recovered (R) individuals, depending on their epidemiological status. The simple SIR formalism allows to capture the main characteristics of disease propagation without the burden of additional equations and parameters. Further, this allows to focus on the effect of size, since the behaviour of a classical SIR metapopulation model without size effects is well understood.

The large city is labelled with the numerical application to Winnipeg in mind: it is identified with the index W . There are a fixed small number n of satellite cities, labelled 1, . . . , n. It is assumed that there is no movement between the smaller cities. This is, to a large extent, a reasonable assumption in the application considered; see Sect. 4. The mathematical analysis is not any more complicated if movement between the satellite cities is added to the model, it is omitted here for simplicity and because of the application considered. All individuals move at the same rate regardless of their epidemiological status. This simplification stems from the short distances involved.

The specificity of the model lies in the assumption that, depending on the size of the population in patches, different types of incidence functions are present (Fig. 1 ). There is a lot of debate about the nature of incidence functions (McCallum et al. 2001) and the interpretation that follows is by no means the only one possible. Here, as in the work of Arino and McCluskey (2010) , it is assumed that mass action incidence is appropriate for smaller populations, while standard (or proportional) incidence better fits larger communities. Since the communities are much smaller than the urban centre, the nature of contacts is different. In the smaller communities, it is quite conceivable that an inhabitant can meet any other inhabitant, which is well described by a mass action incidence function. On the other hand, a large urban centre has several major shopping areas, numerous neighbourhoods, etc. Many inhabitants spend their days in a small number of neighbourhoods, rarely if ever visiting others. This situation is better reproduced by a proportional incidence function. 

The model consists of 3(n + 1) differential equations. Those for the urban centre are given by

while equations for the satellite cities take the form, for i = 1, . . . , n,

Here and throughout the remainder of the text, the index x ∈ {W, 1, . . . , n} is used to alleviate notation; when referring to a city, "city x" refers to any of the urban centre or the smaller cities. The parameter β x is the transmission coefficient in city x; it has units per unit time and per unit time per individual in W and i = 1, . . . , n and represents the rate of infection and the rate of infection per individual, respectively. The parameter γ is the per capita rate of recovery, i.e., 1/γ is the mean of the exponentially distributed time of sojourn of individuals in the infectious class. The constant b x is the birth rate in city x and d is the per capita death rate. The rates of recovery and death are equal in all cities because sanitary conditions are comparable in the cities considered. The birth rate is allowed to differ between locations as it is used to set the equilibrium value of the population in the locations (see below). The parameter m W i is the rate of movement from satellite city i = 1, . . . , n to the urban centre, while m i W is the rate of movement from the urban centre to satellite city i = 1, . . . , n. Unless otherwise indicated, the rates of movement are positive. All other parameters are positive.

The total population in city x is N x = S x + I x + R x . This number is assumed to be large compared to the number n of cities considered, so that an ordinary differential equations approach is justified. System (1) is considered with nonnegative initial conditions such that

To avoid dealing with a trivial case, it is also generally assumed that

Summing all equations in (1) and denoting N = x N x and b = x b x , the total population in the system is governed by

so the total population in the system tends to b/d as t → ∞. It is easy to verify that the positive orthant R 3(n+1) is invariant under the flow of (1), so all state variables are nonnegative and bounded. Some notation will be useful for the remainder of the analysis. If k is a vector, the notation k ≥ 0 indicates that k has all its entries nonnegative, k > 0 means that k ≥ 0 with at least one positive entry; finally, k 0 means that k is entry-wise positive. The same notation is used for matrices. The analysis will also make use of the movement matrix, which is defined as

Some properties of M will be used later on and are summarized in the following result.

Lemma 1 Let c ∈ R + \{0}. The matrix M defined by (2) has the following properties. Proof (1) follows from Theorem 2 in the paper by Arino (2009) . Since −M is an M-matrix, Lemma 6.4.1 in the book by Berman and Plemmons (1994) implies that, for c > 0, −M + cI is a nonsingular M-matrix, proving (2). Since −(M − cI) is a nonsingular M-matrix, Theorem 6.2.3.N 38 in Berman and Plemmons (1994) gives ( 

Consider first the uncoupled system, i.e., take M = 0. The analysis of the system in this case is well-known and summarized here for the reader's convenience; see, e.g., the papers by Korobeinikov and Wake (2002) and Vargas-De-León (2011) . In all cities, whether large or small, i.e., independent of the nature of the incidence function, the total population converges to N * * x = b x /d and the disease free equilibrium (DFE), obtained by solving for (S x , R x ) after setting I x = 0, takes the valuē S x = N * *

x andR x = 0. [The notation N * * x is used to avoid confusion with N * x in the general coupled model as given by (5), which includes mobility.]

Now compute the basic reproduction number R x 0 in city x using next generation matrix method of van den Driessche and Watmough (2002) . In the urban centre,

while for the satellite cities i = 1, . . . , n,

Supposing that I x > 0 gives the endemic equilibria (EEP)

in the large urban centre and

in the satellite cities.

It is then known that if R x 0 < 1, then the DFE is globally asymptotically stable in patch x, while if R x 0 > 1, then the endemic equilibrium is globally asymptotically stable in patch x.

Consider now the coupled system with M = 0. In the remainder of Sect. 3, it is assumed that movement rates from the urban centre to all satellite cities and from all satellite cities to the urban centre are positive. From Lemma 1 (v), it follows that M is irreducible.

From Lemma 1 (ii), −(M−dI) is nonsingular and from Lemma 1 (v), −(M−d I) −1 0. It follows that N * is unique and N * 0. Also, by Lemma 1 (iv), M−dI has all its eigenvalues with real parts less than or equal to −d, so solutions N of the linear system (6) converge to N * .

The proof of Theorem 1, later on, will proceed as in the paper of Li and Shuai (2009) and requires, in particular, the use of the positive invariance of the set

which is easily shown as a consequence of Proposition 1.

Let S = (S W , S 1 , . . . , S n ) T , I = (I W , I 1 , . . . , I n ) T and R = (R W , R 1 , . . . , R n ) T . A disease free equilibrium has I = 0. The following result holds true.

(1) has the unique DFE

which is easier to handle in vector form:

As previously, the matrix M − d I is invertible. Therefore,

It follows that R * = 0 and S * = N * .

The stability of the disease free equilibrium is now considered, using the general basic reproduction number R 0 , i.e., the basic reproduction number for the whole system. Using the method of van den Driessche and Watmough (2002), let

be, respectively, the vectors of new infections and all other flows within the infected classes. Taking partial derivatives of F and V with respect to I W , I 1 , . . . , I n and evaluating at the DFE gives

Let ρ(A) be the spectral radius of a given matrix A.

The DFE (7

The proof below is an extended outline showing in particular that the analysis can proceed as did Li and Shuai (2009) ; see that paper for details.

Proof The expression of R 0 , the local asymptotic stability and instability results are a direct application of Theorem 2 of van den Driessche and Watmough (2002) with the matrices F and V given by (8) and (9).

To prove the global asymptotic stability, proceed as in the proof of Theorem 3.1 in the paper of Li and Shuai (2009) with slight modifications. The matrices F V −1 and V −1 F share the same spectrum (Exercise 1.2.P17 in Horn and Johnson 2013) and

Thus, by Theorem 6.2.7 in the book by Berman and Plemmons (1994) 

, Theorems 1 and 2 in the paper by Berman and Shaked-Monderer (2012)]. Thus,

Differentiating L along solutions of (1) gives

However, since F is diagonal,

Substituting this into (12) and using (11), it follows that

Thus L ≤ 0 when R 0 ≤ 1, so L is a Lyapunov function for (1). Because of the invariance of Γ mentioned earlier, the remainder of the proof of Theorem 3.1 in the paper of Li and Shuai (2009) can be applied, giving the result.

Theorem 2 Suppose that R 0 > 1. Then the system (1) is uniformly persistent and there exists an endemic equilibrium.

Proof Consider system (1) as a nonautonomous system where N(t) is the nonautonomous component, i.e., "forget", for now, that N(t) = S(t) + I(t) + R(t). From Proposition 1, N(t) → N * , so (1) considered as a nonautonomous system is asymptotically autonomous with limit system the one where N(t) has been replaced by N * .

. . ,β n = β n . Substituting this into (1) gives, for all cities x ∈ {W, 1, . . . , n},

Note that m px and m xp are positive if and only if x = 1, . . . , n and p = W or x = W and p = 1, . . . , n, implying that (1d)-(1f) can indeed be written in the form (13). Once this step is accomplished, Proposition 3.2 of Li and Shuai (2009) holds. Indeed, (13) is a simpler version of (1.1) in the paper of Li and Shuai (2009) : here, movement and death rates are equal irrespective of an individual's disease status. Also, the movement matrix M is irreducible. Proposition 3.2 in the paper of Li and Shuai (2009) implies that (13) is uniformly persistent and has an endemic equilibrium (S * , I * ) in the interior of Γ -the dynamics of R is omitted at this stage as it does not influence S and I and can be deduced from N * , S and I.

The result then extends to system (1) by using the theory of asymptotically autonomous systems (Thieme 2000) .

Note that the last theorem in the paper of Li and Shuai (2009) , which establishes the global asymptotic stability of a unique endemic equilibrium, is not applicable here. Indeed, one of the major differences between the model of Li and Shuai (2009) and model (1) here is that movement rates depend on the epidemiological status of individuals in the work of Li and Shuai (2009) .

Finally, remark that because of the nature of F and V , the following result is obtained, which allows under conditions to localize R 0 when the individual R x 0 are known.

Suppose that for all i = 1, . . . , n, N * i = N * * i . Then the general basic reproduction number R 0 satisfies the following inclusion

In particular, if min{R W 0 , R 1 0 , . . . , R n 0 } > 1, then the DFE is unstable and (1) is uniformly persistent, whereas the DFE is globally asymptotically stable if

Proof V has all column sums equal to γ + d, i.e., 1 T V = (γ + d)1 T , with 1 T =  (1, . . . , 1) . From the latter form, it is clear that 1 T V −1 = 1/(γ + d)1 T . The matrices F V −1 and V −1 F share the same spectrum (Exercise 1.2.P17 Horn and Johnson 2013) and so R 0 = ρ(F V −1 ) = ρ(V −1 F). Right multiplication of V −1 by the diagonal matrix F amounts to multiplying each column of V −1 by the successive diagonal entries in F and so

from the hypothesis that N * i = N * * i for all i = 1, . . . , n. The matrix F V −1 is nonnegative by construction (van den Driessche and Watmough 2002), so it follows from Theorem 8.1.22 in the book by Horn and Johnson (2013) that (14) holds true. The remainder of the result is then a straightforward consequence of (14).

The assumption that N * i = N * * i , i = 1, . . . , n, is restrictive. However, from a modelling perspective, it is justified as mobility here represents short term two-way movement rather than migration. The population in the various cities can be interpreted as representing the carrying capacity of these cities. In an unconnected collection of cities, in order to reach this carrying capacity, the birth rate would beb i such that N * i = N * * i . This hypothesis is often used in numerics (and will be used in Sect. 4).

Winnipeg, the capital of the Canadian province of Manitoba, is the urban centre under consideration, while several satellite cities close to Winnipeg are the smaller satellite cities: Portage la Prairie, Selkirk and Steinbach. The situation of Manitoba is ideal for this study: population density is low outside of Winnipeg; the road network connecting Winnipeg to nearby cities is simple and well studied.

Most parameters take the same value in all patches: the cities are geographically close and there is no discernible difference between health systems for the disease under consideration. Disease related parameters are taken to loosely represent influenza. The mean duration of the infectious period is 4.1 days, which is consistent with parameters used in the influenza literature (Arino et al. 2006 ). The disease transmission coefficient β x is allowed to vary in order to obtain different values of R 0 in different locations.

Population counts used come from the 2011 Canadian census and are those of the cities themselves, not their metropolitan areas. See Table 1 for a list of the cities considered, their population, distance from Winnipeg and an estimation of the average daily number of individuals travelling between Winnipeg and these cities. Transportation data originates from the University of Manitoba Transport Information Group (UMTIG) and gives the average daily number of vehicles at several counting stations along major axes between these cities. Estimation of the average daily number of travellers is explained in Appendix 1. The indices used in simulations are shown next to the city names As a consequence, there remains to find values for the birth and movement rates. This is done in two steps: first, the movement rates are computed to reflect the actual movement of individuals between locations estimated from the transportation data. Then the birth rates are computed so that, in simulations, the system remains at an equilibrium with population in all cities equal to their values in the census.

To estimate the movement rates, consider city x and its population N x . Assume that the rate at which individuals leave city x to go to city y is m yx . Thus, ceteris paribus,

where T yx is the number of individuals going from x to y each day. It follows that

This is computed for all pairs (W, i) and (i, W ) of cities considered. Given knowledge of the movement rates and population numbers in the different cities, b is then set so that the latter are conserved given the former. From (5) As a consequence, starting with N(0) = N * allows to conduct numerical simulations at the population equilibrium so that the only effects visible are those due to the disease.

Note that this method can lead to some entries in b taking negative values: to keep the population constant in a community that sees a large movement-induced inflow of individuals and a lower outflow rate, the birth rate might end up negative. This, in turn, could lead to the corresponding S x (t) becoming negative for some t, since the right hand sides of equations (1a) or (1d) would be negative when S x = 0. This is only a potential issue when R 0 > 1, since if R 0 < 1, then the S x are known to converge to N * x (Theorem 1). In any case, this method is only used in numerical simulations in situations where the total population in patches is at equilibrium. Figure 2 shows the sensitivity of the general R 0 to changes in the value of R W 0 and R i 0 , i = 1, 2, 3. In each case, "with disease" means that R x 0 = 1.5, "without disease" means that R x 0 = 0.5. The R x 0 that is made to vary varies from 0.5 to 3. The cases with "disease in Winnipeg" have R W 0 = 1.5, two of the satellite cities with R i 0 = 0.5 and the remaining small city's R i 0 varying; the variation of R W 0 is not considered in this case. The case with "disease in satellites" has R W 0 = 0.5 and two of = 1, 2, 3 with R i 0 = 1.5, with the third one varying. For instance, the first three results (Winnipeg) show the sensitivity of the general R 0 to variations of R W 0 between 0.5 and 3, when the three small cities have R i 0 = 1.5 (first and last boxes) or R 0 = 0.5 (middle bar). Each box and corresponding whiskers are the result of 10,000 simulations.

In Fig. 2 , it can be seen that, with the parameters used, most of the variation of R 0 is driven by Winnipeg; the small cities have little influence, albeit varying: Portage la Prairie can trigger a bifurcation even when R W 0 < 1, while Steinbach has a more limited influence and variations of R 2 0 in Selkirk have almost no effect on the value of the general R 0 . 

Understanding the roles of the satellite cities requires to consider not only their size, but the rates of movement that connect them to Winnipeg. From Table 1 , daily travellers to and from Winnipeg represent approximately 32, 81 and 55 % of the respective populations of Portage la Prairie, Selkirk and Steinbach and these percentages are reflected in the actual movement rates.

As remarked in Fig. 2 , variations of R 2 0 in Selkirk have almost no impact on variations of the general R 0 , regardless of the scenario considered for R x 0 in the other cities. Selkirk is the closest satellite to Winnipeg and has the highest ratio of trips to inhabitants and thus, the highest movement rates. The average time an individual spends in Selkirk (their residence time there) is small. In a sense, this community is barely distinguishable from Winnipeg and additional sensitivity analysis (not shown) indicates that it does not really matter whether one considers it as part of Winnipeg or separated.

Portage la Prairie and Steinbach have comparable populations. However, with the parameters under consideration and with reasonable values of the R x 0 , it can be seen in Fig. 2 that only the former can cause the general R 0 to take values larger than 1 when R W 0 < 1. Intuition indicates that this must be due to the movement rate. If all movement rates were zero, then R 0 would be given by

since the matrix F V −1 is diagonal when M = 0. As movement rates to and from Portage la Prairie are lower, the situation is closer to the uncoupled case and the value of R 1 0 has more impact on the general R 0 . This is confirmed by Fig. 3a , which shows the joint effect of a reduction of movement between Winnipeg and Portage la Prairie and of the value of R 1 0 in Portage la Prairie, when R W 0 = R 2,3 0 = 0.5. On this figure, R 0 < 1 is located above the green curve. When Portage la Prairie is completely isolated from Winnipeg (bottom of the figure), the value of R 1 0 governs the value of the general R 0 . This effect decreases as the movement rates increase.

Although Portage la Prairie can drive the value of the general R 0 , its effect on other aspects of disease dynamics is not as important as that of Winnipeg. Figure 3b shows the attack rate in Winnipeg as a function of the same parameters as those in Fig. 3a . Attack rates are a fundamental concept in epidemiology that are used to quantify the Fig. 3 Plots as functions of R 1 0 in Portage la Prairie and the reduction of movement between Winnipeg and Portage la Prairie. a General R 0 ; the green and red curves show the locations of R 0 = 1 and R 0 = 1.5, respectively. b Attack rate in Winnipeg; the green and red curves show the locations of R 0 = 1 and a 10 % attack rate, respectively burden of an epidemic over its time course. The term "rate" is a misnomer, as most of the time, attack rates refer to the percentage of individuals in a population that have borne the pathogen during the epidemic; however, this is the most commonly used term to describe this concept and it is therefore used here.

In order to obtain attack rates as in Fig. 3b , the following approximation method is used. Numerical solutions to the ODE are computed over a time period. The number of infection-days for a given community x is then obtained by trapezoidal integration of I x (t). This is converted in turn to a number of infections by dividing by the average duration of the infectious period (4.1 days). Finally, this is related to the population N *

x and converted to percentage. In order to obtain comparable results, all attack rates graphs are produced here using simulations for 1 year.

Considering Fig. 3b , one sees that in most cases, even in situations where R 1 0 causes the whole system to go to an endemic equilibrium (below the green curve), the attack rate remains low in Winnipeg. To obtain attack rates over 10 % in a "normal" movement context, R 1 0 must be quite high, well over 2 (red curve). To put things in perspective, most influenza epidemics, whether annual or pandemic, have estimated R 0 values around 1.5. Pandemic planning scenarios prior to the 2009 H1N1 pandemic made assumptions placing the attack rate between 15 and 35 % for mild and severe scenarii, respectively; see Sect. Two in Public Health Agency of Canada (2006) . Interestingly, the higher values of R 0 are obtained in regions of parameter space where the attack rate in Winnipeg is small (lower right corner). Figure 3a , b also illustrate another interesting property of the system: while R 0 appears to depend linearly on R 1 0 , m 1W and m W 1 (green curve in both figures), the attack rate is more complex (red curve in Fig. 3b) . Suppose for instance that R 1 0 were a little larger than 2.5, placing the system in a region where R 0 > 1 but with an attack rate Fig. 4 Plots as functions of R W 0 (Winnipeg) and R 1 0 (Portage la Prairie). a Value of the general R 0 ; the green and red curves show the location of R 0 = 1 and 1.5, respectively. b Overall attack rate; the green and red curves are 5 and 30 %, respectively less than 10 % in Winnipeg. Then reducing the rate of movement between Winnipeg and Portage la Prairie to 50 % of its usual value would result in a marked increase of the attack rate. It is only by bringing movement to less than 10 % of its default value that the attack rate in Winnipeg would be made smaller.

Continuing on the relative roles of the urban centre, Winnipeg, and its least connected satellite, Portage la Prairie, it can be seen in Fig. 4 that the effect of raising the value of R 1 0 is much less than the effect of raising the value of R W 0 , both in terms of effect on the general R 0 (Fig. 4a ) or in terms of the overall attack rate (Fig. 4b) . If only one of R W 0 or R 1 0 is allowed to vary (moving horizontally or vertically in the figure) , then a reduction of the attack rate from 30 to 5 % requires little change to R W 0 (provided R 1 0 is not too large), whereas a lot more reduction of R 1 0 is needed to achieve the same effect. Similarly, reducing the value of the general R 0 to a value less than 1 is only achievable by acting on R 1 0 if R W 0 is already less than 1. Note that the situation is similar if R i 0 , i = 1, 2, 3, is made to vary instead of just R 1 0 (not shown).

The question of "control" is now considered. In this model, there are two separate mechanisms that can be used to control the epidemic: the value of R x 0 can be lowered or mobility can be restrained. In order to evaluate the relative value of such measures, the following experiments are conducted.

In Fig. 5 , R i 0 = 0.5, i = 1, 2, 3, R W 0 ∈ [0.5, 3] and the rate of movement to and from Winnipeg varies from 0 to 100 % of its standard value. Changing the value of R i 0 to, Fig. 5 General R 0 (a) and attack rate in the satellite cities (b) as a function of R W 0 in Winnipeg and the reduction in movement rate between Winnipeg and the satellite cities. The satellite cities have R i 0 = 0.5, i = 1, 2, 3 say, 0.9 (not shown), does not change the situation: the effect of a reduction of the movement rate is almost null, control of the general R 0 is only achievable through control of R W 0 , as can be seen in Fig. 5a . On the other hand, as shown in Fig. 5b , reducing the rate of movement between Winnipeg and the satellite cities to a very low value can help lower the disease attack rate in the satellite cities when R W 0 is large.

In Fig. 6 , R i 0 ∈ [0.5, 3], i = 1, 2, 3 and the rate of movement to and from Winnipeg varies from 0 to 100 % of its standard value. Comparing Fig. 6a , where R W 0 = 0.5 and Fig. 6b , where R W 0 = 0.9, one sees that when the reproduction number in Winnipeg is larger, the amount of effort required to bring the general R 0 to a value less than 1 increases, as the value of R i 0 at which R 0 = 1 is smaller in Fig. 6b . Here again, as in the case with only Portage la Prairie and Winnipeg, decreasing the movement rates results in an increase of R 0 . However, as seen in Fig. 6c , a decrease of the rates of movement has a much more complicated effect on the attack rate in Winnipeg. While a moderate decrease would be detrimental to the attack rate, a radical decrease of movement rates would lower the attack rate.

An SIR metapopulation model for the spread of an infectious disease between a large urban centre and smaller neighbouring satellite cities was considered. Because of the widely different population sizes, different incidence functions were used in the urban centre (standard incidence) and the satellite cities (mass action incidence).

Note that in the present work, death and recovery rates were taken equal in all locations. This is appropriate for the application to Winnipeg and neighbouring satellite cities. However, several communities in Northern Manitoba are in a configuration that closely resembles that of satellite communities: they are accessible only by air or (during the winter) using ice roads. Per capita air travel rates from these communities to and from Winnipeg is much higher than the average travel rate between larger sized Canadian cities. For instance, using the method of computation of travel rates explained earlier and air transportation data from the BioDiaspora Project (see, e.g., Arino and Khan 2014) Several of these communities have no resident doctor. In this context, the rates of recovery and death could be assumed to differ.

The first objective of this work was to study mathematically the effect of coupling together cities with different functional forms for incidence. The total population in each patch was shown to converge to a positive value. An expression for the general basic reproduction number R 0 was derived. The global asymptotic stability of the unique disease free equilibrium when R 0 ≤ 1 was established. The proof followed that in the paper of Li and Shuai (2009) . Interestingly, the Lyapunov function used by Li and Shuai (2009) for a system with all incidence functions of the same type can be adapted easily to the present heterogeneous case. The uniform persistence of the system as well as the existence of an endemic equilibrium were shown by setting other results by Li and Shuai (2009) in an asymptotically autonomous setting. It is suspected that when R 0 > 1, there is a unique globally asymptotically stable endemic equilibrium. However, here the method of proof used by Li and Shuai (2009) was not applicable because of restrictions on the movement matrices there and more work will be needed to prove this. It was finally shown that under some conditions, the general R 0 is bounded below and above by the minimum and maximum values of the individual R x 0 in the cities. The second objective of this work was to investigate a real life situation. For this, Winnipeg and its three main satellite cities were used with parameters appropriate for influenza. One possible mechanism for deriving rates of movement between cities was discussed, using estimates of the average daily number of car trips between Winnipeg and the three smaller cities that were obtained in Appendix 1. It was observed that the large city, Winnipeg, governed for the most part the behaviour of the system. Only Portage la Prairie, which has the smallest rate of movement to and from Winnipeg, was able to trigger an epidemic if all other cities (including Winnipeg) had a low value of R x 0 = 0.5. It was observed that, as far as the satellite cities are concerned, reducing the value of R W 0 was much more important than isolating from Winnipeg, since only almost perfect isolation was able to curb the overall attack rate of the disease in the case that R W 0 > 1. These numerical results provide a glimpse into the dynamics of infectious diseases in a satellite-central hub setting and highlight the major role played by the main urban centre on the course of an epidemic in such a system. They also emphasize the fact that not all satellite cities play equal roles, with the three cities considered here having different capacities to influence the dynamics of infection in the overall system, depending on the nature of their connection to the main urban centre.

Note that when it holds, the bounds given by Proposition 3 allow, for system (1), to answer in the negative an often difficult question in metapopulation models: can movement "make it worse", in the sense that the general R 0 would take values larger than any the individual R x 0 ? Numerical simulations were conducted in a context where this proposition does hold, highlighting one issue with numerics: because the total population is asymptotically constant in each patch and that simulations are conducted at that equilibrium, the effect of having different functional forms for incidence is very small. Numerical simulations (not shown) were conducted which show that, with the parameters used, one obtains virtually the same plots as in Fig. 4 if mass action or proportional incidence are used in all cities.

Finally, note that Figs. 3, 5 and 6 tell a cautionary tale: R 0 is, in the context of mathematical models such as the one here, a bifurcation parameter that allows to distinguish between two dynamical regimes. However, R 0 does not give the whole story. Take for instance Fig. 3 and suppose R 1 0 = 2.5 and the movement rate is at 30 % of standard, i.e., a point slightly below the red curve in Fig. 3a . Suppose that the only control measure available is to act on the movement rates. If one only considers Fig. 3a , then the obvious solution is to encourage more travel. However, Fig. 3b suggests that decreasing the movement rate further would be a much better strategy.

In view of the remarks above, it is clear that the ordinary differential equations framework used here is somewhat limiting and that further work should consider a stochastic approach such as the one used by Arino et al. (2011) ; Arino and Khan (2014) or Lindholm and Britton (2007) . The latter considers a different framework, with an implicit description of movement, but shows that a derivation of some important quantities such as the mean time to extinction of the infection is probably possible in the stochastic equivalent of the model considered here.

Contrary to data on air (Arino and Khan 2014) or train transportation, UMTIG transportation data is not usable directly in the model as it provides flows through points, not trips between locations. Among the vehicles counted at a given counting station, it is not possible to distinguish the ones making a trip between the cities of concern from the ones making either a longer or a shorter trip. Therefore, the numbers obtained here are an approximation.

Observe that the volume of movement on a given route between two locations can be no larger than the minimum of the values obtained at the counting stations on that route between the two locations. So the volume obtained is a loose upper bound. Canadian Vehicle Survey data for 2009 shows that the average vehicle occupancy in Manitoba was 1.65 persons per light vehicle; see Fig. 12 in the report by the Office of Energy Efficiency (2009). Numbers could be further adjusted to take into account this fact. However, due to the already present uncertainty, the raw estimated number of vehicles is used. Lastly, it is assumed that individuals usually try to optimize the length of their trips (whether in distance or in time) and thus avoid complicated alternate routes. Therefore, only the most direct routes between locations are considered.

Only the estimation of the volume between Winnipeg and Portage la Prairie is shown here. The other two proceed similarly. In practice, the traffic inbound from Winnipeg is captured by stations 368 and 2077 (see Fig. 7 ). The last available date for station 368 is 1993. In order to get a more recent estimate, note that traffic through stations 48 and 369 increased an average 31.5 % in the eastbound direction and 37.5 % in the westbound direction, from 1993 to 2011. Adjusting the 1993 estimates for station Fig. 7 Schematization of the location of the vehicle counters used to estimate traffic between Winnipeg (WPG) and Portage la Prairie (PLP). The numbers in disks refer to UMTIG counting stations, those in diamonds are road numbers 368 similarly gives 3,840 vehicles eastbound and 4,015 westbound. Station 2077 has a 2011 estimate of 5,430 vehicles. This is the total of northbound and southbound traffic. Because no further information is available, assume that traffic is split evenly between northbound and southbound vehicles, thus there are 2,715 vehicles travelling in each direction daily. Consider northbound vehicles. They may come from Highway 1 from both directions, or from further south on Highway 240. The latter are captured by station 366, which shows an undirected daily average traffic of 4,680 vehicles, i.e., 2,340 in each direction. Thus, the daily northbound flux of "nonresidents" is estimated to be 375 vehicles. Again, this is split evenly between traffic incoming from Highway 1 westbound and eastbound, giving a total of 188 vehicles coming from Winnipeg. The same volume is estimated to go to Winnipeg each day. Adding the volumes gives a total number of travellers to Winnipeg of 4,028, while there are 4,203 travellers from Winnipeg. To simplify a little further given the already present uncertainty, the average (rounded below) of these two values is used in both directions, giving the number in Table 1 .

Modeling and dynamics of infectious diseases

Using mathematical modelling to integrate disease surveillance and global air transportation data

Effect of a sharp change of the incidence function on the dynamics of a simple disease

Simple models for containment of a pandemic

Some methodological aspects involved in the study by the Bio. Diaspora Project of the spread of infectious diseases along the global air transportation network

Nonnegative matrices in the mathematical sciences

Disease propagation in connected host populations with densitydependent dynamics: the case of the Feline Leukemia Virus

Spread of a novel influenza A (H1N1) virus via global airline transportation

Lyapunov functions and global stability for SIR, SIRS, and SIS epidemiological models

Global stability of an epidemic model in a patchy environment

Endemic persistence or disease extinction: the effect of separation into sub-communities

How should pathogen transmission be modelled?

Tech Rep M141-18/2009E-PDF Public Health Agency of Canada (2006) The Canadian Pandemic Influenza Plan for the Health Sector

Uniform persistence and permanence for non-autonomous semiflows in population biology

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

On the global stability of SIS, SIR and SIRS epidemic models with standard incidence

Threshold dynamics for compartmental epidemic models in periodic environments

Acknowledgments This work was supported in part by NSERC. An earlier version of the model was studied by Lindsay Wessel as part of a MITACS-CDM funded undergraduate research internship at the University of Manitoba. Transportation data was obtained from the University of Manitoba Transport Infrastructure Group (UMTIG). The authors acknowledge suggestions from an anonymous referee and an editor that helped improve the quality of the manuscript.