Persistence time of SIS infections in heterogeneous populations and networks.

For a susceptible-infectious-susceptible 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 (1) a population with heterogeneity in individuals' susceptibility and/or infectiousness; or (2) 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.

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 (van Doorn and Pollett 1993). 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 (Andersson and Djehiche 1998; Ball et al. 2016; van Herwaarden and Grasman 1995) that, denoting by N the typical size of the population, the expected time from endemicity to extinction, , is asymptotically given by an expression of the formwhere the values of A, C depend upon parameters of the model, but not upon N. It is assumed 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 (1971), Andersson and Djehiche (1998) found simple explicit expressions for both A and C 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 et al. (2016) to allow for a general infectious period distribution in place of the exponential distribution assumed by Andersson and Djehiche (1998); they showed that leading-order behaviour is unchanged, so that as before, while the value of C 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 (1995) showed that relationship (1) holds true for a particular susceptible–infectious–removed (SIR) infection model. In this case, however, evaluation of the constant A requires numerical solution of a system of ordinary differential equations, while no method for evaluating C is given.

The system of ordinary differential equations used in van Herwaarden and Grasman (1995) to evaluate A may be regarded as the equations of motion corresponding to a particular Hamiltonian system. More recently, a number of authors (Assaf and Meerson 2010, 2017; Dykman et al. 1994; Elgart and Kamenev 2004; Kamenev and Meerson 2008; Lindley et al. 2014) have applied this Hamiltonian approach to a range of infection models to derive results of the formEquation (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 A 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 (Dorogovtsev et al. 2008).

In contrast to almost all previous work, we are able to find an explicit formula for the constant A in Eq. (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 Sect. 2, we define precisely our heterogeneous population SIS model, and describe how it may be interpreted as approximating a directed network model. In Sect. 3 we recall some general theory that will be required in the sequel. Our main result, Theorem 1, is derived in Sect. 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 Sect. 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, we present an approximate formula for , valid provided is only slightly greater than 1, together with numerical work (Fig. 2) indicating that mean persistence time is again maximised in the homogeneous-population case. In Sect. 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 (Fig. 3) suggesting that this is also true when instead heterogeneity is in infectivities. Finally, in Sect. 7, we discuss the directed network interpretation of our results (Theorem 4) and suggest some directions for further work.

Fig. 2

Contour plot showing the action A (solid contours) as a function of . Fixed parameter values , , . The action is maximised at , with value . Dashed contours show the approximation  computed from formula (25). Dotted contours show a finite-population approximation—see main text for details

Fig. 3

The effect of the infectious period distribution upon the mean persistence time of infection . Fixed parameter values , , , . Solid line (‘Action’) has gradient given by Eq. (18) with , , intercept chosen arbitrarily; dashed line (‘Exponential infectious periods’) computed from Eq. (3) with , ; crosses (‘Heterogeneous susceptibilities’) computed via simulation with , and constant infectious periods; circles (‘Heterogeneous infectivities’) computed via simulation with , and constant infectious periods

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 N individuals divided into k groups, with group i () consisting of individuals. Denote by the proportion of the population belonging to group i, so that . When a group i 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 i 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 i individuals, and represents the susceptibility of group j individuals. (The assumption that the group i to group j 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. Denoting by the number of infected individuals in group j at time t, then the process is a continuous-time Markov chain with transition rates given in Table 1, and the number of susceptible individuals in group j at time  is . We will assume throughout that , and that for all i. Note that our model is a special case of the model of Clancy and Pearce (2013), although we use slightly different notation here. The basic reproduction number is given by the dominant eigenvalue of the matrix M with entries , so that

Table 1

Transition rates for the k-group SIS model

Event	State transition	Transition rate
Infection in group j	I_j→I_j+1	βN∑m=1kλ_mI_mμ_j(N_j-I_j)
Recovery in group j	I_j→I_j-1	γI_j

We now describe how the above model may be interpreted as describing infection spreading through a network. Each of the N 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 Britton et al. (2007) and Chen and Olvera-Cravioto (2013) 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 Dorogovtsev et al. 2008 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 k 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 j’ if they have degrees , and define to be the in and out degrees, respectively, of a group j individual. For , denote by the total number of group j individuals in the population, and by the number of group j individuals who are infectious at time t. Then under the annealed network approximation, the total rate at which group j individuals become infected is given byWhen 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 k-group model with transition rates given in Table 1 by takingfor .

The undirected version of the above annealed network approximation (with ) has been studied by Hindes and Schwartz (2016), via numerical solution of Hamilton’s equations of motion, Eqs. (15,16) below. We note that real-world social networks are often only partially directed, containing both directed and undirected edges (that is to say, bi-directional edges are over-represented compared to what would be expected by chance), see Spricer and Britton (2015). Each individual then has three-dimensional degree , where is the number of undirected edges connected to the individual in question. The above annealed network approximation will remain valid if we replace with , provided, as before, that the resulting network is uncorrelated.

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

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 Q. Suppose that S is made up of an absorbing set of states A (corresponding to absence of disease) and a single transient communicating class C. We denote by the transition rate matrix restricted to C. Then the infection will almost surely die out (i.e. the process will leave C) within finite time, and (Darroch and Seneta 1967) there exists a unique quasi-stationary distribution such that, for any initial state within 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 ofwhere 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 Eq. (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 Ethier and Kurtz (2005); that is, the transition rates are of the formfor some functions , where L is the set of possible jumps from each state and N is some parameter indicating overall size of the system (in our applications, N will be the size of the population). Under mild technical conditions (Ethier and Kurtz 2005, Theorem 11.2.1), the scaled process converges almost surely over finite time intervals, as , to the solution of the ordinary differential equation systemFor our application, we suppose that the system (5) possesses two equilibrium points: a stable endemic equilibrium 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 (Assaf and Meerson 2017) and references therein.

The Hamiltonian of the system is defined to beThis Hamiltonian determines the following two complementary Hamilton–Jacobi partial differential equations: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 bywhere is the endemic equilibrium point of the deterministic system (5), and is the (assumed unique) non-zero equilibrium point of the complementary systemNote that system (5) may be recovered as .

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

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 2k-dimensional system of ordinary differential equations:referred to as the ‘equations of motion’ of the system, and apply an appropriate numerical solver to (10). We then have , where A is the ‘action’ integral,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 Sect. 2 above.

Asymptotic persistence time formulae

Recall the infection model described in Sect. 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,For there is a unique non-zero equilibrium point of the system (12), and it is globally asymptotically stable (Lajmanovich and Yorke 1976). This endemic equilibrium point is given by (Nold 1980)where is the unique positive solution ofThe Hamiltonian (6) corresponding to the process isThe corresponding equations of motion (10) are, for , The non-zero equilibrium point given by (9) satisfiesSetting , then (17) implies thatSubstituting back into Eq. (17), we find that either (corresponding to ) or . The elements of are thusSo 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 A. Our main result is the following.

Theorem 1

Consider the heterogeneous SIS infection model defined in Sect. 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 Eq. (14).

If heterogeneity is in infectivity alone (), then

If heterogeneity is in susceptibility alone (), then

(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

Suppose that for all i, and consider a trajectory along which for . Along such a trajectory, the Hamiltonian simplifies to Since and (except at endpoints of the trajectory) the Hamilton–Jacobi equation reduces to Now along the trajectory under consideration, we have and so Eq. (19) becomes The action is therefore given by as required.

For the SIS model on a finite network, in which each individual u makes contact with each other individual v 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 Wilkinson and Sharkey (2013), Holley and Liggett (1975) and Harris (1976). 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 With the convention that when , take Then and so That is, satisfies the relevant Hamilton–Jacobi 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 For , define the function to be the solution of Setting and substituting from (21) into the definition of R, we find that , and hence the function has a stationary point at Evaluating the function at the point (22), we find and can easily verify that the function (23) does indeed satisfy . Now and , so we once again find that

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 Eq. (3) for total population sizes . The action A is computed from Eq. (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.

Fig. 1

Values of and asymptotic formulae plotted against population size N. Fixed parameter values , , , , , . For these parameter values, , , , , with action  and corresponding action for the homogeneous case . The dots, labelled ‘eigenvalue computation’, are the true values of computed from Eq. (3); the action A is computed from Eq. (18); the homogeneous action is computed as

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 (Marshall et al. 2011). 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 (1), (2) to denote the population under consideration. Recall that  denotes the mean time from quasi-stationarity to disease extinction.In particular, provided heterogeneity is in either infectivity or susceptibility but not both, then is maximised in the homogeneous case.

With heterogeneity in infectivity alone,

With heterogeneity in susceptibility alone,

Consider the case of heterogeneity in infectivity, and suppose that . The function is concave for , and so applying proposition 14.A.3 of Marshall et al. (2011),the final equality coming from the definition (14) of for population 1. The expression is a decreasing function of z, and so from Eq. (14) for population 2 it follows that .

Now define the function . Thenso that for , and hence . That is,The function is concave for , and so again applying proposition 14.A.3 of Marshall et al. (2011),Combining with (24) yieldsand 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 A, 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, Eq. (11). The solid contours in Fig. 2 show the action values A computed in this way. Note that the transformation here amounts to simply re-labelling the groups, so that Fig. 2 is invariant under a rotation of half a turn around the point (1, 1); also, we know from network duality that the action is unchanged under the transformation , so that Fig. 2 is invariant under reflection in the line .

Although we do not have an explicit formula for A when both types of heterogeneity are present, we can obtain an approximate formula valid for close to (and above) 1, as follows. For only slightly larger than 1, the trajectory from to satisfying equations (15,16) may be approximated by a straight line. The integral (11) can easily be evaluated along the straight line connecting to , and we thus find that , whereValues of  are shown as dashed contours in Fig. 2, and we see that for these parameter values does indeed provide a reasonable approximation to the true value A. Note that although the action A is known (from network duality) to be invariant under the interchange of , the approximating formula (25) does not possess this symmetry. Nevertheless, we see that in Fig. 2 the contours corresponding to  do appear symmetrical under reflection in the line , so that the expected symmetry does at least hold approximately here.

For comparison, the dotted contours in Fig. 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 N, formula (1) implies that the action A may be approximated byand the dotted contours show computed values of formula (26). The fact that the dotted contours closely follow the solid contours provides some confirmation both that the action A gives a good approximation to for population sizes above , and that formula (1) does indeed apply to our model.

We see from Fig. 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 Sect. 7 below.

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 s 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 j, and by the total population size. Denoting by the number of group j individuals in infectious stage v at time t, then is a continuous-time Markov chain with transition rates given in Table 2. The number of susceptible individuals in group j is .

Table 2

Transition rates for the k-group, s-stage SIS model

Event	State transition	Transition rate
Infection in group j	I_j1→I_j1+1	βN∑m=1kλ_m∑v=1sI_mvμ_jN_j-∑v=1sI_jv
Transition to next infectious stage	I_j,v-1,I_jv→(I_j,v-1-1,I_jv+1)	sγI_j,v-1 for v=2,3,…,s
Recovery in group j	I_js→I_js-1	sγI_js

Writing and , the corresponding Hamiltonian isIt is immediate from equation (17) of Clancy (2015) that the deterministic endemic equilibrium point is given bywhere is the solution (13) for the model with exponentially distributed infectious periods ().

It is straightforward to show that the elements of are given bywhere is the solution ofFor the SIS model with Erlang-distributed infectious periods in a homogeneous population (), the solution to the relevant Hamilton–Jacobi equation was found in Clancy and Tjia (2018) to beTaking the Legendre transform, we find that for this homogeneous-population model is given byComparing solution (28) for the SIS model with Erlang-distributed infectious periods in a homogeneous population and solution (20) for the SIS model with exponentially distributed infectious periods and heterogeneous susceptibilities, one may now guess the form of the solution for the SIS model with Erlang-distributed infectious periods and heterogeneous susceptibilities, and verify that the relevant Hamilton–Jacobi equation is indeed satisfied. The solution is thus found to beTaking the Legendre transform, we findwhere is the solution ofThe action A in this case is thusas before, and the following result is immediate.

Theorem 3

Theorem 1(ii) remains valid if infectious periods are allowed to follow an Erlang, rather than exponential, distribution. Consequently, Theorem 2(ii) likewise remains valid with Erlang-distributed infectious periods.

Figure 3 illustrates the effect of the infectious period distribution in the case of equal-sized groups (). In constructing this figure, we have assumed (without proof) that asymptotic formula (1) is valid for our model. Consequently, we plot the function , which according to formula (1) should, as N increases, approach a straight line of gradient A and intercept . The dashed line, corresponding to exponentially distributed infectious periods, was computed using the eigenvalue characterisation (3). By network duality, this dashed line may be interpreted as corresponding to either heterogeneous infectivity (, ) or heterogeneous susceptibility (, ). We used Monte Carlo simulation to estimate the mean persistence time  with constant (non-random) infectious periods for the cases of heterogeneous infectivity and heterogeneous susceptibility separately. This infectious period distribution corresponds to an Erlang distribution with s stages in the limit as . An issue that arises is that the time until extinction of infection, starting from quasi-stationarity, is exponentially distributed with mean increasing exponentially in population size, so that to simulate the process to extinction can be very time-consuming. To get around this, we fixed times (the burn-in period) and such that (i) by time the state of the process is approximately quasi-stationary (having started at time zero close to the re-scaled deterministic equilibrium point , and conditioning upon survival to time ); and (ii) by time a substantial proportion of all simulations have reached extinction. We then estimated the mean time to extinction  using the maximum likelihood estimator. That is, denoting by the extinction times of those simulations that went extinct within the time window , and by m the number of simulations that had not gone extinct by time , our estimate isWe have included in Fig. 3 a solid line with gradient equal to the action A computed from formula (18). Note that the intercept of this line was chosen arbitrarily, since we have no way to evaluate the constant C in formula (1) for our model.

We see from Fig. 3 that the dashed line corresponding to exponentially distributed infectious periods does indeed appear to be a straight line of gradient A, providing some confirmation both that the action A gives a reasonable approximation to for population sizes above and that formula (1) is valid for our model. With constant infectious periods, we see that heterogeneous infectivity and heterogeneous susceptibility result in almost identical estimates of , and that these estimates lie close to a straight line of gradient A. It appears that, as with exponentially distributed infectious periods, the value of  is unchanged if we interchange . We therefore conjecture that, similarly to Theorem 3 above, Theorem 1(i) and Theorem 2(i) remain valid with Erlang-distributed infectious periods. The model with constant infectious periods has reduced mean persistence time  compared to the model with exponentially distributed infectious periods, but the difference is in the pre-factor constant C and not the leading-order constant A, in line with the results of Ball et al. (2016) for the homogeneous population case.

Discussion and possible extensions

The main result of this paper, Theorem 1, provides a simple explicit formula for , where is the expected time from endemicity to extinction for an SIS infection model with heterogeneity in either infectivity or susceptibility of individuals, in a population of size N. The only infection model for which such a formula has previously been available is the SIS model in a homogeneous population, either with exponentially distributed infectious periods (Andersson and Djehiche 1998) or with arbitrary infectious period distribution (Ball et al. 2016). Theorem 1 thus represents a significant advance, but many open questions remain.

Firstly, for the SIS model in a homogeneous population, both of Andersson and Djehiche (1998) and Ball et al. (2016) established asymptotic approximations for of the form (1), with explicit formulae for the pre-factor constant C. Our result is less precise than this; we have not shown that an asymptotic formula of the form (1) is valid for our model (although we conjecture, and have presented some numerical evidence, that this is the case), nor have we attempted to evaluate the pre-factor constant C. The asymptotic form (1) has been shown by Assaf and Meerson (2010) to be valid (and formulae given for the constant C) for general 1-dimensional processes of bounded jump size. The technique of Assaf and Meerson (2010) is an extension of the approach employed here, retaining terms beyond the leading order in N. The analysis is considerably more intricate than the leading-order treatment we have restricted ourselves to, and the extension of the approach of Assaf and Meerson (2010) to multi-dimensional processes such as the infection model considered here is the subject of ongoing work.

Secondly, our model as described in Sect. 2 incorporates heterogeneity in both infectivity and susceptibility simultaneously, but we have only been able to provide an explicit asymptotic formula for the cases in which only one of these two heterogeneities is present. When both heterogeneities are present, provided is only slightly above one then  given by formula (25) can be used to approximate , and in Fig. 2 we saw that can indeed provide a reasonable approximation to A. However, even a small error in A can correspond to a very large error in our estimate of , so that an approximate formula such as (25) is of considerably less utility than an exact formula such as (18). In particular, the fact that we have not been able to find an exact formula such as (18) valid when both types of heterogeneity are present severely restricts the class of networks to which our results may be applied under the annealed network approximation—we require either that every individual has the same in-degree, or that every individual has the same out-degree. Nevertheless, for the class of directed networks to which they apply our results represent an interesting step forward, and since network models are of great ongoing interest in infection modelling, we now present our results in a form suited to the network interpretation.

As a preliminary to the statement of our results, we require the concept of convex ordering of random variables, defined as follows (Shaked and Shanthikumar 2007, section 3.A.1). Given two random variables , then is greater than in the sense of convex ordering, denoted , if for all convex functions . If take values in , then an equivalent definition is thatNote that implies ; intuitively, is ‘more variable’ than .

In comparing two populations, we need to define our ‘groups’ slightly differently than in Sect. 2. Specifically, partition the population into groups in such a way that two individuals belong to the same group if they share the same values of both and ; the condition required by Theorem 2 is thus satisfied. (As before, superscripts (1), (2) denote the population under consideration.)

Theorem 4

Consider an SIS infection in a population of N individuals connected by an uncorrelated directed network. Each individual has in-degree and out-degree distributed as , the degrees of distinct individuals being mutually independent, with and for some . Infection transmits along each link from an infectious to a susceptible individual at rate , and when an individual becomes infected it remains so for a time of mean before returning to the susceptible state. Recall that denotes the expected time from quasi-stationarity to extinction of infection.

Suppose that , so every individual has the same in-degree, and that infectious periods are exponentially distributed. Then

, with where is the solution of

for two populations with , , in particular, is maximised when every individual has the same out-degree, .

Suppose that , so every individual has the same out-degree, and that infectious periods follow an Erlang distribution. Then

, with where is the solution of

for two populations with , , , in particular, is maximised when every individual has the same in-degree, .

Note that our asymptotic results (Theorem 1) are exact for the model with transition rates given by Table 1, but approximate under the annealed network interpretation.

A third open issue is to allow for more general infectious period distributions in the case of heterogeneous infectivity (or heterogeneous out-degree, under the network interpretation). We conjecture that Theorem 1(i) remains valid, and hence also Theorem 2(i) and Theorem 4(a), if infectious periods are allowed to follow an Erlang, rather than exponential, distribution. Indeed, Fig. 3 suggests that the mean persistence time  is unchanged when are interchanged. A difficulty here is that, in contrast to the case of heterogeneous susceptibilities, we have not been able to find complete solutions to the relevant Hamilton–Jacobi equations even for the case of exponentially distributed infectious periods, but only, in proving Theorem 1(i), to evaluate along one particular trajectory.

One advantage of simple asymptotic formulae such as provided by Theorem 1, as opposed to the exact formula (3), is that they provide a route to qualitative results such as Theorem 2, that increasing heterogeneity reduces (at least to leading order) the expected persistence time of infection, and in particular that persistence time is maximised in a homogeneous population. Theorem 2 establishes this ordering when heterogeneity is in either infectivity or susceptibility; figure 2 suggests that the result remains true even when both types of heterogeneity are present. It is interesting to compare with the results contained in section 5 of Clancy and Pearce (2013) regarding the effect of such heterogeneities upon the (large population) mean endemic prevalence level . Theorem 7(i) and theorem 10 of Clancy and Pearce (2013) show, respectively, that heterogeneous infectivity alone has no effect upon the endemic prevalence level , whereas heterogeneous susceptibility alone can only decrease , with , corresponding to our Theorem 2(ii) for persistence times. When both types of heterogeneity are combined, theorem 8 of Clancy and Pearce (2013) shows that if infectivity and susceptibility are non-negatively correlated () then the endemic prevalence level cannot be greater than for a homogeneous population with the same value. However, theorem 9 of Clancy and Pearce (2013) together with numerical work shown in figure 3 of Clancy and Pearce (2013) demonstrates that when infectivity and susceptibility are negatively correlated it is possible for the endemic prevalence level to be greater than in the homogeneous case (with values matched). This presents an interesting contrast to our numerical results in Fig. 2, where heterogeneities were found to always decrease (to leading order) the expected persistence time, regardless of whether are positively or negatively correlated.

There is a slightly counter-intuitive aspect to the above results, in that an increase in endemic prevalence level may correspond to a decrease in expected persistence time. This is easily resolved by observing that an increase in prevalence level may be accompanied by a corresponding increase in variability, leading to faster extinction of infection. The effect of heterogeneities upon the variability of the quasi-stationary distribution is studied in section 6 of Clancy and Pearce (2013), via an Ornstein–Uhlenbeck diffusion approximation that leads to a multivariate normal approximation to the quasi-stationary distribution . The variability of this approximating normal distribution is then used as a proxy measure of persistence time. This approach seems reasonable in terms of qualitative comparisons between infection models, and is common in the literature. However, the approach is known to give a very bad numerical approximation to mean persistence time, with incorrect leading-order asymptotic behaviour, due to the failure in the lower tail of the normal approximation to the quasi-stationary distribution (Doering et al. 2005; Clancy and Tjia 2018). The methods of the current paper, in contrast, deal directly with the expected persistence time and yield correct leading-order asymptotic formulae.