1 Introduction

Since the beginning of man, infectious diseases have played an important role in how we interact, migrate, and form societies. Diseases that are highly infectious, can soon develop into an epidemic. A “global epidemic”—often referred to as a pandemic—is a disease that cause infection in a significant portion of a population over a large geographical area. Fortunately, most of the common infectious diseases are not lethal, but the ones that are can result in huge numbers of deaths, often significantly reducing, and even eradicating, the population in certain areas.

Historically, the world has experienced several pandemics, often with high death tolls; The plague of Justinian (541–542AD) struck the Eastern Roman Empire and is estimated to have had a mortality of 50 % to 60 %; the Black Death of the 14th century reduced the world population from an estimated 450 millions to between 350 and 375 millions; the tuberculosis pandemic of the 19th century killed about 25 % of the adult European population; and the Spanish flu (influenza) of 1918 is estimated to have caused the death of between 25 to 50 million people. More recently, the swine flu or the H1N1 influenza (a new variant of the Spanish flu) killed about 300,000 people within the first 12 months in 2009 [1]. In 2002, Hong Kong experienced an outbreak of the severe acute respiratory syndrome (SARS) which, according to the World Health Organization, nearly became a pandemic with 8,422 cases and 916 deaths worldwide between November 2002 and July 2003 (10.9 % fatality) [2]. SARS spread within weeks from Hong Kong to 37 other countries.

Needless to say, there is a huge interest in epidemic outbreaks, spreading and containment which has attracted scientists from a number of different disciplines. Early pioneers in infectious disease modeling were William Hamer and Ronald Ross (Nobel laureate in medicine 1902) who in the early twentieth century applied the law of mass action to explain epidemic behavior of common infectious fevers of childhood and malaria, respectively [3, 4]. The early work by Kermack and MacKendrick [5] on disease dissemination in a homogeneous population also deserves acknowledgment in this context. Since then numerous works have appeared of increasing complexity and detail. One set of models concern spatial dynamics dealing with epidemic spreading over geographic areas [6]. Examples range from forest-fire models [7] to spatial immunization patterns where multiple diseases compete for their hosts [8], and explicit simulation of 85 million people including households, schools and workplaces in Thailand [9], a region of high interest due to the avian H5N1 epidemic.

Another class is metapopulation models. Such models are based on the observation that large human populations are typically divided into a network of smaller subpopulations such as cities or villages. Crudely speaking, a metapopulation consists of a group of connected subpopulations. Between the subpopulations there is an ongoing exchange of people, some of which are carriers of a disease and may infect previously unaffected subpopulations [1014]. The dynamics in each “node” is often assumed to be governed by mass-action kinetics, or sometimes stochastic (relevant for very low population densities). Obviously, the way in which the subpopulations are interconnected has a significant effect on how the disease spreads through the metapopulation. Particular emphasis has therefore, on good grounds, been paid to the impact of airline travel using explicitly, or features of, the world’s airline transportation network [1519]. Recent efforts also try to improve the description of human mobility patterns such as recalling individuals’ geographic origins [20, 21].

A third model class, albeit related to the former, is network models that account for how humans actually meet [2224]. In metapopulation models it is often assumed that the population within a single subpopulation is well mixed, which means that every individual is equally likely to spread the disease to any other member of the population. This is of course not generally true. For example, in large cities any individual only meets a very small fraction of the total population. There is, however, frequent encounters between those belonging to the same social network.

There are also interesting works focusing on the situation where humans are not the main carriers of the disease. Examples include Bubonic plague carried by rats [25] and water-borne cholera epidemic in a South African river network [26].

In this paper, we use a metapopulation description to model disease spreading in a network of cites and villages (nodes) of varying population sizes which are interconnected via a transportation network. The infection model we assume is the well-established susceptible-infected-susceptible model (Sect. 2) where individuals can be in either a susceptible or an infected state and going from one state to the other is controlled by rates. For an isolated city or village, the fraction of infected individuals in the population in the final (stationary) state is fully determined by these rates. However, when the cities and villages in addition are exchanging people with each other, for instance on a transport network, the situation may be much more involved. Under such conditions, the final state is not only determined by the rates, but also depends on structural features of the network [10, 11].

The final state of the metapopulation can be either epidemic or non-epidemic, meaning that a fraction of the population is infected (or not). The separation between these two phases (in parameter space) is known as the epidemic threshold. In this study we go beyond previous works and address the importance of city-size heterogeneity on the epidemic threshold. The size of the city or village is an important aspect because it sets the limit on how large fraction of the population a single individual actually meet within a given time period. The well-mixed approximation works poorly for a large city like New York whereas for a small countryside village it could be adequate. We therefore let the infection dynamics at each node change depending on how many individuals are present at any given time. In previous models, the network’s topology only influenced how people, and thus the disease, spread between the nodes. Here we in addition let the dynamics adapt to the local instantaneous population density, which means that the network structure influence the disease dynamics in two ways. First, how individuals move, and second, how they interact on the nodes. The latter part is in practise achieved by introducing a population dependent infection rate that interpolates between the two limiting cases of “the big city” for which the infection rate is dropping with increasing population, and “the small village” where the infection rate remains independent of population size. Within this setting, we derive new analytic expressions for the epidemic threshold of a general network which interpolates between previously known cases.

This paper is organized as follows; Sect. 2 is devoted to detailing the infection model that will be assumed. Here also the behavior of this model in homogeneous space is discussed. We then present how we model the transportation of individuals between the cities and villages (Sect. 3) before we in Sect. 4 present a unified reaction-diffusion model that combines the two phenomena from the previous two sections. Section 4 also presents analytic results for the epidemic threshold and compare them to simulations. The conclusions from this study are presented in Sect. 5.

2 The Susceptible-Infected-Susceptible Epidemic Model

The model that we will focus on in this work is the classic susceptible-infected-susceptible (SIS) epidemic model [27]. This model is chosen since it is relevant for an important group of infectious diseases. Moreover, it has been studied extensively in the literature, and it can often be given an analytic treatment.

In the standard SIS model there exists two types of individuals—the individuals that are susceptible to an infection (labeled S below), and those individuals that already are infected and can infect others (labeled I). If an infected individual meets an uninfected individual, he will transmit the disease at rate β. On the other hand, an infected individual can recover from the illness at rate μ and thereby again becoming susceptible to the disease. Formally, the reactions of the SIS model can be written in the form

(1a)
(1b)

where we explicitly have indicated the infection (reaction) rate β≥0, and the recovery (healing) rate μ≥0. The SIS model is the prototype for a disease spreading model where individuals that have recovered from a disease will automatically be susceptible to re-infection. On the other hand, if individuals that recover from the disease become immune to it, the classic (prototype) epidemic model describing this situation is the related susceptible-infected-removed (SIR) model [27]. This latter model, however, will not be discussed further in this work.

The first property to note for the dynamics of the SIS model (1a), (1b) is that it conserves the total number of individuals, or in the language of physics, the total number of “particles”. This number we denote by

$$ N=N_S(t)+N_I(t), $$
(2)

where N p (t) is the total number of individuals of type p=I,S at time t. Note that in writing Eq. (2), we have explicitly indicated a time-dependence for N p (t) since the number of individuals of type p may fluctuate in time even if the total number of individuals N in the system remains constant at all times. The second feature of the SIS dynamics is the active role played by the infected individuals (the “I particles”); if they have disappeared completely at time t 0, i.e. N I (t 0)≡0, there is no way in the model to reintroduce the infection and, therefore, the population remains uninfected for all later times tt 0.

2.1 SIS Dynamics in Homogeneous Space

In order to understand the dynamics of the SIS model, we will start by studying it in homogeneous (continuous) space under the assumptions that susceptible and infected individuals are well mixed. Based on reaction (1a), (1b) we consider the model (in discrete time)

(3a)
(3b)

where the second and third term on the right-hand-side of these equations are associated with recovery and infection, respectively. For later convenience, we have in writing Eqs. (3a), (3b) introduced a characteristic population size, N ×, so that the effective infection rate is

(4)

From Eq. (4) it follows that in the limit NN ×, the infection rate β′ is independent of the total population size N. This situation is typical for a smaller community where one, within a rather short period of time, may be “exposed” to almost the entire population of the community. For instance, this may happen during a funeral, a music gathering or a Christmas party. For this reason, we will below refer to this type of interaction being of “the small village” type. In big cities, on the other hand, the probability per unit time to meet, and therefore potentially infect, another individual is reduced as the population of the city increases. In this situation we have chosen to model the effective infection rate with Eq. (4) so that β′∼N −1 in the limit where NN ×. In passing, we note that the “small village” and “big city” limits of Eq. (4) essentially are equivalent to what was called type I and type II dynamics in Ref. [10]. The only difference is that the infection rate (β) used for type I dynamics in Ref. [10] is in our case scaled by the characteristic population size, i.e. β′=β/N ×.

Empirically it is known that the city population distribution is a broad distribution that often can be characterized by an inverse power law—the so-called Zipf’s law [28]. This raises the question as to how one should choose the characteristic population size N ×, which separates the dynamics of a small village from that of a big city. Above we argued that the infection rate in a small village becomes independent of population size due to the likelihood that an infected individual potentially can infect everyone else within finite time. How likely this is depends on the social dynamics in the village, but we find it unlikely to happen if the population of the village is larger than, say, N≈103. Hence, we suggest to choose the characteristic population size as N ×∼104 or at most N ×∼105. Ideally one would have liked to incorporate the actual population size distribution into the dynamics (see Eq. (4)). However, doing so would have made the model more complex and an analytic analysis of its stationary state (see Sect. 3) would have become more challenging. Hence, we have not explored this possibility further in this work.

2.2 Stationary State

Under stationary conditions N I (t) and N S (t) are independent of time t, which we denote by \(\bar{N}_{I}\) and \(\bar{N}_{S}\), respectively. From Eqs. (3a), (3b) it follows that in a stationary state, the following relation must be satisfied

$$ \mu\bar{N}_I = \beta\frac{\bar{N}_S \bar {N}_I }{ N_\times+ N}, $$
(5)

with the additional constraint (see Eq. (2)) that \(\bar{N}_{S} + \bar{N}_{I}=N\). Equation (5) is trivially fulfilled for \(\bar{N}_{I}=0\) (and \(\bar{N}_{S}=N\)), meaning that the disease has died out and will never reappear. However, Eq. (5) is also satisfied if \(\bar{N}_{S}=\mu/\beta'\) with the requirement that \(\bar{N}_{I}=N -\bar{N}_{S}\geq0\). The (non-trivial) condition reads

$$ \bar{N}_I = N \biggl[ 1- \frac{\mu}{\beta} \biggl( 1+ \frac{N_\times}{N} \biggr) \biggr] > 0, $$
(6)

which holds if

$$ \frac{\beta}{\mu} > 1+ \frac{N_\times}{N}. $$
(7)

This means that an “epidemic phase” exists (\(\bar{N}_{I}>0\)) in the long time limit when individuals are infected faster than they recover. In all other cases, the only acceptable solution to Eq. (5) is the non-epidemic solution \(\bar{N}_{I}\equiv0\) and \(\bar{N}_{S}\equiv N\). Since the border between these phases occurs when β/μ=1+N ×/N, this point in parameter space is referred to as the epidemic threshold [27].

3 Migration Within a Metapopulation

In the previous section, we considered the dynamics of the SIS model under the assumption of an isolated, homogeneous well mixed population. It was shown that an epidemic phase may exist. However, humans are not uniformly distributed over the surface of the earth. Instead, throughout history, we have had the tendency to form local communities (“subpopulations”), i.e. areas where the population density is higher than in its surroundings. Today the population density of the world is rather heterogeneous, and the communities range in size from small villages, consisting of a few tens of people, to mega-cities or urban areas where many millions of people are living within a rather small area. A group of such spatially separated populations that are coupled to each other in one way or the other is called a metapopulation (Fig. 1). Within one local community of a metapopulation one may argue that the SIS interactions (1a), (1b) may take place essentially within homogeneous space, so that the epidemic threshold for that community is given by Eq. (7). However, for this equation to be correct one is required to neglect the exchange of people between the various communities of the metapopulation. Today we travel more than ever before; we visit neighboring villages and cities, we travel to further away regional centers, and do international travels to far away destinations for business or pleasure. Hence, the effect of the ever increasing human travel patterns [1517] is that the population centers are coupled, something that in some cases (as we will see below) influences significantly the epidemic threshold.

Fig. 1
figure 1

Schematics of the metapopulation model used in this paper. A real network of subpopulations is mapped onto a complex network. Individuals are distributed over the nodes of the network where the susceptible (light blue men) and infected (dark) individuals interact on the nodes according to the SIS dynamics (Eqs. (3a), (3b)), and travel between the subpopulations (nodes) by diffusion (Color figure online)

In our model, the nodes of the network represent population centers and the links connecting them correspond to the human traveling routes between such centers (Fig. 1). An integrated world will be assumed; by this we mean that any node of the network can be connected to any other node by following the links of the network. Technically this is to say that the network is fully connected and therefore consists of one giant component [29]. We let the network consist of V nodes (or vertices), and let N p (i,t) denote the number of individuals of type p=I,S at population center (or node) i=1,2,…,V at time t. Thus, the population of node i at time t becomes N(i,t)=N S (i,t)+N I (i,t), and by summing over i, the total p-population becomes N p (t)=∑ i N p (i,t). Therefore, the total population of the whole network, independent of the type of individual, can be written as the time-independent constant N=N S (t)+N I (t). Furthermore, the topology of the complex network is described by the so-called adjacency matrix W [29, 30]. The components of this matrix, W ij ≥0, describe the weight or capacity of the traveling route from population center j towards population center i. To keep matters simple, we will assume that the links are symmetric: the capacity of the route ji is equal to that of ij (i.e. W ij =W ji ). Moreover, W ij =0 means that no direct migration (or travel) is possible from j towards i.

To obtain and use a detailed and accurate description of the traveling patterns of man is a complicated and challenging matter [15], and it is outside the scope of this work. Here we will take a simplistic view and use a stochastic migration model where no assumptions are made on how people migrate. Individuals “leaving” population center i will choose a route from center i to j with a probability that is proportional to the capacity W ji . In the next time step, the same procedure is repeated again. This simple migration model is essentially what in physics is called a random walk or a discrete diffusion process on the complex network [3135].

If the susceptible and infected individuals travel according to two independent diffusive processes, it follows from the conservation of “particles” [3234] that the pure migration processes (μ=β=0) satisfy

$$ N_p(i,t+1) = N_p(i,t) + D_p \sum_{j=1}^V T_{ij} N_p(j,t) - D_p N_p(i,t). $$
(8a)

Here D p denotes the “diffusion constant” for p-type individuals, with 0≤D p ≤1, or more precisely, the traveling probability. Moreover, in writing Eq. (8a) we have defined a transfer matrix, T with elements

$$ T_{ij} = \frac{W_{ij}}{\omega_j}, $$
(8b)

where the total outgoing capacity from node j is denoted

$$ \omega_j = \sum _{i=1}^V W_{ij}. $$
(8c)

Under the assumption μ=β=0, the number of S- and I-individuals are separately conserved at any time, and the population sizes in the stationary states of the diffusive migration process are given by

$$ \bar{N}_p(i) \propto\omega_i, \quad p=S,I. $$
(9)

To show that this relation is correct, we start by taking the stationary state limit of Eq. (8a) for which \(N_{p}(i, t+1) = N_{p}(i, t) = \bar{N}_{p}(i)\), then substitute Eq. (9) into the resulting equation, take advantage of Eqs. (8b) and (8c), and finally use that the adjacency matrix is symmetric. The implication of Eq. (9) is that the size of city/village i in the stationary state is proportional to the total migration capacity ω i of that population center, and it is independent of what the diffusion constant might be. As a result of Eq. (9), a separate subpopulation size distribution of villages and cities will not be put explicitly into the model. Instead such information enters implicitly via the distribution of weights ω i , at least it does so for the pure diffusion process. Equation (9) appeals to intuition because larger cities will have transportation systems of larger capacity than smaller cities or villages, in order to serve the larger population. In passing we note that diffusion on complex networks, as described by Eqs. (8a), (8b), (8c), has previously been used to detect network clusters (communities) [30, 3234], and for the cascading failures of complex networks due to overloading [35], to mention a few examples.

4 A Reaction-Diffusion Approach to Epidemic Spreading in a Metapopulation

We will now combine the two phenomena from the previous two sections; infection (reactions) and migration (diffusion). To facilitate future discussion, we introduce the average nodal population density of the network, ρ, defined as the average number of individuals per node (or vertex)

$$ \rho= \frac{N}{V}, $$
(10)

where it is recalled that V is the number of nodes in the network. In terms of population density, Eq. (2) becomes

$$ \rho=\rho_S(t)+\rho_I(t), \qquad\rho_p(t) = \frac{N_p(t)}{V}, \quad p=S,I. $$
(11)

The infection and migration processes of individuals in a metapopulation can now be simulated, given some initial values for N S (i,0) and N I (i,0), a realization of the network, and values for the model parameters. The equations on which such simulations are based are obtained by combining Eqs. (3a), (3b) and (8a), (8b), (8c) so that the relevant set of equations is of reaction-diffusion type. Motivated by the key question in epidemiology studies, we focus on the long time behavior of this reaction-diffusion system. Therefore, our primary interest will be the average nodal density of infected individuals in the stationary state of the system, \(\bar{\rho}_{I}\). The corresponding density of susceptibles is readily obtained as \(\bar{\rho}_{S}=\rho-\bar{\rho}_{I}\) and will thus not be explicitly given. We recall from the preceding discussion that a non-trivial result \(\bar{\rho}_{I}>0\) requires that initially N I (i,0)>0 for at least one node in the network. Therefore, it will from now onward be assumed that N I (i,0)>0 for some i’s. Moreover, if D I is zero, there is no way that an initial infection can be transmitted in our model to other parts of the metapopulation. Since such a “local infection” is not a very interesting and/or realistic situation, we will require that D I >0.

4.1 Numerical Results

Figure 2 shows simulation results (open symbols) for the stationary state of the (scaled) nodal density of infected individuals, \(\bar{\rho}_{I}/\rho\), vs. the average nodal population density ρ/ρ × (with ρ ×=N ×/V) assuming that β/μ=2. The underlying network connecting the subpopulations was assumed to be non-weighted and of scale-free type, p(k)∼k γ, characterized by the exponents γ=3.0 (Fig. 2(a)) and γ=2.5 (Fig. 2(b)). Non-weighted simply means that the population density at given node is proportional its number of incoming links (see Eq. (9)). The simulation results shown in Fig. 2 (as solid lines) were obtained under the assumption of treating ρ p (i,t) as a continuous variable to save computation time, and the networks were generated by the Molloy-Reed algorithm [36]. A stochastic simulation approach similar to that being used in Ref. [10] was also implemented. It was found to give equivalent results except for a higher expenditure of computer time.

Fig. 2
figure 2

The nodal density of infected individuals in the stationary state (\(\bar{\rho}_{I}\)) obtained by simulations (open symbols) or a mean-field approximation (solid lines). The underlying networks were assumed to have a scale-free degree distribution: p(k)∼k γ with γ=3.0 and γ=2.5. The parameters used were D I =1 and β/μ=2, and the various opens symbols labeled in terms of (V, D S ) are: Circles ○ (104, 1), squares □ (102, 1), diamonds ◇ (104, 0) and triangles ▽ (102, 0). We point out that the curves corresponding to the diamonds and triangles (black) fall onto the same line since D S =0. The simulations were based on a (sequential) combination of Eqs. (3a), (3b) and (8a), (8b), (8c) and the analytic results were obtained from Eqs. (22a), (22b) and (24) (Color figure online)

It is clearly seen from Fig. 2 that for sufficiently low nodal population density, ρ, the initial infection will eventually die out and leave the whole population uninfected. Under such conditions, there is no potential risk of the initial infection developing into a large scale epidemic (pandemic). The same simulation also shows that there exists a critical population density, ρ=ρ c , above which an average nodal density will result in a non-zero fraction of the population being infected even in the long time limit. This means that for ρ>ρ c , we are in an epidemic-phase defined by \(\bar{\rho}_{I}>0\). Such a behavior is characteristic for what in physics is known as a phase-transition [37], for instance, observed when water goes from a solid to a liquid phase with increasing temperature.

In passing, it should be mentioned that in the simulations performed to produce Fig. 2, it was assumed that reaction and diffusion were sequential processes, not simultaneous. This is the same assumption made in the original work [10, 11]. However, in a subsequent study [38] it was pointed out that the continuous-time limit of the discrete process used in [10, 11] (and in this work) is not well defined if the reaction-diffusion process were sequential, and the author proposed a formulation for which a continuous-time limit could be obtained. However, doing so renders the formulation more complex with rather similar conclusions. So, we have for the sake of clarity assumed a sequential two-step reaction-diffusion process even if we are well aware of its shortcoming when taking the continuous limit.

4.2 Analytical Treatment of the Model

The critical average nodal density ρ=ρ c is of great importance since it separates the non-epidemic from the epidemic phase, or in other words, defines the model’s epidemic threshold. The purpose of this section is to obtain an analytic expression for this critical population density in terms of the model parameters. To this end, we again follow [10, 11] and impose a mean-field approximation where all nodes of the same degree k are considered equivalent. We introduce the nodal population density

$$ \rho_{p,k}(t) = \frac{N_{p,k}(t)}{V_k}, $$
(12)

where N p,k (t) denotes the total number of p=S,I-type individuals at nodes of degree k in the network, and V k represents the number of nodes of degree k. Just like in the simulations we assume a two-step reaction-diffusion process with the reaction step being the first of the two. Since the reaction step is local to a node, the corresponding equations for the nodal densities ρ p,k (t) are readily obtained from Eqs. (3a), (3b) after a division by V k . The diffusion step, however, needs further discussion because the formalism used in Sect. 3, involving transfer matrices, is no longer optimal. For the migration of individuals of type p into and away from a node of degree k, we obtain the following nodal density

$$ \rho_{p,k}(t+1) = \rho_{p,k}(t) - D_p \rho_{p,k}(t) + k \sum_{k'} p \bigl(k'|k\bigr)\frac{D_p}{k'} \rho_{p,k'}(t), $$
(13)

that follows from the conservation of nodal density of individuals of type p=S or I. The physical interpretation of Eq. (13) is as follows: the first term on the right-hand-side is the density that was already present at nodes of degree k at time t; the second term is associated with the diffusion of a fraction D p of this density away from nodes of degree k; and the last term represents the diffusion into nodes of degree k from neighboring nodes of degrees k′. In writing Eq. (13), we have introduced the conditional probability density, p(k′|k), for a node of degree k being connected to a node of degree k′. If the network is assumed to be uncorrelated, this probability is known to be [39]

$$ p\bigl(k'|k\bigr) = p\bigl(k'\bigr) \frac{k'}{\langle k\rangle }, $$
(14)

where p(k) is the degree distribution and 〈k〉=∑ k kp(k) denotes the average degree. By substituting Eq. (14) into Eq. (13) and performing the sum over k′ in the resulting equation, one is led to

$$ \rho_{p,k}(t+1) = \rho_{p,k}(t) - D_p \rho_{p,k}(t) + \frac{D_p k}{\langle k\rangle } \rho_p(t), $$
(15a)

with

$$ \rho_p(t) = \sum_{k'} p \bigl(k'\bigr) \rho_{p,k'}(t). $$
(15b)

To obtain the final set of coupled equations for ρ S,k (t) and ρ I,k (t) we proceed by first replacing t by t+1 in Eq. (15a) and then replacing the densities ρ p,k (t+1) appearing on the right-hand side of the resulting equation by the right-hand side of Eqs. (3a), (3b) after dividing it by V k . The result is

(16a)
(16b)

where the reaction kernel

$$ \varGamma_k(t) = \frac{ \rho_{S,k}(t) \rho_{I,k}(t) }{ \rho_\times+ \rho_k }, $$
(16c)

and its average

$$ \varGamma(t) = \sum _k p(k)\varGamma_k(t), $$
(16d)

were introduced. In writing Eqs. (16a)–(16d), we have introduced the shorthand notation t ρ S,k =ρ S,k (t+2)−ρ S,k (t) (and similarly for ρ I,k ).Footnote 1

In order to obtain the stationary state of Eqs. (16a) and (16b), we require that \(\rho_{p,k}(t) \equiv\bar {\rho}_{p,k}\) (with p=S,I) is independent of time so that

(17a)
(17b)

where \(\bar{\varGamma}_{k}\) is given by Eq. (16c), but with ρ p,k (t) replaced by \(\bar {\rho}_{p,k}\) and \(\bar{\varGamma }\) is defined as an average over \(\bar{\varGamma}_{k}\) similar to Eq. (16d). Following the strategy from Refs. [10, 11], we then multiply equation (17a) by p(k) and sum over k to obtain

$$ \bar{\rho}_I = \frac{\beta}{\mu} \bar{\varGamma}, $$
(18)

which enables us to simplify Eqs. (17a), (17b) to

(19a)
(19b)

An obvious solution to Eqs. (19a), (19b) is the trivial solution \(\bar{\rho}_{S,k}=\rho\) and \(\bar{\rho}_{I,k}=0\). In the following we will mainly be interested in non-trivial solutions.

Equations (19a), (19b) are valid for a general form of the reaction kernel \(\bar{\varGamma}_{k}\) and have previously been derived in Refs. [10, 11] where it was assumed that the reactions on all nodes of the network were either of “small village type” (ρρ ×) or “big city type” (ρρ ×) (and referred to as reaction type 1 and 2, respectively). However, the results from Refs. [10, 11] do not apply to a mixed reaction type of the form (16c) that, depending on the nodal population density, interpolates between type 1 and 2 reactions. Such dependence on the local density introduces heterogeneous dynamics in the network. In the following, we investigate the effect on the epidemic threshold from this more general reaction type, which we find more adequate for realistic systems.

4.3 Two Special Cases: Migrating and Non-migrating Susceptible Individuals

We will now consider two special cases, namely the situations where the susceptible individuals can or cannot diffuse while the infected individuals always are assumed to diffuse (D I =1). For these two cases we will for simplicity consider D S =0 and D S =1.

Case 1: Non-migrating Susceptible Individuals (D S =0) When the susceptible individuals cannot migrate within the metapopulation, it follows from Eqs. (19a) that

$$ \bar{\rho}_{I,k} = \frac{\beta}{\mu} \bar {\varGamma}_k = \frac{\beta}{\mu} \frac{\bar {\rho}_{S,k} \bar {\rho}_{I,k}}{\rho _\times+ \bar {\rho}_{S,k} + \bar {\rho}_{I,k}}. $$
(20)

In the last transition of Eq. (20) we have taken advantage of the definition (16c) for the reaction kernel. By multiplying Eq. (20) by \((\rho _{\times}+ \bar {\rho}_{S,k} + \bar {\rho}_{I,k})/\bar {\rho}_{I,k}\) and p(k) (the degree distribution), summing over k, and using that \(\rho =\bar{\rho}_{S}+\bar{\rho}_{I}\) (Eq. (11)), leads to

$$ \bar {\rho}_I = \rho\biggl(1- \frac{\mu}{\beta} \biggr) - \frac{\mu}{\beta }\rho_\times, \quad D_S=0. $$
(21)

Formally this expression is only valid when both \(\bar {\rho}_{S}\geq0\) and \(\bar {\rho}_{I}\geq0\). Moreover, since the parameters entering this equation are non-negative, it follows that one must also require β/μ>1. A critical nodal density, ρ=ρ c , is found by solving Eq. (21) under the assumption of \(\bar {\rho}_{I} = 0\). This leads us to write the stationary state of the nodal density of infected individuals in the form

(22a)

where the critical nodal density reads

$$ \rho_c = \frac{\rho_\times}{\frac{\beta}{\mu}-1}, \quad D_S=0, $$
(22b)

which is the same result as for an isolated population (see Eq. (7)). This means that when the susceptible individuals are immobile the epidemic threshold is independent of the network’s topology.

Case 2: Migrating Susceptible Individuals (D S =1) We will now address the situation where both susceptible and infected individuals are allowed to migrate between subpopulations. Under the assumption that D S =D I =1, it follows from Eqs. (19a), (19b) that

$$ \bar {\rho}_{S,k} = \frac{k}{\langle k\rangle }\bar {\rho}_S; \qquad \bar {\rho}_{I,k} = \frac{k}{\langle k\rangle }\bar {\rho}_I, $$
(23)

where the nodal densities of both infected and uninfected individuals are proportional to the node degree. The same result is also found for pure diffusion processes where no infection takes place at all [33, 34]. By introducing relations (23) into Eq. (18), solving the resulting equation for \(\bar {\rho}_{S}\) and using \(\rho= \bar {\rho }_{S}+\bar {\rho}_{I}\), we obtain

$$ \bar {\rho}_I = \rho- \frac{\rho_\times}{\beta/\mu} \biggl[ \sum_k \frac{ p(k) \frac{ k^2 }{ \langle k\rangle ^2 } }{ 1 + \frac{ k }{ \langle k\rangle } \frac{ \rho}{\rho_\times} } \biggr]^{-1}, \quad D_S=1. $$
(24)

This is valid when \(0 < \bar {\rho}_{I} \leq\rho\) and \(\bar {\rho}_{I} =0\) otherwise. Contrary to the D S =0 case, Eq. (24) shows that \(\bar {\rho}_{I}\) is a complicated function depending both on local interaction parameters (β, μ and ρ ×) and global network parameters (p(k)) when susceptible individuals are allowed to migrate. By taking the “small village” and “big city” limits of Eq. (24) one obtains

(25)

where k m denotes the maximum degree of the underlying network (note that β/μ>1). Notice that the quantity \(\frac{k_{m} }{\langle k \rangle } \frac{\rho}{\rho_{\times}}\) is not only large in the “big city” limit where ρ/ρ ×≫1, but potentially also in the “small village” limit (ρ/ρ ×≪1) if the network is sufficiently heterogeneous so k m /〈k〉 is large enough.

The epidemic threshold, ρ=ρ c , is obtained by setting \(\bar {\rho }_{I}=0\) in Eq. (24). In this case, however, we cannot get ρ c in closed form. Rather, it is defined implicitly through

$$ \frac{\beta}{\mu} \frac{\rho_c}{\rho_\times} \biggl[ \sum_k \frac{ p(k) \frac{ k^2 }{ \langle k\rangle ^2 } }{ 1 + \frac{ k }{ \langle k\rangle } \frac{ \rho_c }{\rho_\times} } \biggr] = 1. $$
(26)

Since the summation cannot be evaluated explicitly for a general degree-distribution p(k), one must resort to numerical methods to calculate ρ c for a given set of parameters. Exceptions do however exist. One simple example is a regular grid, i.e. a homogeneous network, where every node has 〈k〉=k links. In this case Eq. (26) yields the same critical density as when the susceptible individuals are immobile (Eq. (22b)), for which there is no dependence on the size of the network (V) or on the number of links each node has (as long as all nodes have the same number of links).

Equations (22a), (22b), and (24) are analytic mean-field expressions for the nodal density of infected individuals in the stationary state. In Fig. 2 these expressions are compared to simulations (discussed in Sect. 4.1) for a few model parameters; They corroborate well with each other.Footnote 2 In particular we point out that the epidemic threshold, the point at which \(\bar {\rho}_{I}\) goes from zero to being non-zero, is well described. Even if the agreement between the simulation and analytic results is satisfactory, it should be stressed that the networks used in the simulations were only characterized by their degree distributions and were otherwise random. Specifically, no clustering or community structures were present [30] and it is left for future research to investigate how well the mean-field results apply to such cases.

4.4 Epidemic Threshold Dependence on Model Parameters

In the previous section we established that the epidemic threshold ρ c is well described by the mean-field results Eqs. (22b) and (26). Here we will discuss in more detail how ρ c depends on the model parameters. In Fig. 3 we present ρ c /ρ × as a function of the rate ratio β/μ for migrating (D S =1) and non-migrating (D S =0) susceptibles for a few different networks. The general trend is that the threshold drops as the infection rate β becomes much larger than the recovery rate μ. Moreover, when D S =1 (or more precisely D S >0), the epidemic threshold also drops when the degree distribution of the network becomes heterogeneous (fat-tailed), or when the size of the network increases.

Fig. 3
figure 3

The epidemic threshold ρ c (scaled by ρ ×), vs. the ratio of the infection and recovery rates, β/μ, for various networks and sizes (V). Equations (22b) or (26) were used to calculate ρ c when D S =0 or D S =1, respectively. The curve labeled by k corresponds to a regular grid where all nodes have degree k. Likewise, the lines labeled by γ correspond to a scale-free degree distribution p(k)∼k γ. For all curves it was assumed that the infected individuals were mobile with D I =1 (Color figure online)

The two overlapping black curves in Fig. 3 correspond to either immobile susceptibles (D S =0) on any network type, or mobile susceptibles (D S =1) on a regular grid where all nodes have k=4 neighbors. The epidemic thresholds for these cases are given by the simple expression ρ c /ρ ×=(β/μ−1)−1 (see Eq. (22b)), and have no dependence on network size.

The two lower curves in Fig. 3 show ρ c for two scale-free networks (p(k)∼k −2.5) of different size. Within a scale-free metapopulation the local population density varies quite dramatically which means that the dynamics in some nodes are of the “small village” type and some (typically much fewer) are of the “big city” type, and the rest somewhere in between. As the network gets bigger, or the “small village” dynamics dominate (\(\rho_{\times}\gg \bar {\rho}_{k}\) for most nodes), the epidemic threshold drops. This behavior can be deduced analytically from Eq. (26) (or Eq. (25)) to giveFootnote 3

$$ \frac{\rho_c}{\rho_\times} \sim \frac{1}{\beta/\mu} \frac{\langle k \rangle ^2}{\langle k^2 \rangle }. $$
(27)

Here we see directly that for a very heterogeneous degree distribution or an infinite network (V→∞), both of which yields 〈k 2〉≫〈k2, the epidemic threshold tends to zero when keeping β/μ finite.

From the general result Eq. (26) we can also deduce the relationship between ρ c and ρ ×. Since this equation is a function of the fraction ρ c /ρ ×, we conclude that they are proportional to each other. The proportionality constant is (β/μ−1)−1 (see Eq. (22b)) for immobile susceptibles (D S =0) and for a regular grid (e.g. k=4 neighbors). For large and heterogeneous networks, or when “small city” dynamics dominate, it is given by Eq. (27).

5 Summary and Conclusions

In this work we have reviewed and applied a generalized reaction-diffusion approach to epidemic spreading within a metapopulation with the overall goal to study the impact of city-size heterogeneity. As the epidemic model we assumed the well-studied susceptible-infected-susceptible (SIS) system. This model exhibits an epidemic and a non-epidemic phase, of which the latter has a nonzero number of infected individuals in the long time limit. The transition between these states, the so-called epidemic threshold, depends non-trivially on the parameters of the model. These parameters are the infection and recovery rates, the cross-over population size separating small village dynamics from large city dynamics, and the topology of the network distributing the subpopulations within the metapopulation. We envision the network to symbolize how cities and villages are interconnected in a transportation network, and human traveling patterns are assumed as simple as possible; People diffuse randomly. In principle, the diffusion rates also enter into the epidemic threshold, but here we only addressed either immobile and mobile susceptibles (D S =0 or D S =1). Instead, we paid special attention to a previously unexplored case, namely, city-size dependence on epidemic spreading. By city size we mean population density. It is an interesting aspect because a disease spreads quite differently in a small city (a village) compared to in a metropolitan city, and the most realistic situation is that all city sizes coexist within a metapopulation. We therefore formally let the local dynamics in each node adapt itself to the instantaneous local population density. In a diffusion model, the size of the city is proportional to the number of its links (we assume an undirected, non-weighted network). This means that the network’s structural features enter the dynamics of the disease in two ways. First, how people migrate, and second, the dynamics in each node. Previous works have assumed that the individuals are infected in the same way in all nodes. Our main result is that we obtain (in the mean-field limit) for a general network an implicit expression for the epidemic threshold expressed in terms of the model parameters. This expression shows good agreement with numerical simulations and its asymptotic behavior interpolates between known limits.

Epidemic spreading is a topic of great social importance that also has attracted the interest of diverse scientific disciplines for decades. The results presented in this paper open for a better understanding of the position of the epidemic threshold and its complex interplay between the local infection parameters and the global parameters characterizing the network topology for a large class of diseases (SIS). Such knowledge may in the future prove critical for reducing the outbreak size, slowing down the speed of spreading and/or spatially confining infectious diseases. Moreover, if epidemic dynamics were better understood, our contemporary societies could earlier, and at better levels of accuracy, predict if a local outbreak has a risk of becoming a global pandemic. If an accurate early warning could be given, appropriate measures could be taken, for example by adjusting the production of vaccines, with the potential of saving numerous lives.