A Re-entrant Phase Transition in the Survival of Secondary Infections on Networks

We study the dynamics of secondary infections on networks, in which only the individuals currently carrying a certain primary infection are susceptible to the secondary infection. In the limit of large sparse networks, the model is mapped to a branching process spreading in a random time-sensitive environment, determined by the dynamics of the underlying primary infection. When both epidemics follow the Susceptible-Infective-Recovered model, we show that in order to survive, it is necessary for the secondary infection to evolve on a timescale that is closely matched to that of the primary infection on which it depends.


Introduction
Superinfections are a major cause of global mortality and morbidity. For example, the WHO estimates 15 million cases worldwide of Hepatitis D, which spreads only amongst carriers of Hepatitis B and greatly worsens their prognosis [25]. There is a need, therefore, to develop a robust understanding of the conditions under which outbreaks of secondary infections are possible. Coevolving infections have been studied previously in the case of symbiotic/antagonistic relationships where infections mutually affect fitness [7,8,11], however, relatively little is known in the case that one infection has a strict obligate relationship with another.
In 2013, Court, Blythe and Allen [6] introduced a model of hierarchical infection referred to as the stacked contact process. Their model concerns the fate of a population of coevolving hosts, spreading as a contact process on a lattice, and parasites, spreading as a contact process restricted to sites currently occupied by hosts. In epidemiological language, the contact processes of [6] correspond to coupled Susceptible-Infective-Susceptible (SIS) epidemics; empty lattice sites are interpreted as susceptible individuals, who may be infected by the primary (host) and then secondary (parasite) infections. Simulations of this model system revealed a surprising feature: the success of the parasites depends non-monotonically on the turnover rate of the host population. Specifically, for the parasite to succeed, it is necessary for the dynamics of the host population to be neither too fast, nor too slow. Later in [19], Lanchier and Zhang rigorously established the main features of the phase diagram for the stacked contact process.
At around the same time, Newman and Ferrario [21] independently proposed a related model in the context of epidemic dynamics in social contact networks. They considered a pair of Susceptible-Infective-Recovered (SIR) epidemics with a strictly obligate relationship such that the secondary infection is only transmitted amongst those who have recovered from the primary. In this formulation, the dynamics of the two diseases are completely separated in time, allowing for analytical treatment of the model using "cavity method" techniques which have been quite successful in the study of epidemics on networks (see, e.g. [20,23]). The introduction of network structure to the population in [21] has the advantage of improving the relevance of the model for human epidemic dynamics, however, by separating the dynamics of the two diseases this model cannot display the curious interaction between infection timescales observed in [6,19].
In this paper we study the dynamics of coevolving SIR superinfections in sparse contact networks. We consider a population of individuals occupying the vertices of an Erdős-Rényi (ER) random graph with mean degree c. A primary infection spreads through the population with infective individuals passing the disease on to their neighbours with rate β 1 , and recovering from the disease with rate ρ 1 . Individuals who are carrying a live primary infection may also play host to a secondary infection, which spreads and recovers with rates β 2 , ρ 2 respectively. See Fig. 1 for an illustration of the possible state transitions. As in [21], our secondary infection is restricted to spread on the subgraph of hosts infected with the primary, however, differently from that paper we consider the more complex case in which this subgraph is evolving in time due to the recovery of primary infections.
As well as arguably improving the realism of the model, moving from lattice to network topologies allows us access to a rigorous branching process approximation-an approach that has previously enjoyed success in approximating SIR-type models in large populations, as seen for instance in [4,24]. By coupling the dynamics of the secondary infection to those of a multi-type branching process, we will be able to characterise the phase diagram of the system. Note that, in the context of large finite networks, when we discuss survival of the Fig. 1 Possible events and their rates in the network superinfection model. Circles represent nodes in the network, with the state of the primary (resp. secondary) infection shown by the colour of the lower-left (resp. top-right) sector; light denotes susceptible, midtone denotes infective, dark denotes recovered Infection Recovery Fig. 2 Phase diagram of the superinfection network model for fixed β 1 /ρ 1 = β 2 /ρ 2 , shown as a function of the relative timescale ϕ = β 1 /β 2 of the infections and the connectivity c of the network. The secondary infection survives with positive probability only in a convex region whose boundary is characterised in our Theorem 1. In this log-log plot, the asymptotic slope of the boundary is − 1 for small ϕ, and 1 for large ϕ, as implied by the scaling laws in (2) infection we mean an asymptotically positive proportion of vertices become infected at some point in time.
The success of the primary infection is controlled by the connectivity (mean node degree) c of the network, and the ratio of the infection and recovery rates α := β 1 /ρ 1 . This parameter is well understood as the basis of the classical single infection process; for fixed α there exists a critical value of c above which the infection survives with positive probability and at or below which we have certain extinction, see e.g. [13]. Note that simultaneously adjusting β 1 and ρ 1 by a multiplicative factor will change the timescale of the disease dynamics, but will not alter the probability of survival since α is unchanged.
Inspired by the results of [6,19], we are interested here in the behaviour possible when the primary and secondary diseases are similarly virulent, but may differ in their in timescales. To this end, we will mainly concentrate on the case that β 2 /ρ 2 = α also. We have made this choice only for simplicity of presentation; the more general case is in fact covered by Lemma 2, and the results are not qualitatively different in other cases. Three parameters then describe success of secondary infection: the connectivity of the underlying graph, c; the ratio between rates of spread and recovery, α; and, crucially, the relative timescales of the two infections, ϕ := β 1 /β 2 . If ϕ 1 then the dynamics of the primary infection are very much faster than those of the secondary; if ϕ 1 then they are much slower. In order for the secondary infection to survive it is perhaps intuitive that it must progress at a rate fast enough compared to the primary infection, else the primary infection will have itself recovered and subsequently ended the secondary infection before it has a chance to spread. Perhaps more surprisingly however we shall also show that the secondary infection should not act too quickly as this too compromises survival potential. Our characterisation of the survival of the secondary infection is illustrated in Fig. 2 and summarised by our main result: Theorem 1 For all α, ϕ > 0 there exists a critical connectivity c such that, in the limit of large network size, for c < c the secondary infection dies out with probability one, and for c > c it survives with positive probability.
Furthermore, the critical connectivity c is found to be the smallest positive solution of the implicit equation where 0 F 1 (a; z) = k z k k!(a) k is a hypergeometric function. In particular, for large and small ϕ we have the scaling behaviour Here we have made use of "big theta" notation, defined as follows: f (x) = (g(x)) as x → ∞ (resp. x → 0) if there exist positive constants L and U and X such that ∀x > X (resp. x < X ) we have The remainder of the article is organised as follows: in the next section we map the early dynamics of our network superinfection model to a certain multitype branching process; in Sect. 3 we compute the long-time behaviour of this process and thus give the proof of Theorem 1; Sect. 4 is for discussion, including illustrative numerical results.

Primary Infection
We begin by recapping the standard branching process approximation to the dynamics of an infection spreading on an Erdős-Rényi random graph [4,17]. Heuristically, the method relies on the fact that for fixed connectivity, short cycles become asymptotically rare in the limit of large graphs, meaning that during the crucial early dynamics of the infection, each susceptible node may have at most one infective neighbour. Let us consider the infection spread as generational; the nth generation being the individuals at graph distance n from the seed vertex that gain the infection at any point in time. In this way the primary infection is modelled as a simple Galton Watson process described by the quantity Z n , giving the number of individuals in the n th generation. The offspring distribution describes the probability p i for an individual to pass the infection on to i others in the next generation. If the offspring distribution has mean μ, the expected number of infected individuals at distance n from the seed is then given by EZ n = μ n . If μ ≤ 1 then the branching process will almost surely go extinct after finitely many generations; if μ > 1 then it may survive; survival of the branching model being simply characterised by the size of the nth generation being non zero for all n.
In the network SIR model, the number of offspring is equated with the number of neighbours (other than the single infected 'parent') that an infective node succeeds in transmitting the disease to before it recovers. There are several sources of randomness: the number of neighbours to potentially infect, the recovery time, and times of infection. We note that whilst the fates of the neighbours of an infected node are not independent (they are jointly exposed to the random time to recovery of the parent) the mean of the offspring distribution can be found simply by multiplying the probability α/(1 + α) to infect any given neighbour before recovery, with the expected number of neighbours to infect, c. From standard branching process theory, we thus deduce that in the limit of large networks the primary infection will have a non-zero chance of survival if and only if μ > 1, that is, if c > 1 + 1/α.
For finite graphs, the coupling between the random graph and branching process model is of course only local. Suppose in a population of size N , in generation n we have m infected individuals, so Z n = m in the branching model. We then have errors coming from the fact that each infective may only be connected to at most N − m susceptibles (not constant for each generation) as well as the fact that each of these may not be unique (and so children in the subsequent generation of the branching process may not be unique). However when m = o √ n the random graph may be coupled to the branching model with high probability; for a proof of this see [10].

Secondary Infection
For the primary infection, to determine if the probability of survival is positive only requires knowledge of two quantities: the expected number of susceptible neighbours an individual has, and the chance any one of those will gain the infection. The difficulty with modelling the secondary infection is that the first of these is dynamic, since the subgraph composed of individuals currently carrying the primary infection changes with time. We account for this additional complexity by introducing a type parameter t, which specifies the time elapsed between the primary and secondary infections. Specifically, if an individual acquires the primary infection at time t 1 and the secondary at time t 2 , then they are said to have type It is clear that at least this much information is required to predict the potential of an individual to transmit the secondary infection to new hosts; for example, the larger t, the more likely an individual is to pass on the primary infection long before it passes on the secondary, by which time the primary infection in the new host may have recovered. We will see in Sect. 3.1 that in fact knowledge of t is enough to completely characterise the distribution of the number and timing of new secondary infections arising from an individual. The progress of the secondary infection is then mapped to that of a multi-type branching process with type space T = [0, ∞).
Where previously survival was predicted by just the mean number of offspring, now the picture is more complicated, and we are required to compute the intensity of production of all types of offspring resulting from all types of parents. This information is captured in the kernel μ(t |t), which is defined by the property that the expected number offspring with types in the interval [a, b] coming from a parent of type t is given by the integral of μ(t |t) over t ∈ [a, b]. This kernel defines a linear operator with the action In words, M[ψ] describes the expected size and composition of the population of offspring arising from a population of parents with types given by ψ.
We say that a kernel μ defined over an interval I is: strictly positive if ∀t, t ∈ I we have μ(t |t) > 0; uniformly positive if ∃ ε > 0 such that ∀t, t ∈ I, μ(t |t) > ε; integrable if μ(t |t)dt dt < ∞. We assume that M can be defined as a linear operator M : C b (T) → C b (T) over the space of continuous bounded functions on the compact intervalT = [0, ∞] equipped with the supremum norm, and in particular that μ has vanishing mass as t goes to infinity. One then has the following general result: Note that inf t,t ∈[0,T ] μ(t |t) > 0 and so the kernel is strictly positive and we refer to Note that operators M (T ) andM (T ) share eigenvalues so we may equivalently consider the top eigenvalueλ (T ) ofM (T ) . Since μ is continuous and integrable, for all ε > 0 there exists T such that where · · · is the operator norm induced by the infinity norm on C b (T).
We have already observed that the principal eigenvalue λ of M can be separated from the rest of the spectrum by a closed curve. Hence, by Kato and hence λ (T ) > 1 and {Z (T ) n } survives with positive probability. The untrimmed process satisfies Z n ≥ Z (T ) n and hence also survives with positive probability.
To prove our main result about the survival of the secondary infection, we must explicitly identify the operator M, analyse its spectrum, and compute the scaling behaviour when the timescales of the infections are well separated.

Production Kernel
The form of the kernel μ(t |t) may be found by considering when a type t parent will have a type t offspring. For this to happen, the parent must pass on the primary infection at some time s (measured from the moment they first acquired it), and then pass on the secondary infection at time s + t . The primary and secondary infections in the parent, and the primary infection in the child, must all survive long enough for this process to complete. We find it useful to break the calculation into two cases, depending on whether the primary infection is transmitted before or after the parent acquires the secondary; that is, depending on the order of s and t. The case s < t is illustrated in Fig. 3(i). To achieve a type t offspring in this case: the transmission time s > 0 of the primary must occur before t but after t − t (which may be negative); the secondary must be transmitted s + t − t time units after it was acquired in the parent at time t; the primary infection in the parent must not recover in the time between s and t; and none of the three active infections may recover in the window of time between t and s + t . Putting these contributions together, we reach where (· · · ) + denotes the positive part, and the prefactor of c comes from the expected number of neighbours to which the infection may be transmitted.
Similarly, the case s ≥ t is illustrated Fig. 3(ii). Here transmission of the primary may occur any time after t, with the secondary being transmitted t time units later. Both infections in the parent must survive until time s, after which all three infections must survive for at least t time units. The resulting expression is Combining the two cases and evaluating the integral gives We are now ready to state our result about the spectrum of the production operator resulting from this kernel.
is a hypergeometric function.
Proof From part 1 of Lemma 1, to determine that λ is the top eigenvalue of M, it is sufficient to exhibit a non-negative function ψ such that λψ = Mψ. We begin a search for such a function by considering the successive action of M starting from the initial state ψ 0 (t) = δ 0 (t), corresponding to a single seed infected individual who acquires the primary and secondary infections at the same instant. Defining the series we observe that each iterate ψ n is a member of a family, Ψ , of functions that can be written as a certain positive sum of exponentials: We look for an eigenfunction of M that lies in Ψ . The eigenvalue equation λψ = M[ψ] is thus reduced to a statement about the coefficients {a k }. Specifically, we find where b k = β 1 β 2 (β 1 + (k + 1)ρ 1 + ρ 2 ) Equating coefficients in (12) determines where the {a k } are found to satisfy This recursive equation specifies a solution up to a multiplicative constant: where (· · · ) k denotes the Pochhammer symbol. Combining this result with (13), yields the implicit equation (9) for λ given in the statement.

Bounds on the Ratio of Hypergeometric Functions
Recall that the survival of the primary infection is dependent only on its birth-death ratio α and the connectivity of the underlying graph c, while the secondary infection additionally depends on its relative speed when compared primary, ϕ := β 1 /β 2 . As per the discussion in Sect. 2.2, we specialise to the case that β 1 /ρ 1 = β 2 /ρ 2 = α. Then the implicit eigenvalue equation (9) can be rewritten in terms of the parameters α and ϕ to give where γ = (1 + ϕ)(1 + 1/α) and denotes the hypergeometric ratio Our strategy to prove the scaling relations claimed in Theorem 1, will be to replace this function by suitably simple upper and lower bounds with the same asymptotic behaviour. Fortunately, there is a substantial literature on topic that we may draw on.

Lemma 3
For a > 0 write j a for the smallest positive root of J a , the Bessel function of the first kind. Then a(a + 2) < j 2 a < 4(a + 1)(a + 2), and for all z ∈ (0, j a ) we have Proof Ismail and Muldoon [15] list many different bounds on j a , including those in (18) coming from formulas (6.7) and (6.22) in that article. For the second part, it is well-known [1] that the Bessel functions of the first kind may be expressed as This function has previously been studied by Ifantis and Siafarikas [14], who proved various inequalities including their formulas (1.2) and (2.17) which imply the lower and upper bounds of (19).

Proof of Theorem 1
Proof As argued previously, in the limit of large Erdős-Rényi random graphs with mean degree c, the survival probability of the secondary infection coincides with that of a multitype branching process {Z n } with production kernel given by Eq. (8). From Lemma 2 and Theorem 1 we establish that Z n has non-zero probability to survive indefinitely if and only if λ > 1, where λ is the largest real number satisfying Noticing that λ appears only in ratio with c, it follows that the condition for the possibility of survival may be rewritten in terms of the critical connectivity c such that for c > c we have λ > 1. Rearranging Eq. (21) we straightforwardly find that c is the smallest positive solution to which is precisely Eq. (1), as required.
To quantify the scaling behaviour of c for large and small ϕ, we recall the definition of "big theta" notation: f (x) = (g(x)) as x → ∞ (resp. x → 0) if there exist positive constants L and U and X such that ∀x > X (resp. x < X ) we have Two sufficient conditions are easy to check: We will use the bounds in Lemma 3 to exhibit functions with appropriate finite limits that sandwich c . Specifically, recalling γ = (1 + ϕ)(1 + 1/α), let First we check the upper bound. From (22) and the lower bound of unity in Eq. (19) of Lemma 3, we have that For the lower bound, we note that the upper bound on given in Lemma 3 implies a lower bound on c as the smallest positive l satisfying the equation In fact there is only one solution: The lower bound l(ϕ) < c given in (23) follows immediately from this and the lower bound on j 2 γ given in Eq. (18) of Lemma 3. It remains to check that the upper and lower bounds both have the desired scaling in large and small ϕ. We begin with u(ϕ), which has easily determined limits both of which are finite and positive, implying u(ϕ) = (ϕ) for large ϕ and u(ϕ) = (1/ϕ) for small ϕ. For the lower bound we use these results to obtain It follows from the defintion of l(ϕ) and finiteness of these limits that l(ϕ) has the same scaling form as u(ϕ) for both large and small arguments. Since u and l sandwich c , the desired scaling is confirmed.

Discussion
Theorem 1 provides an exact but implicit formula for the region in which survival of the secondary infection is possible (in the limit of infinitely large graphs), and establishes the scaling behaviour of the boundary of this region for large and small values of the parameter ϕ = β 1 /β 2 which controls the relative timescales of the two infections. Knowledge of this scaling behaviour is enough to prove that, for fixed α and c, the survival of the secondary infection is confined to a bounded region of ϕ values-this is the reentrant phase transition of our title. Figures 4 and 5 show the results of numerical simulations of both the branching process and the network model to illustrate this phenomenon. It is interesting to note that the simulations of the network process and the limiting branching process are not in perfect agreement. Viewing the mean outbreak size over 1000 runs of the model we see in Fig. 4 that, while we have agreement with the branching process for small values of ϕ, large outbreaks still seem to be possible beyond the point predicted by the branching process. Moreover, considering the individual simulation results it seems that this unexpected tail is comprised of a few very large outbreaks; while outbreaks of any size are rare for large values of ϕ, when they do happen they reach most of the graph. By considering the infection spread in a closed connected community we start to encounter finite size effects. Recall that the branching approximation is only valid when the number of infected is relatively small compared to the size of the graph. As the outbreak becomes large the approximation breaks down, a problem exacerbated by the two levels of infection we study. Furthermore in a more highly connected environment we may have the existence of transmission routes for the secondary infection to primary infected cousins as well as direct descendants allowing opportunity for the secondary infection to progress before direct primary progression. Similar finite size scaling effects have been observed in other coevolving infection models; see [7] for example. Comparing the average outbreak size with individual realisations demonstrates an interesting choice of risk vs reward in the strategy of a secondary infection, due to the different locations of the maxima of the curves shown in the top panel of Fig. 4. The values of ϕ for which outbreaks are most likely to occur (blue curve) are in the lower end of the survival window, corresponding to smaller total outbreak sizes (black stars). Conversely, larger values of ϕ have potential for much larger outbreaks, but come with a higher risk of rapid extinction. Looking at this another way, in nature we should expect survival probability to be a strongly selected characteristic, and hence to find that the majority of secondary infections reach only a minority of primary hosts.
The work presented here could easily be extended to a host of other random graph models, for example by building on techniques of [2,5,18]. It may also be interesting to explore the application of the model (or variants) to other areas, including: the successive invasions of different species necessary to rebuild a diverse ecosystem in a damaged habitat; the evolution of hyperparasitism (that is, parasites that live on other parasites); radicalisation, and the incremental spread of increasingly extreme political views through social media.