Persistence time of SIS infections in heterogeneous populations and networks
Abstract
For a susceptibleinfectioussusceptible (SIS) infection model in a heterogeneous population, we present simple formulae giving the leadingorder 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 fadeout; 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 statespace Markov chain models, the expected persistence time for an infection that has become endemic in the population (i.e. starting from quasistationarity) 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 timeconsuming, 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
(1) 
where the values of depend upon parameters of the model, but not upon . We assume here that the process is supercritical, so that longterm endemicity is possible.
For the classic susceptibleinfectioussusceptible (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 supercriticality (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 leadingorder 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. Predating the above work, van Herwaarden and Grasman [26] showed that relationship (1) holds true for a particular suceptibleinfectiousremoved (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
(2) 
Equation (2) is not as precise as relationship (1), but does at least give the leadingorder 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 realworld infections. For instance, for a number of infections (eg SARS) it has been hypothesised that there exists a subgroup of ‘superspreaders’ 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 highdimensional) 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 largepopulation 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 homogeneouspopulation 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 continuoustime 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
Event  State transition  Transition rate 

Infection in group  
Recovery in group 
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 indegree and outdegree 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 halfedges, with stubs pointing inwards and stubs pointing outwards. Inwardpointing stubs are then paired with outwardpointing stubs throughout the population, to create links between individuals. In order that this process can produce a valid network, with no left over halfedges, 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 socalled ‘annealed’ network approximation is valid for an ensemble of such networks. This is a meanfield 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 nonzero 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
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
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 continuoustime Markov process on finite statespace 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 quasistationary distribution such that, for any initial state within ,
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
(3) 
where is the eigenvalue of with largest real part. The time to extinction from quasistationarity is exponentially distributed with mean .
Although may be computed exactly from equation (3), this can become impractical when the statespace 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 densitydependent process in the sense of chapter 11 of [15]; that is, the transition rates are of the form
(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
(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 diseasefree 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
(6) 
This Hamiltonian determines the following two complementary HamiltonJacobi partial differential equations:
(7) 
Each of these HamiltonJacobi 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 HamiltonJacobi equations (7), the leadingorder asymptotic behaviour of the mean time to extinction is given by
(8) 
where is the endemic equilibrium point of the deterministic system (5), and is the (assumed unique) nonzero equilibrium point of the complementary system
(9) 
Note that system (5) may be recovered as .
When (as is usually the case) it is not possible to find an analytical solution to either of the HamiltonJacobi 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:
(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,
(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,
(12) 
For there is a unique nonzero equilibrium point of the system (12), and it is globally asymptotically stable [20]. This endemic equilibrium point is given by [24]
(13) 
where is the unique positive solution of
(14) 
The Hamiltonian (6) corresponding to the process is
The corresponding equations of motion (10) are, for ,
(15)  
(16) 
The nonzero equilibrium point given by (9) satisfies
(17) 
Setting , then (17) implies that
Substituting back into equation (17), we find that either (corresponding to ) or . The elements of are thus
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 quasistationarity to disease extinction, and that is defined to be the unique positive solution of equation (14).

If heterogeneity is in infectivity alone (), then
(18) 
If heterogeneity is in susceptibility alone (), then
(Note: under the network interpretation, the assumption corresponds to every individual having the same indegree, whereas corresponds to every individual having the same outdegree.)
Proof.

Suppose that for all , and consider a trajectory along which for . Along such a trajectory, the Hamiltonian simplifies to
Since and (except at endpoints of the trajectory) the HamiltonJacobi equation reduces to
(19) The action is therefore given by
as required.

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 quasistationarity, , 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
With the convention that when , take
(20) Then
and so
That is, satisfies the relevant HamiltonJacobi equation. The action is then given by
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
and so any stationary point satisfies
(21) For , define the function to be the solution of
Setting and substituting from (21) into the definition of , we find that , and hence the function has a stationary point at
(22)
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 overestimate 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 quasistationarity to disease extinction.

With heterogeneity in infectivity alone,

With heterogeneity in susceptibility alone,
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],
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
so that for , and hence . That is,
(24) 
The function is concave for , and so again applying proposition 14.A.3 of [22],
Combining with (24) yields
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 equalsized 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 relabelling 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 quasistationarity to disease extinction in a population of size , formula (1) implies that the action may be approximated by
(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 continuoustime Markov chain with transition rates given in table 2. The number of susceptible individuals in group is .
Event  State transition  Transition rate 

Infection in group  
Transition to next infectious stage  for  
Recovery in group 