## Abstract

The percolation paradigm is widely used in spatially explicit epidemic models where disease spreads between neighbouring hosts. It has been successful in identifying epidemic thresholds for invasion, separating non-invasive regimes, where the disease never invades the system, from invasive regimes where the probability of invasion is positive. However, its power is mainly limited to homogeneous systems. When heterogeneity (environmental stochasticity) is introduced, the value of the epidemic threshold is, in general, not predictable without numerical simulations. Here, we analyse the role of heterogeneity in a stochastic susceptible–infected–removed epidemic model on a two-dimensional lattice. In the homogeneous case, equivalent to bond percolation, the probability of invasion is controlled by a single parameter, the transmissibility of the pathogen between neighbouring hosts. In the heterogeneous model, the transmissibility becomes a random variable drawn from a probability distribution. We investigate how heterogeneity in transmissibility influences the value of the invasion threshold, and find that the resilience of the system to invasion can be suitably described by two control parameters, the mean and variance of the transmissibility. We analyse a two-dimensional phase diagram, where the threshold is represented by a phase boundary separating an invasive regime in the high-mean, low-variance region from a non-invasive regime in the low-mean, high-variance region of the parameter space. We thus show that the percolation paradigm can be extended to the heterogeneous case. Our results have practical implications for the analysis of disease control strategies in realistic heterogeneous epidemic systems.

## 1. Introduction

One of the main challenges in epidemiology is to explain the huge variability observed in the outcome of epidemics. Such variability affects the invasion and persistence of disease in a host population, the rate of disease progress, the final fraction of the population affected, and even whether or not the disease spreads at all (Gilligan 2002). Much of the observed variation is owing to the inherent stochasticity of the epidemic process, acting at different spatial scales. At the small scale, *demographic* stochasticity intervenes in the interaction between host and pathogen: pathogen transmission between two individuals is in general a random event. Epidemic spread also depends upon epidemiological parameters that vary across the population owing to stochastic *environmental* factors, whereby hosts can differ in properties such as susceptibility or infectivity. When the spatial structure of the epidemic is explicit, the interplay at different scales between environmental stochasticity (or *heterogeneity*) and demographic stochasticity (henceforth, simply stochasticity) contributes to the variability of the outcome reflected in the generation of spatio-temporal patterns of infected hosts (Durrett & Levin 1994; Shaw 1994; Zadoks & van den Bosch 1994). Among botanical epidemics in particular, such spatio-temporal variability is evident not only between different epidemics of the same host–pathogen system, but also within the same epidemic; for example through the occurrence of patches of disease within a single field, with some areas exhibiting high and others low incidences (Gibson *et al.* 1999).

For a broad class of host–pathogen systems, transmission only occurs between neighbouring hosts. This class includes many plant pathogens, such as soil-borne pathogens transmitted between neighbouring plants (Bailey *et al.* 2000; Otten *et al.* 2004), or, at a larger scale, pathogens spreading within a mosaic of neighbouring susceptible fields or even farms (Truscott & Gilligan 2001; Stacey *et al.* 2004). The class also comprises animal pathogens that spread in host populations living in a fixed habitat, as has been shown, e.g. for plague in populations of gerbils inhabiting a network of connected burrows (Davis *et al.* 2008). Finally, at the microscopic scale, some diseases are known to spread among contiguous cells (cf. neural ganglion cells in Pérez-Reche *et al.* 2010). Nearest-neighbour interactions lead naturally to consideration of *percolation* (Stauffer & Aharony 1991). Percolation theory provides tools and concepts to scale up from the small-scale connectivity properties of a system (e.g. whether or not two neighbours are connected) to its large-scale properties, such as the size of clusters of connected sites. In epidemiology, percolation has been used to study the ‘infection path’ followed by the pathogen through the system (Mollison 1977; Grassberger 1983), and the connectivity of habitable sites in the landscape (Otten *et al.* 2004). The most relevant result is the existence of a *spatial threshold* for epidemic invasion. The threshold separates a non-invasive, safe regime (where the pathogen can never invade the population) from an invasive regime (where the pathogen has a non-zero probability to invade). However, including the effect of host heterogeneity in the *percolation paradigm* has proved to be not straightforward. With a few notable exceptions (Kuulasmaa 1982; Cox & Durrett 1988; Sander *et al.* 2002) and recent interest in the field of complex networks (Kenah & Robins 2007; Miller 2007, 2008; Meester & Trapman 2008), the problem has been mostly neglected.

In this paper, we propose a significant advance towards the solution of the problem of incorporating host heterogeneity in spatial systems. We study the effects of heterogeneity on a stochastic SIR (susceptible–infected–removed) model on a lattice (Murray 2002). The model is chosen because an exact mapping onto bond percolation holds in the homogeneous case (Grassberger 1983), where an epidemiological parameter, the *transmissibility* *ψ* (the probability to transmit infection along a bond), is equivalent to the bond probability in percolation. The value of *ψ* is thus the only parameter controlling the resilience of the system to invasion, and, in particular, the invasion threshold *ψ*_{c} coincides with the bond-percolation threshold *p*_{c}^{bond}. This also implies that in a few relevant cases the value of the invasion threshold is known analytically (Stauffer & Aharony 1991). Thus, the model provides an ideal starting point to study the effects of inherent heterogeneity in spatially explicit systems. Our aim is to understand how heterogeneity changes the invasion threshold of the system, and whether the threshold can be still predicted without the need for numerical simulations or not. One of the main challenges in the study of a heterogeneous system is to find if it can be described by means of an equivalent homogeneous mean-field system characterized by effective parameters. It is known (Cox & Durrett 1988) that the effect of heterogeneity depends crucially on the presence of *correlations* in *ψ* along different bonds. In the uncorrelated case, the system is equivalent to an effective homogeneous system with *ψ* = 〈*ψ*〉, and the threshold is given by 〈*ψ*〉 = *p*_{c}^{bond} (Sander *et al.* 2002). In the correlated case, no such equivalence holds, and the value of the invasion threshold cannot be inferred from 〈*ψ*〉. The threshold is known to shift from 〈*ψ*〉 = *p*_{c}^{bond} (in homogeneous and heterogeneous uncorrelated systems) to 〈*ψ*〉 = *p*_{c}^{site} (site percolation threshold), when the variance of the transmissibility *σ*_{ψ}^{2} in correlated systems is maximal (Kuulasmaa 1982; Kuulasmaa & Zachary 1984; Cox & Durrett 1988). This result points to *σ*_{ψ}^{2} as a good measure of heterogeneity in correlated systems, but it poses several questions. Is the invasion threshold an increasing function of the variance of *ψ*, and in general, how much information about the threshold is provided by *σ*_{ψ}^{2}? More specifically, do we need to know the whole distribution of *ψ* to estimate the threshold?

These questions are systematically addressed in the paper. Our main result is the existence of a phase diagram in the two-dimensional space (〈*ψ*〉, *σ*_{ψ}^{2}), with a (lattice-dependent) phase boundary separating invasive and non-invasive regimes. We also show that for all the correlated systems studied, the position of the threshold is determined with excellent approximation by the two parameters 〈*ψ*〉 and *σ*_{ψ}^{2}. Higher (and in general smaller) moments of the distribution *h*(*ψ*) shift the threshold value for 〈*ψ*〉 within the interval [*p*_{c}^{bond}, *p*_{c}^{site}] only by a small amount. For the square lattice, for example, the shift is at most 3 per cent of the interval length, which can be considered negligible for most applications. Our phase diagram extends the predictive power of the percolation paradigm for SIR epidemics to heterogeneous correlated systems. These results also have relevant implications for the implementation of control: in this respect, a heterogeneous correlated system is qualitatively different from its homogeneous counterpart, because the parameter space in which to search for an optimal control strategy is now two-dimensional. While the task of finding an optimal control strategy becomes non-trivial, we nevertheless show that the problem can be solved. This can open up new possibilities for the implementation of control in spatial epidemic systems.

## 2. Model

We give here a brief description of the model and the problems involved; a more detailed version is contained in the first section of the electronic supplementary material. We consider a stochastic SIR process on a square lattice (for concreteness) with *N* sites, where an epidemic spreads by the nearest-neighbour transmission of a pathogen, starting from a single infected site (primary infection). Upon infection, a host (donor) *i* is infectious (I) for a time *τ*_{i} (infectious period) before recovering or being removed (R). During this time, it transmits the pathogen to its susceptible neighbours (recipients) *j* according to a Poisson process with rate *β*_{ij} (infection rate). The probability that the pathogen is transmitted from *i* to *j* within the time *τ*_{i} is the *transmissibility*:
2.1

The system is *heterogeneous* when *β*_{ij} and *τ*_{i} are random variables. We consider them to be independent and identically distributed (iid) random variables, drawn from the distributions *f*(*β*_{ij}) and *g*(*τ*_{i}), respectively. In this case, the *ψ*_{ij} are also random (but not necessarily independent) variables, with a probability distribution *h*(*ψ*_{ij}) that can be calculated from *f*(*β*_{ij}) and *g*(*τ*_{i}).

In a finite system, the process stops after a finite time with only S and R hosts left. The crucial question we address in this paper is how the distributions of the infection rates *f*(*β*_{ij}) and infectious periods *g*(*τ*_{i}) affect the probability that a *major outbreak* or *invasion* occurs (Diekmann & Heesterbeek 2000), when the final set of R sites spans the entire lattice. We also consider whether the resilience of the system to epidemic invasion can be described by a few effective control parameters. The evaluation of the probability of invasion for the process *P*_{inv} (*f*(*β*_{ij}),*g*(*τ*_{i})) and the identification of the relevant parameters controlling invasion are the main objectives of this paper.

If the system is *homogeneous* (constant parameters: *τ*_{i} = *τ*, *β*_{ij} = *β*), the identification of the relevant parameters is straightforward. The final state of the SIR process can be mapped onto a bond percolation problem (Grassberger 1983; Stauffer & Aharony 1991), where the probability that a bond is ‘open’ is equal to *ψ* (see the electronic supplementary material for details). This leads to the following definition of a critical (threshold) value for the transmissibility,
2.2
where *p*_{c}^{bond} is the bond-percolation threshold (*p*_{c}^{bond} = ½ for the square lattice; Stauffer & Aharony 1991). In an infinite system, the critical value separates an *invasive* region in the parameter space (*P*_{inv}(*β*, *τ* ) > 0 for *ψ*_{c}> *p*_{c}^{bond}) from a *non-invasive*, safe region (*P*_{inv}(*β*, *τ*) = 0 for *ψ*_{c}< *p*_{c}^{bond}). The same holds for a *bond-heterogeneous* process (*τ*_{i} = *τ* constant, *β*_{ij} random), for which the *ψ*_{ij} are iid random variables. Since different bonds are open or closed independently, the process is still equivalent to homogeneous bond percolation, and equation (2.2) still holds, once *ψ* is replaced by 〈*ψ*〉_{β} (the average of *ψ*_{ij} from equation (2.1) with respect to *f*(*β*), Sander *et al.* 2002). Hence, there still exists a single control parameter 〈*ψ*〉_{β}, independent of the moments of the distribution *f*(*β*) other than the first. The solution is more challenging in the case of a *site-heterogeneous* process (*β*_{ij} = *β* constant, *τ*_{i} random). The *ψ*_{ij} = 1 − e^{−βτi} = *ψ*_{i} here are host-dependent and are not iid random variables). The equivalent percolation problem is *correlated*, because bonds outgoing from the same site are not open or closed independently (see the electronic supplementary material for proof and an explicit example). As a result, the relationship between the quantity 〈*ψ*〉_{τ} (the average of *ψ*_{ij} with respect to *g*(*τ*)) and the resilience of the system is less straightforward.

Indeed, the average transmissibility at the invasion threshold can take any value in the range
2.3
with bounds corresponding to the bond- and site-percolation thresholds *p*_{c}^{bond} and *p*_{c}^{site}. The system is always in the non-invasive regime for 〈*ψ*〉_{τ} < *p*_{c}^{bond}, and always in the invasive regime for 〈*ψ*〉_{τ} > *p*_{c}^{site} (whether the system is in the invasive regime or not for 〈*ψ*〉_{τ} = *p*_{c}^{bond} and 〈*ψ*〉_{τ} = *p*_{c}^{site} depends on the topology of the lattice and is still an open problem: e.g. Grimmett 1999); we do not discuss these cases here). However, in the middle interval *p*_{c}^{bond} < 〈*ψ*〉_{τ} < *p*_{c}^{site} the system can be in either regime, depending on the particular distribution *g*(*τ*). Thus, the average transmissibility 〈*ψ*〉_{τ} alone is no longer a good control parameter to describe the resilience of the system to invasion. The width of the interval changes with the topology of the system: in general, *p*_{c}^{bond} = *p*_{c}^{site} on a tree graph, while *p*_{c}^{bond}< *p*_{c}^{site} for lattices with dimension *d* ≥ 2 (Grimmett 1999). For the square and triangular lattices, the interval (*p*_{c}^{bond}, *p*_{c}^{site}) is equal to (½, 0.593 … ) and (0.347 … , ½), respectively.

The challenge is to find alternative (or additional) control parameters that allow the prediction of the value of the epidemic threshold in a site-heterogeneous system. A hint comes from the fact that the critical value *p*_{c}^{bond} is obtained when the *ψ*_{ij} are uncorrelated, while the value *p*_{c}^{site} is obtained when the *ψ*_{ij} outgoing from the same site are *maximally* correlated. In the former case, the *ψ*_{ij} are drawn from a distribution with null variance, in the latter case from a distribution with maximal variance (see the electronic supplementary material for details). This leads to the hypothesis that the value of the threshold increases with the variance of *ψ*. This hypothesis is investigated below.

In what follows, we start from the analysis of the influence of correlated heterogeneity on the invasion threshold in a model system (the *mixed* case, where both the correlated and uncorrelated kinds of heterogeneity are present, is considered in the electronic supplementary material). The model is a two-host system, where heterogeneity is owing to the existence of two classes of susceptible hosts, randomly arranged on a lattice. It is used here to test how the invasion threshold is influenced by heterogeneity, and to demonstrate to what extent a description based on a homogeneous approximation using the average transmissibility fails. The variance of the transmissibility is shown to play an important role in this case. We also demonstrate that the resilience of the system to epidemic invasion is increased when the spatial pattern of the two classes of hosts is regular. The study of the two-host system leads us to the main result of the paper that links the value of the invasion threshold to the variability of a generic heterogeneous system. We proceed in two steps. In the first step, we investigate the dependence of the threshold on *σ*_{ψ}^{2} for specific families of bimodal distributions for *ψ* (corresponding to different families of two-host systems), and show that for a given 〈*ψ*〉 there is a transition from the invasive (low-variance) to the non-invasive (high-variance) regime for a critical value *σ*_{ψ}^{2} = *σ*_{ψ,c}^{2}. This demonstrates the existence of a *phase diagram* in the plane (〈*ψ*〉, *σ*_{ψ}^{2}), where a phase boundary separates the two regimes. The second step is to show that the influence of higher-order moments of the distribution for *ψ* on the location of the phase boundary is negligible. We find that the position of the phase boundary for different families of distributions (both bimodal and other than bimodal, with common values of 〈*ψ*〉 and *σ*_{ψ}^{2} but differing in all the higher moments) is approximately the same *for all distributions* tested. This means that the pair of control parameters (〈*ψ*〉, *σ*_{ψ}^{2}) can determine if a system is in the invasive regime or not, up to small corrections owing to higher moments.

## 3. A testing ground: the two-host system

In this model system, hosts belonging to two different classes, A and B, are arranged on the square lattice, occupying a fraction *ρ*_{A} and *ρ*_{B} = 1 − *ρ*_{A} of the lattice sites, respectively. This topology is analogous, for example, to agricultural systems where different plant species (or varieties) coexist in the same field (Garrett & Mundt 1999; Cook *et al.* 2007), or at a larger scale, where some fields are treated with a pesticide and some others are not. Two different scenarios are considered below: (i) site heterogeneity with A and B hosts randomly placed on the lattice sites, and (ii) periodic arrangements of hosts. A third scenario, mixed heterogeneity, is discussed in the electronic supplementary material. In the case of site heterogeneity, each class is characterized by its own transmissibility, *ψ*_{A} and *ψ*_{B}: this occurs when two different infectious periods, *τ*_{A} and *τ*_{B}, are associated with different hosts, and the rate *β* is host-independent. The distribution, mean and variance of *ψ* are
3.1a
3.1b
3.1c

Here and in the rest of the paper, *δ*(*x*) is the Dirac delta-function with the property ∫_{−∞}^{+∞} *f*(*x*)*δ*(*x* − *x*_{0})d*x* = *f*(*x*_{0}). The (invasive or non-invasive) behaviour of the system depends on the values of the parameters *ρ*_{A}, *ψ*_{A} and *ψ*_{B}. In figure 1, the phase diagram is shown in the parameter plane (*ψ*_{A}, *ψ*_{B}) for different values of *ρ*_{A}. The phase boundary separating the invasive and non-invasive regimes lies between the straight lines corresponding to the lower and upper bounds, defined by the equations 〈*ψ*〉 = *p*_{c}^{bond} and 〈*ψ*〉 = *p*_{c}^{site}, with 〈*ψ*〉 given by equation (3.1*b*). The point *ψ*_{A} = *ψ*_{B} corresponds to a homogeneous system for which the phase boundary coincides with the lower bound. As heterogeneity is increased (*ψ*_{A} ≠ *ψ*_{B}), the boundary approaches the upper bound. The degree of heterogeneity of the system can be described quantitatively by the variance of *h*(*ψ*), given by equation (3.1*c*). It increases in the direction indicated by the arrows (see figure 1). The threshold for invasion increases monotonically with the variance. This behaviour is studied in more detail in the following section. It is worth noting that in the model considered here (and in this paper in general) the transmissibility between donor and recipient depends on the donor only: recent work (Miller 2007, 2008; Meester & Trapman 2008) has considered the more general case where the transmissibility depends both on the donor and the recipient (i.e. both host infectivity and host susceptibility are taken into account).

The results discussed so far, and the theorems giving bounds for heterogeneous SIR processes (Kuulasmaa 1982; Cox & Durrett 1988), hold for systems where the transmissibility of any host is a random variable which does not depend on the host position in space. As soon as this condition is dropped, the bounds do not hold any longer: an example is a two-host system where A and B hosts follow a regular spatial pattern. Such a system can be compared with its random counterpart with the same values of *ψ*_{A}, *ψ*_{B} and *ρ*_{A}. Comparison of the invasive regions of two periodic systems, one in the form ABAB, the other AABAAB (see figure 2), show that the periodic systems are more resilient to invasion.

## 4. Dependence of the epidemic threshold on the variance of *ψ*

When the distribution *h*(*ψ*) contains bond heterogeneity only, the resilience of a heterogeneous system to epidemic invasion is completely determined by the average transmissibility 〈*ψ*〉. If site heterogeneity is also present (e.g. owing to a distribution of infectious periods), the resilience of the system cannot be determined by 〈*ψ*〉 alone. In principle, knowledge of the whole distribution *h*(*ψ*) is needed. Systems with different distributions characterized by the same 〈*ψ*〉 can be in different regimes (invasive or non-invasive) depending upon moments of higher order, 〈(*ψ* − 〈*ψ*〉)^{k}〉, *k* = 2,3, … Below, we demonstrate that only the second moment plays an important role in systems with site heterogeneity.

As follows from equation (2.3), in the case of site heterogeneity the threshold 〈*ψ*〉_{c} varies between the lower bound *p*_{c}^{bond}, corresponding to zero variance in *ψ*, and the upper bound, *p*_{c}^{site}, corresponding to the maximal variance for a given value of 〈*ψ*〉,
which is achieved for a particular two-host system with *h*(*ψ*) = 〈*ψ*〉*δ*(1 − *ψ*) + (1 − 〈*ψ*〉)*δ*(*ψ*) (cf. equation (3.1*a*) and the electronic supplementary material). Therefore, when the average transmissibility falls within the window *p*_{c}^{bond} ≤ 〈*ψ*〉 ≤ *p*_{c}^{site}, the system can be in different regimes (invasive or non-invasive) depending on the variance of the distribution, 0 ≤ *σ*_{ψ}^{2} ≤ *σ*_{ψ,max}^{2}. We carried out a numerical investigation to find the phase boundary for the SIR process in the (〈*ψ*〉, *σ*_{ψ}^{2}) plane for several different families of distributions where both 〈*ψ*〉 and *σ*_{ψ}^{2} could be adjusted (details about the numerical approach and the distributions used are in the electronic supplementary material). The dependence of the probability of invasion *P*_{inv} on 〈*ψ*〉 and *σ*_{ψ}^{2} for the family of bimodal distributions (equation (3.1*a*)) with *ψ*_{B} = 0 is shown in figure 3. In general, the probability of invasion *P*_{inv} is an increasing function of 〈*ψ*〉 and a decreasing function of *σ*_{ψ}^{2} (see figure 3*a*). The phase transition from the invasive to the non-invasive regime can be achieved by either increasing 〈*ψ*〉 or decreasing *σ*_{ψ}^{2} as control parameters. In the first case (red line in figure 3*a*,*b*), the transition occurs at a value 〈*ψ*〉_{c} such that *p*_{c}^{bond} ≤ 〈*ψ*〉_{c} ≤ *p*_{c}^{site}. In the second case (green line in figure 3*a*,*c*), it occurs for a *σ*_{ψ,c}^{2} such that 0 ≤ *σ*_{ψ,c}^{2} ≤ *σ*_{ψ,max}^{2}. The bound *σ*_{ψ}^{2} = 0 (where *P*_{inv} > 0) corresponds to a homogeneous process equivalent to bond percolation, while *σ*_{ψ}^{2} = *σ*_{ψ,max}^{2} (where *P*_{inv} = 0) corresponds to a maximally correlated process equivalent to site percolation.

The behaviour of the phase boundary between the two regimes (corresponding to the steep transition in the surface in figure 3) can be studied by projecting the surface onto the (〈*ψ*〉, *σ*_{ψ}^{2}) plane (see figure 4*a*). In the phase diagram, the phase boundary partitions the space into a non-invasive (low-〈*ψ*〉, high-*σ*_{ψ}^{2}) and an invasive (high-〈*ψ*〉, low-*σ* _{ψ}^{2}) region. As predicted, the phase boundary goes from the point (*p*_{c}^{bond}, 0) (homogeneous system, bond percolation case) to (*p*_{c}^{site}, *p*_{c}^{site} (1 − *p*_{c}^{site}); maximally heterogeneous system, site percolation case). A system with *p*_{c}^{bond} ≤ 〈*ψ*〉 ≤ *p*_{c}^{site} can be in the invasive or in the non-invasive regime depending on the value of *σ*_{ψ}^{2}.

By changing the average and variance of *ψ* in the family of bimodal distributions with *ψ*_{B} = 0, we are also changing the higher moments that depend on the same parameters and are not shown in figure 4*a*. The position of the phase boundary in figure 4*a* can, in general, also depend on higher moments, and thus on the particular family of distributions chosen. We investigated this effect by comparing phase diagrams for several different families of distributions. Figure 4*b* displays the phase boundaries obtained for some of the families (see the electronic supplementary material for details). The distance between the phase boundaries corresponding to different families turned out to be small for all the examples that we tested. Considering points of the boundaries with the same value of *σ*_{ψ}^{2}, the difference between the corresponding values of 〈*ψ*〉 is always ≤0.03*Δ*_{p}, where *Δ*_{p} = *p*_{c}^{site} − *p*_{c}^{bond} (cf. figure 4*b*). This shows that higher moments have a negligible effect on the position of the phase boundary. This result was confirmed by simulations on regular lattices different from the square (see electronic supplementary material). In the case of the square lattice, the phase boundary is also well approximated by a straight line joining the two points (*p*_{c}^{bond}, 0) and (*p*_{c}^{site}, *p*_{c}^{site} (1 − *p*_{c}^{site})).

We conclude that the resilience of a site-heterogeneous system to epidemic invasion is mainly determined by the first two moments of the distribution for the transmissibility, 〈*ψ*〉 and *σ*_{ψ}^{2}. Changing only the higher moments of the distribution gives corrections to the position of the phase boundary that are negligible for practical applications (such as control, see §5). In other words, by going from a homogeneous to a heterogeneous system, the number of relevant control parameters increases from one to two. In the electronic supplementary material, we give a heuristic explanation of this result based on the zero-function technique introduced by (Kuulasmaa 1982).

Our discussion is focused on the position of the invasion threshold. We do not address here the question of the dependence of *P*_{inv} on heterogeneity when the system is in the invasive regime (see Kuulasmaa 1982; Kenah & Robins 2007; Miller 2007). In the electronic supplementary material, we discuss this dependence in the case of the Bethe lattice, where an exact analytical solution can be found.

## 5. Discussion and conclusions

We have developed a method to analyse the effects of heterogeneity on spread of epidemics in spatial systems. The epidemic spread for an SIR model is defined by transmissibilities between hosts. The transmissibility, in general, depends both on bond and site properties. It is convenient to distinguish between two different types of heterogeneity in transmissibility that can occur in real systems: (i) bond heterogeneity, associated with variability of links between different hosts and (ii) site heterogeneity, determined by variability amongst hosts. The first type, bond heterogeneity, is uncorrelated when the transmissibilities between different pairs of hosts are independent of each other. Site heterogeneity is necessarily *correlated*, because all the transmissibilities from a particular donor host are dependent on the properties of the host. The analysis of SIR epidemics in systems with bond heterogeneity is relatively straightforward. Site heterogeneity is less obvious and, as it follows from our analysis, the two first moments of the transmissibility are crucial for the evaluation of the probability of invasion. One interesting finding is that in systems with correlated heterogeneity it is possible to reduce significantly the probability of invasion by increasing the variance in transmissibility and keeping the mean value unchanged. This feature is a consequence of the behaviour evident in the phase diagram for invasion in the (〈*ψ*〉, *σ*_{ψ}^{2}) space (figure 4). This fact can be used for an effective control of epidemics, whereby both average and variance of the transmissibility can be adjusted to bring the system to the non-invasive regime.

As an example, we analyse its use for the implementation of control strategies in realistic systems with natural variability. With respect to conventional approaches for homogeneously mixed systems, the novelty lies in the fact that we are allowed to use both 〈*ψ*〉 and *σ*_{ψ}^{2} as parameters for the optimization of the control strategy. Consider for example a homogeneous population with constant transmissibility *ψ*_{0} = 〈*ψ*〉 > *p*_{c}^{bond}, hence in the invasive regime (for example, a mosaic of crops, each representing a host unit). In the phase diagram (figure 5), the system is represented by a point on the *σ*_{ψ}^{2} = 0 axis on the right side of the phase boundary (approximated by a straight line). Suppose that a control agent (a protectant) is applied before the arrival of infection in order to bring the system into the non-invasive regime. The effect of the agent on a host is to decrease the transmissibility from *ψ*_{0} to *ε**ψ*_{0} ≤ *ψ*_{0}, where *ε* is a parameter varying from 0 (complete protection) to 1 (no treatment), depending, e.g. on the amount applied. A second adjustable parameter is the fraction *ρ* of hosts to be treated, giving full (*ρ* = 1) or partial (*ρ* < 1) coverage. We define a particular choice of *ρ* and *ε* as a *control strategy*. The minimal condition for a strategy to be successful is that the point representing the treated system lies on the phase boundary (where *P*_{inv} = 0), because further application of the control agent can only push the point further into the non-invasive regime.

We consider three strategies: (i) a ‘homogeneous’ treatment (only 〈*ψ*〉 is changed), and (ii, iii) two ‘heterogeneous’ treatments (both 〈*ψ*〉 and *σ*_{ψ}^{2} are changed). In the *homogeneous* treatment (figure 5*a*), all the hosts are treated with the same amount of control agent, so that *σ*_{ψ}^{2} = 0 and 〈*ψ*〉 is decreased to the value *ε**ψ*_{0}, reaching the boundary for *ε* = *p*_{c}^{bond}/*ψ*_{0}. In a *heterogeneous* treatment (figure 5*b*,*c*) only a fraction *ρ* < 1 of the hosts is treated, leaving the rest unchanged. The initial homogeneous system becomes a two-host system, with average and variance of *ψ* given by equations (3.1*b*,*c*), where host A is untreated (*ψ*_{A} = *ψ*_{0}) and host B is treated (*ψ*_{B} = *ε**ψ*_{0}). The effect of the control agent on the system is both to decrease the average 〈*ψ*〉 and to increase the variance *σ*_{ψ}^{2}, reaching the phase boundary at a point 〈*ψ*〉 > *p*_{c}^{bond}, *σ*_{ψ}^{2} > 0. The condition that the point lies on the phase boundary given by *σ*_{ψ}^{2} = *p*_{c}^{site} (1 − *p*_{c}^{site})/(*p*_{c}^{site} − *p*_{c}^{bond})〈*ψ*〉 corresponds to a condition for *ρ* and *ε* to lie on a *constraint curve* in the parameter space:
5.1

The constraint curve has endpoints in (1, *ε*_{max}) (corresponding to the homogeneous strategy, figure 5*a*) and (*ρ*_{min}, 0) (maximally heterogeneous strategy, figure 5*c*), where *ρ*_{min} and *ε*_{max} depend on *ψ*_{0}. We remark that in a mean-field approach, where only the average transmissibility would be taken into account, there would be no difference between a heterogeneous strategy that decreases the average to 〈*ψ*〉 = *ρ**ε**ψ*_{0} + (1 − *ρ*)*ψ*_{0} = (1 − *ρ*− *ρ**ε*)*ψ*_{0}, and a homogeneous strategy decreasing the transmissibility of every host by a factor 1 − *ρ*− *ρ**ε*. The difference becomes relevant in spatially explicit systems.

All the strategies are equivalent as far as the cost of control is not considered. If we introduce a cost function *c*(*ρ*, *ε*) (Pinch 1993), the optimal strategy can be defined as the point (*ρ*, *ε*) corresponding to the global minimum of *c*(*ρ*, *ε*) along the curve given by equation (5.1). The optimal strategy depends crucially on the form of the cost function. We analyse here the factorized form:
5.2
where (i) we make the assumption that the total cost is proportional to the coverage *ρ*, while (ii) with the parameter *α* > 0 we allow for a non-linear dependence of the cost on the *treatment efficacy*, defined as the fractional decrease of transmissibility 1 − *ε* for a treated host. The non-linear dependence can be obtained, for example, by assuming that the efficacy of a given amount *Q* of control agent is 1 − *ε* = *Q*^{1/α}, and that the cost of *Q* is proportional to *Q*. A non-linear law for *ρ* can also be considered, but the following results do not change qualitatively. For our choice of the cost function in equation (5.2), the minimum of *c*(*ρ*, *ε*) always coincides with one of the two endpoints of the constraint curve. However, finding which one of the two strategies is optimal is not trivial, and depends on *α* and the initial transmissibility *ψ*_{0}. We can eliminate one of the two parameters through *g*(*ρ*, *ε*) = 0, and consider e.g. *c*(*ρ*, *ε*) = *c*(*ρ*, *ε*(*ρ*)) as a function of *ρ* only. This way, it can be shown (see figure 6*a*) that no local minima for *c*(*ρ*, *ε*) exist as a function of *ρ* (the same occurs for *ε*). We can then compare the values of the cost function at the two endpoints:
and

The homogeneous strategy (1, *ε*_{max}) (figure 5*a*) is optimal when *c*(1, *ε*_{max}) < *c*(*ρ*_{min}, 1), that is, for
5.3

Similarly, the maximally heterogeneous strategy (*ρ*_{min}, 0) (figure 5*c*) is optimal when *α* < *κ*(*ψ*_{0}). For *α* = *κ*(*ψ*_{0}), the two strategies are both optimal with the same cost. In general, *κ*(*ψ*_{0}) is an increasing function of *ψ*_{0} in the interval [*p*_{c}^{bond}, 1] with *κ*(*p*_{c}^{bond}) = 1, and *κ*(1) ≃ 0.77 on the square lattice (see figure 6*b*). Hence, even with a simple choice for the cost function, the optimal control strategy displays a non-trivial dependence on the initial transmissibility of the system. These results are relevant for plant populations, and more generally for any system where the spatial structure of the population is known.

## Acknowledgements

F.M.N., S.N.T., F.J.P.R. and C.A.G. wish to thank BBSRC for funding (Grant no. BB/E017312/1). C.A.G. also acknowledges support of a BBSRC Professorial Fellowship. F.M.N. wishes to thank M. Parry for useful discussions about probability distributions on the unit interval, and T. Mang for a careful reading of the manuscript. The authors thank an anonymous referee for providing analytical insights.

## Footnotes

- Received June 14, 2010.
- Accepted June 24, 2010.

- © 2010 The Royal Society