A sharp threshold for disease persistence in host metapopulations

Abstract

A sharp threshold is established that separates disease persistence from the extinction of small disease outbreaks in an S→E→I→R→S type metapopulation model. The travel rates between patches depend on disease prevalence. The threshold is formulated in terms of a basic replacement ratio (disease reproduction number), ℛ₀, and, equivalently, in terms of the spectral bound of a transmission and travel matrix. Since frequency-dependent (standard) incidence is assumed, the threshold results do not require knowledge of a disease-free equilibrium. As a trade-off, for ℛ₀>1, only uniform weak disease persistence is shown in general, while uniform strong persistence is proved for the special case of constant recruitment of susceptibles into the patch populations. For ℛ₀<1, Lyapunov's direct stability method shows that small disease outbreaks do not spread much and eventually die out.

Keywords:

• Persistence
• Extinction
• Basic reproduction number

1. Introduction

Many compartmental epidemic models focus on how infectious diseases develop with time and ignore the other category that shapes reality, space. Spatial disease spread can be included with distance as a continuous variable leading to partial differential equation models as in 1, 7, 17, 18, to integro-differential equations as in 25, 26, or to equations of mixed type 21. (For the spatial spread of epidemic outbreaks see 20, 21, 31 and the references therein.) Alternatively, space can be included as a discrete variable leading to (spatially explicit) metapopulation models as in 24 that consist of a (possibly very large) system of ordinary differential equations (ODEs). The underlying idea is to subdivide the spatial region under consideration into a number of discrete patches, which may represent districts, cities, countries, ponds (for the study of amphibian decline), etc. The disease is carried from one patch to another by individuals traveling between patches. As a simplifying assumption, individuals only change disease status when they are on a patch and not when they travel. (See 5 for a model with disease transmission during travel.) Work on S→I→S type models 2, 23, 33–36 and on models for influenza 13, 14 indicates that a patchy environment and travel between patches can influence disease spread in a complicated way and, depending on parameters values that include characteristics of the specific disease, can enhance or stifle disease spread.

Our aim is to add to the threshold analysis of the type metapopulation model in 23 by discussing dynamic disease persistence and the existence of an endemic equilibrium. We also want to extend the threshold analysis to travel rates that depend on the disease prevalence on the patches and to recruitment rates that depend in addition on the patch population densities (see section 2). The variable travel rates take into account behavioral changes that may occur as individuals adjust to such factors as the severity of the disease or travel restrictions. We prove the existence of a sharp threshold that separates disease persistence (sections 4 and 5) from extinction of small disease outbreaks (section 6). This threshold can be formulated in terms of a basic replacement ratio (basic parasite reproduction number), ℛ₀, or, equivalently, in terms of the spectral bound of a transmission and travel matrix (section 3). By contrast with 23, our results (except for Theorem 6.3) do not depend on the existence of a disease-free equilibrium because our methods fully exploit the crucial choice that disease transmission is modeled by frequency-dependent (standard) incidence. The threshold condition for density-dependent (mass action) incidence would involve the patch population densities at a unique disease-free equilibrium (see the Discussion). As a trade-off, our approach struggles with the fact that, under frequency-dependent incidence, the disease can drive both host and parasite to extinction if it induces fatalities or reduces fertility 9, 29, 38. This is why, under general recruitment, we can only establish uniform weak persistence of the disease if (section 4). For the same reason, we formulate disease persistence in terms of the frequencies rather than the densities of the infectives on the patches. So disease persistence does not necessarily imply persistence of the causative agent (parasite, pathogen).

For uniform strong disease persistence, we assume constant recruitment into the patch populations (section 5) as in 23. Under this strong assumption, not only the disease but also the host and the parasite persist, and there exists an endemic equilibrium (a question left unanswered in 23). Alternatively, the assumption could be made that the disease neither induces fatalities nor reduces fertility as has been done in 34, 35 for an S→I→S type model.

The threshold condition is sharp as we show that small disease outbreaks do not spread much and eventually die out if (section 6). Since the disease-free dynamics are unclear, we do not use a linearized stability analysis, but Lyapunov's direct stability method.

The epidemic metapopulation model introduced in 24 not only keeps track of the patch where an individual is currently located (as in the model analyzed in this paper) but also of the patch on which an individual is born and usually resides. For m epidemiological classes on n patches such models lead to a system of mn ² rather than mn ODEs. Under the assumption that travel is independent of disease prevalence, they are analyzed for the S→I→S case in 3, for the case in 4, and for a model of SARS that includes quarantine in 22. An explicit expression for the basic reproduction number is given, and numerical simulations for the S→I→S model indicate that this number acts as a threshold between extinction and persistence of the disease. It remains to be explored whether the methods we develop in this paper will also work for these even more complex models with frequency-dependent incidence.

2. The model

We consider a host metapopulation that is geographically distributed over n patches (districts, countries). The disease divides each patch population into four classes (or compartments): susceptible individuals (represented by letter S), exposed individuals (infected but not yet infective, represented by letter E), infective individuals (represented by letter I), and removed (or recovered) individuals (represented by letter R). We include the possibility that recovered individuals can become susceptible again, so our model will be of type. We assume that disease transmission and transition between disease classes only occurs on patches and not during travel. See 5 for a model with disease transmission during travel.

Let denote the respective numbers of susceptible, exposed, infective, and removed individuals and N _i (t) denote the total population size in patch i at time t (see table 1). We introduce the following notation for the various epidemiological classes,

The letter

is used for an epidemiological class. The 3n-tuple J represents the state of the non-susceptible part of the host population. The dynamics of the host population and the disease are described by the following system, for i=1, …, n,

Individuals are recruited into the local population of patch i at a rate . This may happen by immigration or birth, or, in the case of a sexually transmitted disease, by entering the sexually active part of the patch population. Individuals on patch i and in the class die at a per capita rate and move to the next epidemiological class at a per capita rate . Disease transmission is modeled by frequency-dependent (standard) incidence, with κ _i being the per capita infection rate on patch i (consult the Discussion for the consequences of this choice). Individuals in epidemiological class C travel from patch k to patch i at a per capita rate . Without loss of generality we can assume that

where . More generally than in 23, the recruitment rate may depend on the state of the population for some of our results. Furthermore individuals may adjust their travel rates to the disease status of the population (table 2).

Table 1. Model variables.

t	time
i	patch number
S i (t)	number of susceptibles in patch i at time t
E i (t)	number of exposed individuals in patch i at time t
I i (t)	number of infectives in patch i at time t
R i (t)	number of removed individuals in patch i at time t
N i (t)	population size of patch i at time t

Table 2. Parameters and parameter functions.

	per capita mortality rate of class C in patch i
	per capita transition rate in patch i from class C to the next
κ i	per capita infection rate on patch i
	recruitment rate of patch i
	per capita travel rate of class C from patch k to patch i

The following assumptions are made throughout the paper.

Assumption 2.1

The functions are defined for , . The functions , , are defined for . These functions are nonnegative and are locally Lipschitz continuous on their respective domains. All parameters are positive with the possible exception of some or all of which may be 0. Furthermore,

• (a) ,

• (b) ,

• (c) the matrix is irreducible,

• (d) the matrix is irreducible.

All norms are sum-norms, . So ‖J‖ is the total number of non-susceptibles and ‖(S, J)‖ is the size of the total population. S≫0 means that , i.e. the vector S has all its coordinates positive.

Assumption 2.1 (a) is satisfied, e.g. if all recruitment rates are bounded. Assumption 2.1 (b) guarantees that each local population survives in the absence of the disease. By Assumption 2.1 (c) and (d), every patch can be reached from every other patch by susceptible and infected individuals.

For N _i =0, we define . The expression is then a locally Lipschitz continuous function of . We define the total population size as

Notice that if all vectors are positive (see table 3).

Table 3. Vector and other notation.

	J=(E, I, R)
	generic epidemiological class
	total size of epidemiological class C
	total population size

Theorem 2.2

For all , there exists a unique solution , R(t)) of (2) with initial data that is defined for all t≥0 and takes values in . Further there exists some c>0 such that, for all non-negative solutions of (2), and for all t≥0.

Proof

By Theorem [30; Theorem A.4], there exists a unique solution with values in that is defined on some interval [0, b) with . If b<∞, then . We add all equations in system (2),

Set , then μ>0 and

By Assumption 2.1 (a), we can choose c>0 such that

By (6), whenever t∈[0, b) and . This implies for all t∈[0, b). So b=∞, and .   ▪

3. The basic replacement ratio, ℛ₀

To establish a threshold parameter that separates disease persistence from the extinction of small disease outbreaks, the following linear operator B on plays a crucial role: for i=1, …, n,

where . The operator B can be represented by a quasi-positive matrix which is irreducible by Assumption 2.1 (d) and can be rewritten as

with appropriate linear operators on ℝ ⁿ . The positive operators D ^I and are represented by diagonal matrices with all diagonal entries being positive. The diagonal entries of D ^I are the per capita infection rates κ _i , while the diagonal entries of are the per capita transition rates from the exposed to the infectious disease stage. The operator P is represented by a matrix with off-diagonal entries and diagonal entries . For Q, the entries are the same as those for P except that the superscript E is replaced by the superscript I.

Let s(B) denote the spectral bound of a linear (bounded) operator B, i.e. the largest real part of its eigenvalues, and let r(B) denote its spectral radius, i.e. the largest absolute value of its eigenvalues. The operators P and Q have the form D−M, where D is represented by a diagonal matrix with all diagonal entries being positive and M is represented by with for P and for Q. Thus −P is represented by a quasi-positive matrix. Further, the column sums of −P are negative, which implies that s(P)>0 and r(D ⁻¹ M)<1 (6, Theorem 6.13], [30, Theorem A.44], 32). With 𝕀 denoting the identity operator,

and P ⁻¹ (and also Q ⁻¹) is a positive operator represented by a matrix with strictly positive diagonal.

Proposition 3.1

With B given by (8), s(B) and have the same sign.

Proof

Let s(B)>0. By the Perron–Frobenius theory (see, for example 27; A or 30, section A.8, there exists an eigenvector (v, w)≫(0, 0) of B with B (v, w)≫(0, 0). By (8),

Since P ⁻¹ and Q ⁻¹ are positive operators represented by matrices with strictly positive diagonals,

We substitute one inequality into the other, . We set . Then . This implies that the spectral radius r of strictly exceeds 1. Now let s(B)=0. By the Perron–Frobenius theory, there exists an eigenvector (v, w)≫(0, 0) of B with B (v, w)=(0, 0). By (8),

Proceeding similarly as before, . We set . Then . Since z≫0, the spectral radius r of is one. Similarly we show that s(B)<0 implies r<1.   ▪

As in 23,

can be interpreted as the basic reproduction number (basic replacement ratio) of the disease 6, 32. ℛ₀ will turn out to be the threshold parameter that separates disease persistence (if from the extinction of small outbreaks of the disease (if ). An alternative threshold parameter is s(B), the spectral bound of the travel and transmission operator B given by (7).

4. Uniform weak disease persistence if ℛ₀>1

Since our choice of frequency-dependent incidence makes it possible that the disease drives the host (and the causative disease agent) to extinction 9, 29, 38, we formulate disease persistence in terms of the frequencies rather than the densities of infective individuals on the patches.

Definition 4.1

The disease is said to be uniformly weakly persistent if there exists some ε>0 such that

for all nonnegative solutions of (2) with and S(0)≫0.

We will show that the disease is uniformly weakly persistent if . This condition is sharp, as we will show in section 6 that, if , small disease outbreaks die out. Uniform weak disease persistence is a stronger concept than instability of the disease-free state which is proved in 23 for . In our model, which is more general than the one in 23, there is even no clear candidate for the disease-free state.

We will use the following Lemma several times.

Lemma 4.2

There exist such that the following holds for all nonnegative solutions of (2) with ‖S(0)‖>0: If then i=1, …, n.

Proof

We obtain the following differential inequality from the first subsystem in (2),

For δ₁>0 small enough, Assumption 2.1 (c) provides a positive operator A represented by an irreducible n×n matrix (α _ik ) such that

By Theorem 2.2, we can choose c>0 and t ₂≥1 large enough such that satisfies the differential vector inequality

Here c>0 can be chosen independently of the solution while t ₂ may depend on the solution. We integrate,

Since A is an irreducible nonnegative matrix, all entries of e ^A are positive. Let ξ be their minimum. Then ξ>0 and for all t≥t ₂. Set .   ▪

Theorem 4.3

If the disease is uniformly weakly persistent.

Proof

Suppose the statement is false. Choose an arbitrarily small ε>0. By Definition 4.1, there exists a solution of (2) such that and S(0)≫0, but

Then S(t)≫0 for all t≥0 and and for all t>0. Shifting forward in time we can assume that

By Theorem 2.2, there exists some c>0 such that . Shifting forward in time again and increasing c, we can assume that for all t≥0. For , we define

Then

We complete the proof in three steps.

Step 1 There exists some such that

To prove step 1, we obtain the following differential inequality from the second subsystem in (2),

We add over i and set , ,

This implies

For R¯, the proof is similar.

Step 2 There exists some δ₂>0 such that , i=1, …, n.

To prove step 2, we derive the following differential inequality from the first subsystem in (2) and (11), with ,

By Assumption 2.1 (b), there exist some such that

Let . We choose ε>0 so small that, by Step 1, there exists some t ₁>0 such that for all t≥t ₁. So for all t≥t ₁. Choose ε>0 so small that . By (12),

This implies that

Step 2 now follows from Lemma 4.2.

Step 3 The contradiction.

Combining step 1 and step 2 and shifting forward in time, we can assume that , , and for i=1, …, n and C=E, I, R. The positive constants and do not depend on ε. Then and, from the second and third subsystem in (2),

By step 2 and by Theorem 2.2, again shifting forward in time if necessary, we can assume that for all t≥0, i=1, …, n with δ>0 and c>0 not depending on ε. Since the functions are continuous, as ‖J‖→0. So, for any η∈(0, 1) we can choose ε>0 small enough such that, by step 1,

By (14), we have a differential vector inequality

with a linear operator B _ε that is associated with an irreducible quasi-positive matrix by Assumption 2.1 (d). Moreover as , where B is the operator in (8). Since , s(B)>0 by Proposition 3.1. The eigenvalues depend continuously on the operator (or representing matrix) and the spectral bound is an eigenvalue, therefore for sufficiently small ε>0. By the Perron–Frobenius theory (see, for example, 27; A or 30; section A.8), we can choose a vector v≫0 such that where * denotes the dual operator. Then

Since v≫0, first for t=0 and then for all t≥0. Since , as t→∞. Without loss of generality we can choose v such that v _j ≤1 for j=1, …, 2n, which implies as t→∞. This contradiction to Theorem 2.2 completes the proof.   ▪

5. Uniform strong disease persistence and existence of endemic equilibria

Typically, uniform weak implies uniform strong persistence (and the existence of an endemic equilibrium) if the dynamical system has a compact attractor. Unfortunately, Theorem 4.3 has the proviso that S(0)≫0 (see Definition 4.1), which must be built into the state space. So the existence of a compact attractor would imply strong host persistence. However, since our model assumes standard (frequency-dependent) incidence, it is possible that the disease drives the host to extinction 9, 29, 38. A suitable assumption that rules this out can be found in 23, namely that all recruitment rates Λ _i are positive constants. We adopt it for this section. Under the assumption of constant recruitment, it is easy to see that the host is uniformly strongly persistent.

Lemma 5.1

If all Λ _i are positive constants, then, for i=1, …, n, S _i (t)>0 for all t>0, and there exist constants δ _i >0 such that for all nonnegative solutions of (2).

Proof

By the first subsystem in (2),

By the fluctuation method (12, [30; Proposition A.22]),

By Theorem 2.2, there exist , independent of the solution, such that

   ▪

If , the parasite is uniformly strongly persistent as well, at least in total abundance. Since S(t)≫0 for t>0 by Lemma 5.1, the subsequent persistence results do not need the solutions of (2) to satisfy S(0)≫0.

Theorem 5.2

If and all Λ _i are positive constants, there exists some ε>0 such that

for all nonnegative solutions of (2) with .

Proof

We apply [30; Theorem A.32]. Let

By Lemma 5.1, the solution takes its values in X for t>0. Define by

In the language of 30; section A.5, the semiflow Φ induced by the solutions of (2) is uniformly weakly ρ-persistent by Theorem 4.3. The compactness condition in 30, section A.5 follows from Theorem 2.2 and Lemma 5.1. By [30; Theorem A.32], Φ is uniformly strongly ρ-persistent. The claim of this theorem for C=I now follows from Lemma 5.1 and the fact that implies for all t>0. Once we have the statement for C=I, it follows for C=E by Lemma 5.1 and also for C=R.   ▪

Corollary 5.3

If and all Λ _i are positive constants, there exists an equilibrium state of (2) with for C=E, I, R.

Proof

We apply [37; Theorem 1.3.7]. Define X as in the proof of Theorem 5.2 and . Then X ₀ is convex and relatively open in X. By Theorems 2.2 and 4.3 and Lemma 5.1, the semiflow induced by the solutions of (2) is uniformly persistent with respect to in the language of 37. By [37; Theorem 1.3.7], (2) has an equilibrium in X ₀. It easily follows that not only for C=I, but also for C=E, R.   ▪

Under a further irreducibility assumption concerning the migration of exposed or infective individuals, the parasite persists uniformly strongly on each patch.

Theorem 5.4

Let and all Λ _i be positive constants. Further assume that for each the matrix is irreducible. Then there exists some ε>0 such that

for all solutions with Further there exists an (endemic) equilibrium of (2) in (0, ∞)⁴ⁿ .

Proof

Let the state space X be as in Theorem 5.2. Fix and define by . By Theorem 2.2 and Lemma 5.1, the compactness condition of [30;Theorem A.34] is satisfied. Notice that every total orbit of Φ is associated with a solution of (2) that is defined for all times and takes value in X. By our irreducibility assumption, whenever for all . The claim for C=I now follows from [30; Theorem A.34]. For C∈{E, R}, modify . For C=S, the statement has already been shown in Lemma 5.1. The existence of an equilibrium in (0, ∞)⁴ⁿ follows from [37; Theorem 1.3.7].   ▪

6. Extinction of small outbreaks if ℛ₀<1

In order to illustrate that the condition is sharp for disease persistence, we show that small disease outbreaks do not spread much and eventually die out if . Recall that J=(E, I, R) and all norms are sum norms. So ‖J‖ is the total number of individuals that are not susceptible.

Theorem 6.1

If , the following local stability results hold for nonnegative solutions of (2) .

• For every ε>0 there exists some δ>0 such that for all t≥0 provided .

• Further there exists some δ₀>0 such that as t→∞ whenever .

Proof

We use Lyapunov's direct stability method. The second and third subsystems of (2) imply the inequalities

which can be written as

with the linear operator B _J being given by the right-hand sides of the system (15). By Assumption 2.1 (d), B _J is represented by a quasi-positive irreducible matrix for J=0. Since , s(B _J )<0 for J=0 by Proposition 3.1. Since the eigenvalues of a matrix depend continuously on the entries of the matrix and s(B _J ) is an eigenvalue, we find some δ₁>0 and a quasi-positive irreducible matrix such that and whenever and . By the Perron–Frobenius theory, there exists an eigenvector such that . The eigenvector v serves to construct the local Lyapunov function . Let ε>0. Choose Let and

We claim that τ=∞ if δ is chosen small enough. Suppose τ<∞. Then for all t∈[0, τ] and . So

and

So is a decreasing function of t∈[0, τ] and in particular I(t)), . Since v≫0, there exists some c≥1, only dependent on v, such that

Let and . Then, for t∈[0, τ],

We integrate this inequality,

with . So, for t∈[0, τ],

We choose δ so small that . But then we find some such that for all , a contradiction to (17). So τ=∞ and for all t≥0. To obtain the first part of our stability statement, we choose δ>0 small enough that . For the second statement we choose δ₀>0 such that . Let . Then for all t≥0. By (18), which now holds for all t≥0, as t→∞. By (19) and the fluctuation method,

The preceding result does not preclude that both parasite and host go extinct, which is a possibility as our model uses standard (frequency-dependent) incidence 9, 29, 38. The next result shows that, if , the disease dies out and the host persists if the initial numbers of exposed, infective or removed individuals are small in comparison to the initial number of susceptible individuals.

Theorem 6.2

If , there exist such that the following holds for all nonnegative solutions of (2) : for any there exists some δ>0 such that

Proof

Let η>0 to be chosen later. By Theorem 6.1, there exists some δ>0 such that for all t≥0 whenever . Let . We add the first subsystem in (2),

Then

with . By Assumption 2.1 (b), there exist such that

We choose . Then

Let and . Let t be a time at which if such a time exists. Then

We now choose η small enough (in dependence of ε) so that the expression in the parentheses above is positive. Then whenever . So for all t≥0 and, by Lemma 4.2, >0, i=1, …, n. By Theorem 6.1,

By (20) and (21),

If η>0 is small enough (recall for all t≥0), we obtain from Lemma 4.2 that , i=1, …, n, for some that does not depend on the solution.   ▪

We now assume as in section 5 that the recruitment rates Λ _i are positive constants for i=1, …, n. Then the disease-free system, and

has a globally stable equilibrium . Indeed, (22) can be rewritten as with the operator A having a spectral bound s(A)<0. By the variation of constants formula,

Theorem 6.3

If and all Λ _i are positive constants, there exists some with the following property: for any there exists some δ>0 such that

Proof

By Theorems 6.1 and 6.2, the first subsystem (the one for S) of system (2) is asymptotically autonomous with limit system (22). The claim follows from the global stability of S ^⋄ for (22) and results on asymptotically autonomous systems in 15 (or 28).   ▪

If the travel matrices are independent of J, the disease-free equilibrium is globally asymptotically stable if [23; Theorem 2.2].

7. Discussion

The sharp threshold results in this paper that separate disease persistence (if ) from the extinction of small disease outbreaks (if ) crucially depend on modeling disease transmission by frequency-dependent (standard) incidence. If frequency-dependent incidence were replaced by density-dependent (mass action) incidence, ℛ₀ would involve the population densities on all patches at a disease-free equilibrium, which would need to be unique and globally attracting for the disease-free population dynamics. As a trade-off, for frequency-dependent incidence, it is possible that the disease drives the host and the causative disease agent to extinction 9, 29, 38 while host persistence is typically automatic for density-dependent incidence. For a complex model like this, it seems very hard to establish a threshold for host persistence unless one makes assumptions such as the one made in section 5 and 23, constant recruitment into each patch population, or the one made in 34, 35, no disease fatalities and no fertility reduction. Strangely enough, this problem also has repercussions for disease persistence: we can only prove uniform weak disease persistence in general, and host persistence is required in order to go from uniform weak to uniform strong disease persistence under the assumptions that we have imposed.

We refer to 10, 16 and the references therein for the discussion as to whether and when frequency-dependent or density-dependent incidence is more adequate. The choice of one over the other not only has consequences for whether or not the disease can eradicate the host, but also for the occurrences of undamped oscillations: they seem to occur more often for density-dependent incidence than for frequency-dependent incidence 8, 11, 19.

Acknowledgements

The authors thank two anonymous referees for helpful comments. This study was partially supported by NSF grant DMS 0314529 (T.D. and H.R.T.) and by NSERC and MITACS (P.v.d.D).