Persistence time of SIS infections in heterogeneous populations and networks

# Persistence time of SIS infections in heterogeneous populations and networks

Damian Clancy
Department of Actuarial Mathematics and Statistics
Maxwell Institute for Mathematical Sciences
Heriot-Watt University
Edinburgh
EH14 4AS
UK
d.clancy@hw.ac.uk
###### Abstract

For a susceptible-infectious-susceptible (SIS) infection model in a heterogeneous population, we present simple formulae giving the leading-order asymptotic (large population) behaviour of the mean persistence time, from an endemic state to extinction of infection. Our model may be interpreted as describing an infection spreading through either (i) a population with heterogeneity in individuals’ susceptibility and/or infectiousness; or (ii) a heterogeneous directed network. Using our asymptotic formulae, we show that such heterogeneity can only reduce (to leading order) the mean persistence time compared to a corresponding homogeneous population, and that the greater the degree of heterogeneity, the more quickly infection will die out.

Keywords: Stochastic epidemic models; large deviations; endemic fade-out; directed configuration model; superspreaders

## 1 Introduction

In modelling endemic infections, a quantity of particular interest is the persistence time until infection dies out from the population. For discrete state-space Markov chain models, the expected persistence time for an infection that has become endemic in the population (i.e. starting from quasi-stationarity) may be found as an eigenvalue of the transition rate matrix. However, for large populations and for more complicated models, numerical computation of this exact solution can be very time-consuming, and may also suffer from numerical instability. Moreover, it is not straightforward to use this eigenvalue characterization to investigate, for instance, the effect of population heterogeneities upon the expected persistence time. Approximation methods are therefore required. For a number of infection models, it has been shown that, denoting by the typical size of the population, the expected time from endemicity to extinction, , is asymptotically given by an expression of the form

 τ ∼ C√Nexp(AN) (1)

where the values of depend upon parameters of the model, but not upon . We assume here that the process is super-critical, so that long-term endemicity is possible.

For the classic susceptible-infectious-susceptible (SIS) model of Weiss and Dishon [27], Andersson and Djehiche [1] found simple explicit expressions for both and in terms of the basic reproduction number (the expected number of secondary cases caused by a typical primary case in an otherwise susceptible population), under the assumption of super-criticality (that is, ); specifically, and , assuming time is scaled such that individual infectious periods are of mean 1. This was extended by Ball, Britton and Neal [4] to allow for a general infectious period distribution in place of the exponential distribution assumed by [1]; they showed that leading-order behaviour is unchanged, so that as before, while the value of depends upon the infectious period distribution, and may be straightforwardly evaluated provided this distribution is known. Pre-dating the above work, van Herwaarden and Grasman [26] showed that relationship (1) holds true for a particular suceptible-infectious-removed (SIR) infection model. In this case, however, evaluation of the constant  requires numerical solution of a system of ordinary differential equations, while no method for evaluating  is given.

The system of ordinary differential equations used in [26] to evaluate  may be regarded as the equations of motion corresponding to a particular Hamiltonian system. More recently, a number of authors [2, 3, 13, 14, 19, 21] have applied this Hamiltonian approach to a range of infection models to derive results of the form

 limN→∞lnτN = A. (2)

Equation (2) is not as precise as relationship (1), but does at least give the leading-order behaviour of in the large population () limit. Evaluation of  generally requires numerical solution of the equations of motion, and consequently much of the research effort has focused upon developing efficient numerical procedures.

We shall apply the Hamiltonian approach to approximate the expected persistence time from endemicity, , for an SIS model incorporating heterogeneity in individual infectivities and susceptibilities. Such heterogeneity is a common feature of real-world infections. For instance, for a number of infections (eg SARS) it has been hypothesised that there exists a subgroup of ‘super-spreaders’ within the population, being individuals of higher infectivity than the rest. Heterogeneous susceptibilities may arise, for instance, through individuals having differing histories of prior exposure to infection or vaccination. Alternatively, our model may be interpreted as a model for infection spreading on an uncorrelated (that is, with no correlations between degrees of neighbouring individuals) directed network [12].

In contrast to almost all previous work, we are able to find an explicit formula for the constant  in equation (2), at least provided the heterogeneity is in either infectivity, or susceptibility, but not both. As well as being much quicker and easier to evaluate than the solution to a (typically high-dimensional) system of ordinary differential equations, a further advantage of such an explicit formula is that it may be used to establish qualitative results about the effects of model assumptions. Specifically, we investigate the effect of increasing heterogeneity upon the persistence time of infection in the population.

The remainder of the paper is structured as follows. In section 2, we define precisely our heterogeneous population SIS model, and describe how it may be interpreted as approximating a directed network model. In section 3 we recall some general theory that will be required in the sequel. Our main result, theorem 1, is derived in section 4, establishing explicit asymptotic formulae for in the large-population limit, provided that heterogeneity is in either infectivity or susceptibility, but not both. Using these explicit formulae, we go on in section 5 to demonstrate that the greater the level of heterogeneity in either infectivity or susceptibility, the more rapidly extinction of infection will occur (on average, to leading order in a large population). In the case that both heterogeneities are present simultaneously, numerical work (figure 2) suggests that mean persistence time is again maximised in the homogeneous-population case. In section 6 we demonstrate that if heterogeneity is in susceptibilities, our asymptotic formula for , and hence also our conclusion that greater heterogeneity reduces mean persistence time, remain valid for more general infectious period distributions than the classic exponential distribution. We present numerical evidence (figure 3) suggesting that this is also true when instead heterogeneity is in infectivities. Finally, in section 7, we discuss the directed network interpretation of our results (theorem 4) and suggest some directions for further work.

## 2 The SIS infection model in a heterogeneous population or directed network

We first formulate our model in terms of a population divided into a fixed number of groups, and then describe how the same model may be interpreted as modelling an infection spreading on a directed network.

Consider a closed population of individuals divided into groups, with group  () consisting of individuals. Denote by the proportion of the population belonging to group , so that . When a group  individual becomes infected, it remains so for a time distributed as an exponential random variable with mean (assumed for simplicity to be the same for each group). During this infectious period, the group  infective makes contacts with each individual in each group  at the points of a Poisson process of rate , where is some overall measure of infectiousness, represents the infectivity of group  individuals, and represents the susceptibility of group  individuals. (The assumption that the group  to group  infection rate factorises in this way is sometimes referred to as ‘separable mixing’.) Without loss of generality, we scale the , values so that . These Poisson processes and infectious periods are all mutually independent. If a contacted individual is susceptible, then it becomes infected (and infectious); if the contacted individual is already infected then the contact has no effect. Thus the process is a continuous-time Markov chain with transition rates given in table 1, and the number of susceptible individuals in group  at time  is . We will assume throughout that , and that for all . Note that our model is a special case of the model of [8], although we use slightly different notation here. The basic reproduction number is given by the dominant eigenvalue of the matrix with entries , so that

 R0 = βγk∑i=1λiμifi.

We now describe how the above model may be interpreted as describing infection spreading through a network. Each of the individuals in the population is assigned an in-degree  and out-degree  according to some joint probability mass function on . These degrees are assigned independently to distinct individuals, but the in and out degree need not be independent for a single individual. To each individual we attach ‘stubs’, or half-edges, with stubs pointing inwards and stubs pointing outwards. Inward-pointing stubs are then paired with outward-pointing stubs throughout the population, to create links between individuals. In order that this process can produce a valid network, with no left over half-edges, we clearly require that . We do not concern ourselves with the precise mechanism by which stubs are paired off (see [5, 6] for relevant discussion); rather, we shall simply assume that the resulting network is uncorrelated, so that the so-called ‘annealed’ network approximation is valid for an ensemble of such networks. This is a mean-field approximation for heterogeneous networks, and may be described as follows (see [12] for more comprehensive discussion).

Denote by the rate at which infection transmits along each link from an infectious individual to a susceptible individual. Suppose for simplicity that there are a finite number  of pairs having non-zero probability, and define a bijective function that assigns a unique number to each of the possible pairs. We say an individual is of ‘group ’ if they have degrees , and define to be the in and out degrees, respectively, of a group  individual. For , denote by the total number of group  individuals in the population, and by the number of group  individuals who are infectious at time . Then under the annealed network approximation, the total rate at which group  individuals become infected is given by

 κNE[din](k∑m=1dout(m)Im)din(j)(Nj−Ij).

When individuals become infected, they remain so for an exponentially distributed time of mean  before returning to the susceptible state.

This network model may be approximated by the -group model with transition rates given in table 1 by taking

 β = κE[dout], fj = p(c−1(j)), μj = din(j)/E[din], λj = dout(j)/E[dout],

for .

The undirected version of the above annealed network approximation (with ) has been studied by [17], via numerical solution of Hamilton’s equations of motion, equations (15,16) below.

In the next section we present some relevant general theory, before going on in section 4 to apply these general methods to the model described above.

## 3 General theory regarding persistence time from endemicity

Consider an infection modelled by a continuous-time Markov process on finite state-space  with transition rate matrix . Suppose that  is made up of an absorbing set of states  (corresponding to absence of disease) and a single transient communicating class . We denote by the transition rate matrix restricted to . Then the infection will almost surely die out (i.e. the process will leave ) within finite time, and [10] there exists a unique quasi-stationary distribution such that, for any initial state within ,

 qx = limt→∞Pr(X(t)=x|X(t)∈C) for x∈C.

That is, provided the infection does not die out, it will settle to the endemic distribution . The distribution may be found as the unique solution of

 qQC = −(1/τ)q with ∑x∈Cqx=1, (3)

where is the eigenvalue of with largest real part. The time to extinction from quasi-stationarity is exponentially distributed with mean .

Although may be computed exactly from equation (3), this can become impractical when the state-space is large, and it is not straightforward from (3) to establish qualitative results. Approximation methods are therefore valuable, and in particular, methods from Hamiltonian statistical mechanics may be used to study the leading order asymptotic (large population) behaviour of , as follows.

Suppose that is a density-dependent process in the sense of chapter 11 of [15]; that is, the transition rates are of the form

 P(X(t+δt)=x+l∣X(t)=x) = NWl(xN)+o(δt) for x∈S, l∈L, (4)

for some functions , where is the set of possible jumps from each state and is some parameter indicating overall size of the system (in our applications, will be the size of the population). Under mild technical conditions ([15], Theorem 11.2.1), the scaled process converges almost surely over finite time intervals, as , to the solution of the ordinary differential equation system

 dydt = ∑l∈LlWl(y). (5)

For our application, we suppose that the system (5) possesses two equilibrium points: a stable endemic equilbrium point with all components strictly positive, and an unstable disease-free equilibrium point at . We next summarise some key results from the Hamiltonian approach, in a form suited to our application. Detailed justifications and extensions of the method may be found in the review paper [3] and references therein.

The Hamiltonian of the system is defined to be

 H(y,θ) = ∑l∈LWl(y)(eθTl−1). (6)

This Hamiltonian determines the following two complementary Hamilton-Jacobi partial differential equations:

 H(y,∂V∂y)=0 and H(∂U∂θ,θ)=0. (7)

Each of these Hamilton-Jacobi equations is a way of expressing the eigenvector equation (3) while retaining only leading order terms in the limit (see Appendix for a brief outline of the derivations).

If we can solve either of the Hamilton-Jacobi equations (7), the leading-order asymptotic behaviour of the mean time to extinction  is given by

 limN→∞lnτN = V(0)−V(y∗)=U(0)−U(θ∗), (8)

where is the endemic equilibrium point of the deterministic system (5), and is the (assumed unique) non-zero equilibrium point of the complementary system

 dθdt = −∂H∂y∣∣∣y=0. (9)

Note that system (5) may be recovered as .

The solutions , to equations (7) are related via the Legendre transform; that is,

 U(θ)=supy{yTθ−V(y)}, V(y)=supθ{θTy−U(θ)},

see [23].

When (as is usually the case) it is not possible to find an analytical solution to either of the Hamilton-Jacobi equations (7), they may be solved numerically using the method of characteristics. That is, we write down the following -dimensional system of ordinary differential equations:

 dydt=∂H∂θ=∑l∈LlWl(y)eθTl,dθdt=−∂H∂y=−∑l∈L∂Wl∂y(eθTl−1),⎫⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪⎭ (10)

referred to as the ‘equations of motion’ of the system, and apply an appropriate numerical solver to (10). We then have , where is the ‘action’ integral,

 A=∫∞−∞θTdydtdt=−∫∞−∞yTdθdtdt, (11)

the integral in each case being evaluated along a trajectory from to . Note that .

Having set out the general Hamiltonian approach, we will now apply this technique to the infection model described in section 2 above.

## 4 Asymptotic persistence time formulae

Recall the infection model described in section 2, with transition rates given in table 1. In the large population limit, the scaled infection process converges almost surely, over finite time intervals, to the deterministic process satisfying the system of ordinary differential equations (5); that is,

 dyidt = β(k∑j=1λjyj)μi(fi−yi)−γyi for i=1,2,…,k. (12)

For there is a unique non-zero equilibrium point of the system (12), and it is globally asymptotically stable [20]. This endemic equilibrium point is given by [24]

 y∗i = μifiD(λ,μ)1+μiD(λ,μ) for i=1,2,…,k, (13)

where is the unique positive solution of

 βγk∑j=1μjfjλj1+μjD(λ,μ) = 1 (14)

The Hamiltonian (6) corresponding to the process is

 H(y,θ) = β(k∑j=1λjyj)(k∑i=1μi(fi−yi)(eθi−1))+γk∑i=1yi(e−θi−1).

The corresponding equations of motion (10) are, for ,

 dyidt = β(k∑j=1λjyj)μi(fi−yi)eθi−γyie−θi, (15) dθidt = −βλik∑j=1μj(fj−yj)(eθj−1)+β(k∑j=1λjyj)μi(eθi−1)−γ(e−θi−1). (16)

The non-zero equilibrium point given by (9) satisfies

 βλik∑j=1μjfj(eθ∗j−1)+γ(e−θ∗i−1) = 0 for i=1,2,…,k. (17)

Setting , then (17) implies that

 e−θ∗i = 1+λiB for i=1,2,…,k.

Substituting back into equation (17), we find that either (corresponding to ) or . The elements of are thus

 θ∗i = −ln(1+λiD(μ,λ)) for % i=1,2,…,k.

So far, we have allowed for heterogeneities in both infectivity and susceptibility simultaneously. If we restrict to only one type of heterogeneity, then it becomes possible to find an explicit formula for the action . Our main result is the following.

Theorem 1. Consider the heterogeneous SIS infection model defined in section 2, with transition rates given in table 1, and suppose . Recall that denotes the mean time from quasi-stationarity to disease extinction, and that is defined to be the unique positive solution of equation (14).

1. If heterogeneity is in infectivity alone (), then

 limN→∞lnτN = k∑i=1filn(1+λiD(1,λ))−γβD(1,λ). (18)
2. If heterogeneity is in susceptibility alone (), then

 limN→∞lnτN = k∑i=1filn(1+μiD(1,μ))−γβD(1,μ).

(Note: under the network interpretation, the assumption corresponds to every individual having the same in-degree, whereas corresponds to every individual having the same out-degree.)

Proof.

1. Suppose that for all , and consider a trajectory along which for . Along such a trajectory, the Hamiltonian simplifies to

 H(y,θ) = β(k∑j=1λjyj)(k∑i=1(fi−yi)(11+λiz−1))+γk∑i=1yiλiz = γz(k∑j=1λjyj)(1−βγk∑i=1fiλi1+λiz+βγk∑i=1yiλi1+λiz).

Since and (except at endpoints of the trajectory) the Hamilton-Jacobi equation reduces to

 1−βγk∑i=1fiλi1+λiz+βγk∑i=1λi1+λiz∂U∂θi = 0 (19)

Now along the trajectory under consideration, we have

 dUdz = −k∑i=1λi1+λiz∂U∂θi

and so equation (19) becomes

 dUdz = γβ−k∑i=1fiλi1+λiz ⇒U(θ(z))−U(θ(0)) = (γβ)z−k∑i=1∫z0fiλi1+λixdx ⇒U(θ)−U(0) = (γβ)z−k∑i=1filn(1+λiz) along θi=−ln(1+λiz).

The action is therefore given by

 A=U(0)−U(θ∗) = k∑i=1filn(1+λiD(1,λ))−γβD(1,λ),

as required.

2. For the SIS model on a finite network, in which each individual makes contact with each other individual at rate , it is known that, provided infectious periods are exponentially distributed, the decay parameter of the process is unchanged under transposition of the matrix of infection rates . This follows from the property of ‘network duality’, see [28, 18, 16]. In our context, this implies that the mean time to extinction from quasi-stationarity, , is identical if we interchange the roles of . Hence part (ii) of the theorem follows immediately from part (i).

We can confirm this as follows. With , the Hamiltonian may be written as

 H(y,θ) = γk∑i=1(eθi−1)(βγ(k∑j=1yj)μi(fi−yi)−yie−θi).

With the convention that when , take

 V(y) = k∑i=1yi(1+lnyi−ln(βγμi))−(k∑i=1yi)ln(k∑i=1yi) (20) +k∑i=1(fi−yi)ln(fi−yi).

Then

 ∂V∂yi = ln⎛⎜ ⎜⎝yiβγμi(fi−yi)(∑jyj)⎞⎟ ⎟⎠ for i=1,2,…,k,

and so

 H(y,∂V∂y) = 0

That is, satisfies the relevant Hamilton-Jacobi equation. The action is then given by

 A=V(0)−V(y∗) = k∑i=1filn(1+μiD(1,μ))−γβD(1,μ).

As expected, we recover the formula for the case of heterogeneous infectivity, but with the roles of interchanged.

Having found the solution for the case , we can find the corresponding function as the Legendre transform of . For , we have

 ddyi(yTθ−V(y)) = θi−ln⎛⎜ ⎜⎝yiβγμi(fi−yi)(∑jyj)⎞⎟ ⎟⎠,

and so any stationary point satisfies

 yi = βγμifi(∑jyj)eθi1+βγμi(∑jyj)eθi. (21)

For , define the function to be the solution of

 βγk∑j=1μjfje−θj+μjQ(μ,θ) = 1

Setting and substituting from (21) into the definition of , we find that , and hence the function has a stationary point at

 yi = μifieθiQ(μ,θ)1+μieθiQ(μ,θ). (22)

Evaluating the function at the point (22), we find

 U(θ) = k∑i=1filn(1+μieθiQ(μ,θ))−γβQ(μ,θ) (23)

and can easily verify that the function (23) does indeed satisfy . Now and , so we once again find that

 A = U(0)−U(θ∗)=k∑i=1filn(1+μiD(1,μ))−γβD(1,μ).

Although we did not actually need to find the functions in order to prove theorem 1(ii), we include them because knowledge of these functions can be of assistance in generalising and extending our results. We will demonstrate this in theorem 3 below.

Figure 1 illustrates theorem 1 in the case of groups with heterogeneity in infectivity (the graph for the corresponding case with heterogeneity in susceptibility is identical, by network duality). The exact value of is computed from equation (3) for total population sizes . The action  is computed from equation (18). For comparison, we also show the action  computed for the homogeneous population SIS model with the same value for . We see that formula (18) gives a good approximation to for population sizes from around upwards. We can also see that if we were to ignore heterogeneity and use the homogeneous population result, we would drastically over-estimate the persistence time of infection. We demonstrate this point in theorem 2 below, as well as comparing heterogeneous populations of greater or lesser degrees of heterogeneity.

## 5 The effect of increasing heterogeneity

Using the formulae of the previous section, we are now in a position to investigate the effect of increasing heterogeneity upon the persistence time of infection. First, in order to compare different levels of heterogeneity, we recall the definition of majorization [22]. For any , denote by the ordered components of . Then for , we say is majorized by , denoted , if and for . An equivalent definition is that for all convex functions . Intuitively, is ‘more heterogeneous’ than . More generally, given a probability vector (with components summing to 1) , then is -majorized by , written , if there exists a permutation such that and with and for .

Theorem 2. Consider two populations, with , , and each having the same group structure , where we use superscripts to denote the population under consideration. Recall that  denotes the mean time from quasi-stationarity to disease extinction.

• With heterogeneity in infectivity alone,

 λ(1)≺fλ(2) ⇒ limN→∞lnτ(1)N≥limN→∞lnτ(2)N.
• With heterogeneity in susceptibility alone,

 μ(1)≺fμ(2) ⇒ limN→∞lnτ(1)N≥limN→∞lnτ(2)N.

In particular, provided heterogeneity is in either infectivity or susceptibility but not both, then is maximised in the homogeneous case.

Proof. Consider the case of heterogeneity in infectivity, and suppose that . The function is concave for , and so applying proposition 14.A.3 of [22],

 k∑i=1fiλ(2)i1+λ(2)iD(1,λ(1)) ≤ k∑i=1fiλ(1)i1+λ(1)iD(1,λ(1))=γβ,

the final equality coming from the definition (14) of for population 1. The expression is a decreasing function of , and so from equation (14) for population 2 it follows that .

Now define the function . Then

 dψdz = ∑ifiλ(1)i1+λ(1)iz−γβ,

so that for , and hence . That is,

 k∑i=1filn(1+λ(1)iD(1,λ(2)))−γD(1,λ(2))β ≤ A(1). (24)

The function is concave for , and so again applying proposition 14.A.3 of [22],

 k∑i=1filn(1+λ(2)iD(1,λ(2))) ≤ k∑i=1filn(1+λ(1)iD(1,λ(2))).

Combining with (24) yields

 A(2) ≤ A(1),

and the result follows.

Part (ii) of the theorem follows immediately by interchanging the roles of .

Figure 2 illustrates theorem 2, as well as showing the effect of allowing heterogeneity in both infectivity and susceptibility simultaneously, for the case of equal-sized groups (). The constraints on the elements of in this case reduce to , and so we plot the action as a function of . We choose to keep fixed, with the value of being varied in order to achieve this. With both heterogeneities present, we have no explicit formula for the action , and instead compute it by first solving the equations of motion (15,16) numerically using the Matlab bvp4c command, and then integrating the numerical solution along the trajectory, equation (11). The solid contours in figure 2 show the action values  computed in this way. Note that the transformation here amounts to simply re-labelling the groups, so that figure 2 is invariant under a rotation of half a turn around the point ; also, we know from network duality that the action is unchanged under the transformation , so that figure 2 is invariant under reflection in the line .

For comparison, the dotted contours in figure 2 were computed by solving the eigenvalue equation (3) numerically for and , and assuming (without proof) that asymptotic formula (1) is valid for our model. Denoting by the mean time from quasi-stationarity to disease extinction in a population of size , formula (1) implies that the action  may be approximated by

 (ln(τ500√500)−ln(τ400√400))/100, (25)

and the dotted contours show computed values of formula (25). The fact that the dotted contours closely follow the solid contours provides some confirmation both that the action  gives a good approximation to for population sizes above , and that formula (1) does indeed apply to our model.

We see from figure 2 that the action decreases as we move away from the point , not only along the lines and , as ensured by theorem 2, but in any direction. That is, heterogeneity in infectivity, or susceptibility, or any combination of the two, reduces the value of compared to the homogeneous case. We discuss this further in section 7 below.

## 6 Generalising the infectious period distribution

So far, we have made the conventional assumption that individuals’ infectious periods are exponentially distributed. This is purely a mathematical convenience, not motivated by biological realism. Realism can be greatly improved by allowing infectious periods to follow an Erlang distribution, using the ‘method of stages’. That is, when an individual becomes infected, it passes through infectious stages, remaining in each stage for an exponentially distributed time of mean , before returning to susceptibility. As before, we denote by the (constant) number of individuals in group , and by the total population size. Denoting by the number of group  individuals in infectious stage  at time , then is a continuous-time Markov chain with transition rates given in table 2. The number of susceptible individuals in group  is .