Royal Society Publishing

Heterogeneity in susceptible–infected–removed (SIR) epidemics on lattices

Franco M. Neri, Francisco J. Pérez-Reche, Sergei N. Taraskin, Christopher A. Gilligan

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 [1]. 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 [24]. 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 [5].

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 [6,7], or, at a larger scale, pathogens spreading within a mosaic of neighbouring susceptible fields or even farms [8,9]. 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 [10]. Finally, at the microscopic scale, some diseases are known to spread among contiguous cells (cf. neural ganglion cells in [11]). Nearest-neighbour interactions lead naturally to consideration of percolation [12]. 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 [13,14], and the connectivity of habitable sites in the landscape [7]. 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 [1517] and recent interest in the field of complex networks [1821], 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 [22]. The model is chosen because an exact mapping onto bond percolation holds in the homogeneous case [14], 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 pcbond. This also implies that in a few relevant cases the value of the invasion threshold is known analytically [12]. 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 [16] 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 〈ψ〉 = pcbond [17]. 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 〈ψ〉 = pcbond (in homogeneous and heterogeneous uncorrelated systems) to 〈ψ〉 = pcsite (site percolation threshold), when the variance of the transmissibility σψ2 in correlated systems is maximal [15,16,23]. 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 [pcbond, pcsite] 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:Embedded Image 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 [24], 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 Pinv (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 [12,14], 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,Embedded Image 2.2where pcbond is the bond-percolation threshold (pcbond = ½ for the square lattice; [12]). In an infinite system, the critical value separates an invasive region in the parameter space (Pinv(β, τ ) > 0 for ψc> pcbond) from a non-invasive, safe region (Pinv(β, τ) = 0 for ψc< pcbond). 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(β), [17]). 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 rangeEmbedded Image 2.3with bounds corresponding to the bond- and site-percolation thresholds pcbond and pcsite. The system is always in the non-invasive regime for 〈ψτ < pcbond, and always in the invasive regime for 〈ψτ > pcsite (whether the system is in the invasive regime or not for 〈ψτ = pcbond and 〈ψτ = pcsite depends on the topology of the lattice and is still an open problem: e.g. [25]); we do not discuss these cases here). However, in the middle interval pcbond < 〈ψτ < pcsite 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, pcbond = pcsite on a tree graph, while pcbond< pcsite for lattices with dimension d ≥ 2 [25]. For the square and triangular lattices, the interval (pcbond, pcsite) 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 pcbond is obtained when the ψij are uncorrelated, while the value pcsite 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 = σψ,c2. 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 [26,27], 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 ψ areEmbedded Image 3.1aEmbedded Image 3.1bEmbedded Image 3.1c

Here and in the rest of the paper, δ(x) is the Dirac delta-function with the property ∫−∞+∞ f(x)δ(xx0)dx = f(x0). 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 〈ψ〉 = pcbond and 〈ψ〉 = pcsite, with 〈ψ〉 given by equation (3.1b). 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.1c). 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 [1921] 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).

Figure 1.

Effect of heterogeneity on the phase diagram of a two-host system on a square lattice. Each frame corresponds to a different value of ρA (concentration of A hosts; (a) ρA = 0.05; (b) ρA = 0.3; (c) ρA = 0.5). The non-invasive and invasive regimes are in green and orange, respectively. The red line marks the phase boundary, and the black solid lines are the upper and lower bounds for the phase boundary (〈ψ〉 = pcbond, pcsite, respectively, from equation (2.3)). Along the phase boundary, the variance of h(ψ) increases from 0 (at the homogeneous point ψA = ψB = pcbond) in the direction indicated by the arrows. As the level of heterogeneity ρA is increased (from left to right), the phase boundary moves away from the lower bound and approaches the upper bound, the minimum distance being reached at the edges (ψA = 1,ψB = 0) and (ψA = 0,ψB = 1), where the variance of h(ψ) is maximal.

The results discussed so far, and the theorems giving bounds for heterogeneous SIR processes [15,16], 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.

Figure 2.

Comparison of the invasive regions of random and periodic two-host systems on a square lattice with given values of ψA, ψB, and ρA: (a) ABAB systems with Embedded Image; (b) AABAAB systems with Embedded Image. The black solid lines are the upper and lower bounds for the phase boundary of the random system (compare figure 1). The phase boundary for the random system (red line) always lies below that for the periodic system (blue line), so that the region between the two curves (shaded) is non-invasive for the periodic system and invasive for the random one, and the periodic system is in general more resilient to invasion. The orange and green regions are invasive and non-invasive, respectively, for both systems.

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 pcbond, corresponding to zero variance in ψ, and the upper bound, pcsite, corresponding to the maximal variance for a given value of 〈ψ〉,Embedded Image which is achieved for a particular two-host system with h(ψ) = 〈ψδ(1 − ψ) + (1 − 〈ψ〉)δ(ψ) (cf. equation (3.1a) and the electronic supplementary material). Therefore, when the average transmissibility falls within the window pcbond ≤ 〈ψ〉 ≤ pcsite, the system can be in different regimes (invasive or non-invasive) depending on the variance of the distribution, 0 ≤ σψ2σψ,max2. 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 Pinv on 〈ψ〉 and σψ2 for the family of bimodal distributions (equation (3.1a)) with ψB = 0 is shown in figure 3. In general, the probability of invasion Pinv is an increasing function of 〈ψ〉 and a decreasing function of σψ2 (see figure 3a). 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 3a,b), the transition occurs at a value 〈ψc such that pcbond ≤ 〈ψcpcsite. In the second case (green line in figure 3a,c), it occurs for a σψ,c2 such that 0 ≤ σψ,c2σψ,max2. The bound σψ2 = 0 (where Pinv > 0) corresponds to a homogeneous process equivalent to bond percolation, while σψ2 = σψ,max2 (where Pinv = 0) corresponds to a maximally correlated process equivalent to site percolation.

Figure 3.

Probability of invasion Pinv for a process with site heterogeneity on a square 1000 × 1000 lattice. (a) Pinv as a function of 〈ψ〉 and σψ2 . Each point in the plane (〈ψ〉, σψ2) corresponds to a choice of the parameters ψA and ρA in a family of bimodal distributions for ψ (equations (3.1)) with ψB = 0. The value 〈ψ〉 is chosen to lie in the interval (pcbond, pcsite) and the variance σψ2 varies from 0 to the maximal value 〈ψ〉 (1 − 〈ψ〉). The phase boundary between invasive (orange) and non-invasive (green) regions is shown below the surface. (b) Pinv as a function of 〈ψ〉, with constant variance σψ2 = 0.1 (same as red curve in (a)). (c) Pinv as a function of σψ2, with constant average 〈ψ〉 = 0.55 (same as green curve in (a)). Note that the horizontal axis in (c) is inverted.

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 4a). 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 (pcbond, 0) (homogeneous system, bond percolation case) to (pcsite, pcsite (1 − pcsite); maximally heterogeneous system, site percolation case). A system with pcbond ≤ 〈ψ〉 ≤ pcsite can be in the invasive or in the non-invasive regime depending on the value of σψ2.

Figure 4.

(a) Phase diagram in the (〈ψ〉, σψ2 ) plane for a system with site heterogeneity. The points of the phase boundary separating the invasive and non-invasive regions (red points) were calculated numerically for the family of bimodal probability distributions for ψ given by equation (3.1a) with ψB = 0. The accessible points lie below the curve σψ2 = 〈ψ〉 (1 − 〈ψ〉). The values of σψ2 are exact; the error on the values of 〈ψ〉 is about 0.01Δp, Δp = pcsitepcbond. (b) Points of the phase boundary on the same lattice for different families of distributions: bimodal with ψB = 0 (red circles), β distribution (green triangles), truncated Gaussian (blue crosses). The difference along the 〈ψ〉 axis between points for different distributions is less than 0.03Δp. The points lie approximately on a straight line (dashed black line).

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 4a. The position of the phase boundary in figure 4a 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 4b 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 = pcsitepcbond (cf. figure 4b). 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 (pcbond, 0) and (pcsite, pcsite (1 − pcsite)).

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 [15].

Our discussion is focused on the position of the invasion threshold. We do not address here the question of the dependence of Pinv on heterogeneity when the system is in the invasive regime (see [15,18,19]). 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 = 〈ψ〉 > pcbond, 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 Pinv = 0), because further application of the control agent can only push the point further into the non-invasive regime.

Figure 5.

Effect of different control strategies on a square-lattice homogeneous system with transmissibility ψ0 = 0.7, represented by a black circle on the (〈ψ〉, σψ2) plane. The red circle represents the system after treatment. The phase boundary between the invasive (orange) and non-invasive (green) regions is approximated by a straight line. The effect of treatment on a host is to decrease its transmissibility to a value εψ0, 0 ≤ ε< 1. (a) Homogeneous strategy (ρ = 1, ε = εmax = pcbond /ψ0). All the hosts are treated. The final system is homogeneous with average transmissibility pcbond . (b) Heterogeneous strategy (ρ = 0.4, ε = 0.35). A fraction ρ of the hosts is treated. The system reaches the phase boundary at a point with 〈ψ〉 > pcbond, σψ2 > 0. (c) Maximally heterogeneous strategy (ρ = 0.24, ε = 0). A smaller fraction of the hosts is treated, their transmissibility being decreased to 0. The variance of the treated system is maximal for the initial condition ψ0.

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 5a), 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 ε = pcbond/ψ0. In a heterogeneous treatment (figure 5b,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.1b,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 〈ψ〉 > pcbond, σψ2 > 0. The condition that the point lies on the phase boundary given by σψ2 = pcsite (1 − pcsite)/(pcsitepcbond)〈ψ〉 corresponds to a condition for ρ and ε to lie on a constraint curve in the parameter space:Embedded Image 5.1

The constraint curve has endpoints in (1, εmax) (corresponding to the homogeneous strategy, figure 5a) and (ρmin, 0) (maximally heterogeneous strategy, figure 5c), 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(ρ, ε) [28], 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:Embedded Image 5.2where (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 − ε = Q1/α, 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 6a) 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:Embedded Image andEmbedded Image

Figure 6.

Optimal control strategy for a square-lattice system with initial constant transmissibility ψ0, with treatment parametrized by ρ and ε, and cost function given by equation (5.2). (a) The cost c(ρ, ε) as a function of ρ for ψ0 = 0.7 and α = 0.5 (blue line), α = κ (ψ0)≃1.16 (magenta), α = 2 (yellow). The α = κ (ψ0) line has a maximum around 0.4. (b) The values of ψ0 and α determine whether the optimal strategy is the homogeneous one (1,εmax) (grey region) or the maximally heterogeneous one (ρmin, 0) (yellow region). Along the curve α = κ (ψ0) (red line), both strategies are optimal.

The homogeneous strategy (1, εmax) (figure 5a) is optimal when c(1, εmax) < c(ρmin, 1), that is, forEmbedded Image 5.3

Similarly, the maximally heterogeneous strategy (ρmin, 0) (figure 5c) 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 [pcbond, 1] with κ(pcbond) = 1, and κ(1) ≃ 0.77 on the square lattice (see figure 6b). 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.

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

References

View Abstract