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.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 N the typical size of the population, the expected time from endemicity to extinction, τ , is asymptotically given by an expression of the form where the values of A, C depend upon parameters of the model, but not upon N .
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 A and C in terms of the basic reproduction number R 0 (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, R 0 > 1); specifically, A = (1/R 0 )− 1 + ln R 0 and C = R 0 √ 2π/(R 0 − 1) 2 , 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 A = (1/R 0 ) − 1 + ln R 0 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 [26] showed that relationship (1) holds true for a particular suceptible-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 [26] to evaluate A 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 lim Equation ( 2) is not as precise as relationship (1), but does at least give the leading-order behaviour of τ in the large population (N → ∞) 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 [12].
In contrast to almost all previous work, we are able to find an explicit formula for the constant A 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 ln τ 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 ln τ , 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.

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 (i = 1, 2, . . ., k) consisting of N i individuals.Denote by Recovery in group j I j → I j − 1 γI j Table 1: Transition rates for the k-group SIS model.
proportion of the population belonging to group i, so that i f i = 1.When a group i individual becomes infected, it remains so for a time distributed as an exponential random variable with mean 1/γ (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 j = 1, 2, . . ., k at the points of a Poisson process of rate βλ i µ j /N , where β is some overall measure of infectiousness, λ i represents the infectivity of group i individuals, and µ j 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 λ i , µ j values so that i λ i f i = j µ j f j = 1.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 {I(t) = (I 1 (t), I 2 (t), . . ., I k (t)) : t ≥ 0} is a continuous-time Markov chain with transition rates given in table 1, and the number of susceptible individuals in group j at time t ≥ 0 is S j (t) = N j − I j (t).We will assume throughout that β, γ > 0, and that f i , λ i , µ i > 0 for all i.Note that our model is a special case of the model of [8], although we use slightly different notation here.The basic reproduction number R 0 is given by the dominant eigenvalue of the matrix M with entries m ij = βλ i µ j f j , so that 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 d in and out-degree d out according to some joint probability mass function p d in , d out on Z 2 + .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 d in stubs pointing inwards and d out 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 E d in = E [d out ].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 k of d in , d out pairs having non-zero probability, and define a bijective function c d in , d out : Z 2 + → {1, 2, . . ., k} that assigns a unique number to each of the possible d in , d out pairs.We say an individual is of 'group j' if they have degrees d in , d out = c −1 (j), and define d in (j), d out (j) to be the in and out degrees, respectively, of a group j individual.For j = 1, 2, . . ., k, denote by N j the total number of group j individuals in the population, and by I j (t) 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 by When individuals become infected, they remain so for an exponentially distributed time of mean 1/γ 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 taking 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.

General theory regarding persistence time from endemicity
Consider an infection modelled by a continuous-time Markov process {X(t) : t ≥ 0} on finite state-space S ⊂ Z k 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 Q C 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 [10] there exists a unique quasi-stationary distribution q = {q x : x ∈ C} such that, for any initial state within C, That is, provided the infection does not die out, it will settle to the endemic distribution q.The distribution q may be found as the unique solution of where −(1/τ ) is the eigenvalue of Q C 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 X(t) is a density-dependent process in the sense of chapter 11 of [15]; that is, the transition rates are of the form for some functions W l : R k → R + , where L is the set of possible jumps from each state x ∈ S 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 ( [15], Theorem 11.2.1), the scaled process X(t)/N converges almost surely over finite time intervals, as N → ∞, to the solution y(t) of the ordinary differential equation system For our application, we suppose that the system (5) possesses two equilibrium points: a stable endemic equilbrium point y * with all components strictly positive, and an unstable disease-free equilibrium point at y = 0. 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 This 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 N → ∞ (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 lim where y * is the endemic equilibrium point of the deterministic system (5), and θ * is the (assumed unique) non-zero equilibrium point of the complementary system Note that system (5) may be recovered as dy dt = ∂H ∂θ θ=0 .The solutions U (θ), V (y) to equations ( 7) are related via the Legendre transform; that is, 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 2kdimensional 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 lim N →∞ (ln τ )/N = A, where A is the 'action' integral, the integral in each case being evaluated along a trajectory from (y * , 0) to (0, θ * ).Note that Having set out the general Hamiltonian approach, we will now apply this technique to the infection model described in section 2 above.

Asymptotic persistence time formulae
Recall the infection model {I(t) : t ≥ 0} described in section 2, with transition rates given in table 1.In the large population limit, the scaled infection process I(t)/N converges almost surely, over finite time intervals, to the deterministic process y(t) satisfying the system of ordinary differential equations (5); that is, For R 0 > 1 there is a unique non-zero equilibrium point y * of the system (12), and it is globally asymptotically stable [20].This endemic equilibrium point y * is given by [24] where D(λ, µ) is the unique positive solution of The Hamiltonian (6) corresponding to the process I(t) is The corresponding equations of motion (10) are, for i = 1, 2, . . ., k, The non-zero equilibrium point θ * given by ( 9) satisfies Setting B = (β/γ) j f j µ j 1 − e θ * j , then (17) implies that Substituting back into equation ( 17), we find that either B = 0 (corresponding to θ = 0) or B = D(µ, λ).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 A. 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 R 0 > 1. Recall that τ denotes the mean time from quasi-stationarity to disease extinction, and that D(λ, µ) is defined to be the unique positive solution of equation ( 14).
(i) If heterogeneity is in infectivity alone (µ = 1), then lim (Note: under the network interpretation, the assumption µ = 1 corresponds to every individual having the same in-degree, whereas λ = 1 corresponds to every individual having the same out-degree.)Proof.
(i) Suppose that µ i = 1 for all i, and consider a trajectory {θ(z Along such a trajectory, the Hamiltonian simplifies to Since z > 0 and j λ j y j > 0 (except at endpoints of the trajectory) the Hamilton-Jacobi equation H ∂U ∂θ , θ = 0 reduces to Now along the trajectory under consideration, we have and so equation (19) becomes The action is therefore given by (ii) For the SIS model on a finite network, in which each individual u makes contact with each other individual v at rate β uv , 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 {β uv }.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 λ = 1, the Hamiltonian may be written as With the convention that y ln y = 0 when y = 0, take Then and so That is, V (y) 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 V (y) for the case λ = 1, we can find the corresponding function U (θ) as the Legendre transform of V (y).For i = 1, 2, . . ., k, we have and so any stationary point satisfies For θ ∈ R k , define the function Q(µ, θ) to be the solution of Setting R = (β/γ) j y j and substituting from (21) into the definition of R, we find that R = Q(µ, θ), and hence the function y T θ − V (y) has a stationary point at Evaluating the function y T θ − V (y) at the point (22), we find and can easily verify that the function (23) does indeed satisfy H ∂U ∂θ , θ = 0. Now Q(µ, 0) = D(1, µ) and Q(µ, θ * ) = 0, so we once again find that Although we did not actually need to find the functions U (θ), V (y) 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 k = 2 groups with heterogeneity in infectivity (the graph for the corresponding case with heterogeneity in susceptibility is identical, by network duality).The exact value of (ln τ )/N is computed from equation (3) for total population sizes N = 100, 150, . . ., 650.The action A is computed from equation (18).For comparison, we also show the action A 0 = (1/R 0 ) − 1 + ln R 0 computed for the homogeneous population SIS model with the same value for R 0 .We see that formula (18) gives a good approximation to (ln τ )/N for population sizes from around N = 300 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.
(ii) With heterogeneity in susceptibility alone, In particular, provided heterogeneity is in either infectivity or susceptibility but not both, then lim N → ∞ (ln τ ) /N is maximised in the homogeneous case.

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 k = 2 equal-sized groups (f 1 = f 2 = 0.5).The constraints on the elements of λ, µ in this case reduce to λ 1 + λ 2 = µ 1 + µ 2 = 2, and so we plot the action as a function of (λ 1 , µ 1 ) ∈ (0, 2) 2 .We choose to keep R 0 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, equation (11).The solid contours in figure 2 show the action values A computed in this way.Note that the transformation (λ 1 , µ 1 ) → (2 − λ 1 , 2 − µ 1 ) here amounts to simply re-labelling the groups, so that figure 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 (λ 1 , µ 1 ) → (µ 1 , λ 1 ), so that figure 2 is invariant under reflection in the line λ 1 = µ 1 .
For comparison, the dotted contours in figure 2 were computed by solving the eigenvalue equation ( 3) numerically for N = 400 and N = 500, and assuming (without proof) that asymptotic formula ( 1) is valid for our model.Denoting by τ N 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 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 A gives a good approximation to (ln τ )/N for population sizes above N = 400, 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 (λ 1 , µ 1 ) = (1, 1), not only along the lines λ 1 = 1 and µ 1 = 1, 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 lim N →∞ (ln τ )/N compared to the homogeneous case.We discuss this further in section 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 (sγ) −1 , before returning to susceptibility.As before, we denote by N j the (constant) number of individuals in group j, and by N = y is e −θis − 1 .
It is immediate from equation ( 17) of [7] that the deterministic endemic equilibrium point is given by where y * i is the solution (13) for the model with exponentially distributed infectious periods (s = 1).
It is straightforward to show that the elements of θ * are given by where D s (µ, λ) is the solution of For the SIS model with Erlang-distributed infectious periods in a homogeneous population (k = 1), the solution U (θ) to the relevant Hamilton-Jacobi equation was found in [9] to be Taking the Legendre transform, we find that V (y) for this homogeneouspopulation model is given by Comparing solution (27) 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 V (y) for the SIS model with Erlangdistributed infectious periods and heterogeneous susceptibilities, and verify that the relevant Hamilton-Jacobi equation is indeed satisfied.The solution is thus found to be Taking the Legendre transform, we find where Q s (µ, θ) is the solution of The action A in this case is thus as 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 k = 2 equal-sized groups (f 1 = f 2 = 0.5).In constructing this figure, we have assumed (without proof) that asymptotic formula (1) is valid for our model.Consequently, we plot the function 1  2 ln N + ln τ , which according to formula (1) should, as N increases, approach a straight line of gradient A and intercept ln C. 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 (λ = (5/3, 1/3), µ = (1, 1)) or heterogeneous suceptibility = (1, 1), µ = (5/3, 1/3)).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 s → ∞.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 t 0 (the burn-in period) and t max such that (i) by time t 0 the state of the process is approximately quasi-stationary (having started at time zero close to the re-scaled deterministic equilibrium point N y * , and conditioning upon survival to time t 0 ); and (ii) by time t max 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 T 1 , T 2 , . . ., T r the extinction times of those simulations that went extinct within the time window (t 0 , t max ), and by m the number of simulations that had not gone extinct by time t max , our estimate is We have included in figure 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 figure 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 (ln τ )/N for population sizes above N = 200 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 [4] for the homogeneous population case.7 Discussion and possible extensions The main result of this paper, theorem 1, provides a simple explicit formula for lim N →∞ (ln τ )/N , 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 [1] or with arbitrary infectious period distribution [4].Theorem 1 thus represents a significant advance, but many open questions remain.
Firstly, for the SIS model in a homogeneous population, both of [1,4] 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 [2] to be valid (and formulae given for the constant C) for general 1-dimensional processes of bounded jump size.The technique of [2] 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 it is not clear how the approach of [2] may be extended to multi-dimensional processes.
Secondly, our model as described in section 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.In particular, this 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.
In comparing two populations, we need to define our 'groups' slightly differently than in section 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 d (ii) for two populations with κ (1) = κ (2) , γ (1) = γ (2) , d out ≤ cv d out ⇒ A out ≥ A (ii) for two populations with κ (1) = κ (2) , γ (1) = γ (2) , d in in ; in particular, A in is maximised when every individual has the same in-degree, P d in = µ = 1.
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.We conjecture that theorem 1(i) remains valid, and hence also theorem 2(i), if infectious periods are allowed to follow an Erlang, rather than exponential, distribution.Indeed, figure 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 U (θ), V (y) to the relevant Hamilton-Jacobi equations even for the case of exponentially distributed infectious periods, but only, in proving theorem 1(i), to evaluate U (θ) 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 [8] regarding the effect of such heterogeneities upon the (large population) mean endemic prevalence level y * = k i=1 y * i .Theorem 7(i) and theorem 10 of [8] show, respectively, that heterogeneous infectivity alone has no effect upon the endemic prevalence level y * , whereas heterogeneous suceptibility alone can only decrease y * , with µ (1) ≺ f µ (2) ⇒ y * (1) ≥ y * (2) , corresponding to our theorem 2(ii) for persistence times.When both types of heterogeneity are combined, theorem 8 of [8] shows that if infectivity and susceptibility are non-negatively correlated ( k i=1 λ i µ i f i ≥ 1) then the endemic prevalence level cannot be greater than for a homogeneous population with the same R 0 value.However, theorem 9 of [8] together with numerical work shown in figure 3 of [8] 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 R 0 values matched).This presents an interesting contrast to our numerical results in figure 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 [8], via an Ornstein-Uhlenbeck diffusion approximation that leads to a multivariate normal approximation to the quasi-stationary distribution q.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 [11,9].The methods of the current paper, in contrast, deal directly with the expected persistence time and yield correct leading-order asymptotic formulae.

A Appendix
The Hamilton-Jacobi equations (7) amount to two alternative ways of expressing the eigenvector equation (3), retaining only terms of leading order in N in the limit as N → ∞.We briefly outline the derivation of each equation, referring the reader to [3] and references therein for full details.
For a process with transition rates of the form (4), equation (3) may be written as l∈L q x−l W l x − l N − q x W l x N = −(τ N ) −1 q x for x ∈ C. (29) Suppose (ansatz) that q x = exp (−N V (x/N ) + o(N )) for some function V (•).Writing y = x/N , then V y − l N = V (y) − (l T /N ) ∂V ∂y + o(1/N ) and hence Similarly, W l y − l N = W l (y) + o(1/N ).Retaining only terms of leading order in N , and assuming that τ is sufficiently large for the right hand side of equation (29) to be neglected, then equation (29) reduces to l∈L W l (y) exp l T ∂V ∂y − 1 = 0.

P
condition f(1) = f(2) = f 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 d in , d out , the degrees of distinct individuals being mutually independent, with E d in = E [d out ] = µ and d in , d out ≤ d max for some d max ∈ N. 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 1/γ before returning to the susceptible state.Recall that τ denotes the expected time from quasi-stationarity to extinction of infection.(a) Suppose that P d in = µ = 1, so every individual has the same in-degree, and that infectious periods are exponentially distributed.(d out = i) ln (1 + iD out ) A out is maximised when every individual has the same out-degree, P (d out = µ) = 1.(b)Suppose that P (d out = µ) = 1, so every individual has the same outdegree, and that infectious periods follow an Erlang distribution.

Table 2 :
Transition rates for the k-group, s-stage SIS model.total population size.Denoting by I jv (t) the number of group j individuals in infectious stage v at time t, then {I jv (t) : j = 1, 2, ..., k, v = 1, 2, ..., s, t ≥ 0} is a continuous-time Markov chain with transition rates given in table2.The number of susceptible individuals in group j is S j (t) = N j − s v=1 I jv (t).Writing y = {y iv