## Abstract

The endosymbiont *Wolbachia* infects a large number of insect species and is capable of rapid spread when introduced into a novel host population. The bacteria spread by manipulating their hosts' reproduction, and their dynamics are influenced by the demographic structure of the host population and patterns of contact between individuals. Reaction–diffusion models of the spatial spread of *Wolbachia* provide a simple analytical description of their spatial dynamics but do not account for significant details of host population dynamics. We develop a metapopulation model describing the spatial dynamics of *Wolbachia* in an age-structured host insect population regulated by juvenile density-dependent competition. The model produces similar dynamics to the reaction–diffusion model in the limiting case where the host's habitat quality is spatially homogeneous and *Wolbachia* has a small effect on host fitness. When habitat quality varies spatially, *Wolbachia* spread is usually much slower, and the conditions necessary for local invasion are strongly affected by immigration of insects from surrounding regions. Spread is most difficult when variation in habitat quality is spatially correlated. The results show that spatial variation in the density-dependent competition experienced by juvenile host insects can strongly affect the spread of *Wolbachia* infections, which is important to the use of *Wolbachia* to control insect vectors of human disease and other pests.

## 1. Introduction

It has been estimated that up to 60 per cent of all insect species may be infected with the endosymbiotic bacteria *Wolbachia* [1]*. Wolbachia* are vertically transmitted from mothers to offspring, and can spread rapidly when introduced into a novel host population. The bacteria increase in frequency by manipulating their hosts' reproductive system, often using a mechanism known as cytoplasmic incompatibility (CI) [2,3]. CI causes incompatibility between the sperm of *Wolbachia-*infected males and the eggs of uninfected females, leading to failure of embryonic development unless the egg also carries the bacterium [2,4]. The presence of *Wolbachia* carrying males reduces the proportion of viable offspring that are uninfected, allowing the bacterium to spread. Therefore, *Wolbachia* epidemiology is closely linked to the host's demography and patterns of contact between individual hosts.

Currently, *Wolbachia* is being actively investigated as a potential tool to help control certain mosquito-borne diseases. In the mosquito *Aedes aegypti*, which is a vector of important human pathogens, including the dengue, yellow fever and the chikungunya viruses, *Wolbachia* infection can inhibit the development and transmission of human pathogens [5–8]. The bacteria can also reduce mosquito lifespan, which can have a strong impact on the pathogen transmission rate [9]. Proof of principle field releases of *Wolbachia* into wild populations of *A. aegypti* has recently been conducted [10,11]. However, mosquitoes have a complex life cycle in which the adult stage is preceded by juvenile stages that live in completely different microhabitats. Mosquito population dynamics depend on the availability of larval breeding habitat, which can be highly variable in space and time [12–18]. Density-dependent population regulation is also thought most likely to occur in the larval stage [19]. Models that represent these demographic processes are needed to understand better the dynamics of *Wolbachia* infection in these mosquito populations.

Reaction–diffusion models have been used to analyse the spatial dynamics of *Wolbachia* infection [20–22]. This approach assumes that the fitness advantage of *Wolbachia-*carriers in a particular location depends on the local infection frequency, the strength of CI as measured by offspring mortality in incompatible matings, and the fitness costs incurred through *Wolbachia* infection. If the spatial movement of the host population is assumed to be governed by a diffusion process in one dimension (see [21] and the electronic supplementary material) and if the abundance of hosts is assumed constant in space and time, the reaction–diffusion model can be solved to give a simple formula for the equilibrium speed (*v,* measured per host generation) of a travelling wave of *Wolbachia* infection,
1.1where *σ* is the standard deviation of the distance between the birthplace of parents and their offspring and *s*_{h} is the fraction of offspring resulting from an incompatible mating that fails to develop [21,23]. It is assumed that *Wolbachia* infection reduces host fecundity, by a proportion *s*_{f}, and . The quantity represents the trade-off between the fitness costs of *Wolbachia* infection and the fitness advantage afforded by CI. Equation (1.1) shows that if a *Wolbachia* infection that has become established in one location will spread spatially.

While reaction–diffusion models allow the derivation of the equilibrium travelling wave solution, they do not take into account important aspects of the demography of the host population that may affect spread [21,23–25]. In particular, the assumption of constant host population abundance restricts their ability to explore how host population dynamic processes interact with the dynamics of *Wolbachia* infection. For example, the introduction of *Wolbachia* can strongly affect the size of the local host population due to the action of CI, which can have implications for the spread of the infection. In many insect host populations, including mosquito species that are a target for *Wolbachia*-based vector control interventions, host age structure and density-dependent population regulation are important determinants of their dynamics. Non-spatial models of *Wolbachia* dynamics have demonstrated that these processes can strongly influence *Wolbachia* spread. For example, the introduction of insects infected with *Wolbachia* into a novel host population perturbs juvenile host density and the level of juvenile density-dependent competition, which affects the rate of insect introduction that is necessary to trigger spread [26,27].

In this study, a metapopulation model is developed to explore how the demography and dynamics of the host population affects *Wolbachia* spread. A key difference between this modelling approach and the reaction–diffusion model is that the host population size is dynamic and regulated in each subpopulation of the metapopulation (which we shall refer to as a patch) by juvenile density-dependent competition*.* The model also accounts for host age-structure, age-specific effects of *Wolbachia* infection on the host's demography and imperfect maternal transmission of *Wolbachia*. We first assess how well the analytic solution of the reaction–diffusion model equilibrium predicts properties of the spatial dynamics of the metapopulation model, including the equilibrium speed of the travelling wave of *Wolbachia* infection and the conditions for spatial spread to occur. We then explore the influence of spatial heterogeneity in larval carrying capacity on the spread of *Wolbachia*. To produce results that are most relevant to the application of *Wolbachia* to controlling mosquito-borne diseases, the model is parametrized to represent the demography of a mosquito population. Realistic features of the mosquito habitat are incorporated, including spatially correlated stochastic variation in larval habitat quality.

## 2. Model development

An age-structured metapopulation model was developed to explore the spatial dynamics of *Wolbachia* infection in an insect host population. The host population is subdivided into an array of patches connected by migration. A demographic model describes the host population dynamics in each patch. Values of the model parameters are as given in table 1, unless otherwise specified.

### 2.1. Host demography and *Wolbachia* infection

The representation of host age-structure and the effects of *Wolbachia* infection on host demography is similar to that used by Hancock *et al.* [26]. A detailed description of this part of the model is provided in the electronic supplementary material and only the main elements are summarized here.

The insect life cycle is divided into larval and adult stages. Let *L*(*t,l*) be the numbers of larvae at time *t* which have been in the larval stage for time *l*. Let the probability of a larva survives until time *l* be *θ*_{L}(*t,l*) which because larval mortality may be density-dependent is a function of time. Larval mortality is described by a power function, , where is the total larval density at time *t* and *α* and *β* are constants. The parameter α scales the quality of the habitat while higher values of the parameter *β* produce a steeper response to changing density (which we shall refer to as strong density-dependence) [19]. The duration of the larval stage is *T*_{L}.

Let *A*(*t,a*) be the total number of adults of age *a* at time *t.* The probability that an adult survives until age *a* is defined as *θ*_{A}(*a*) which depends on age alone. Let be the total number of adults at time *t* (obtained by integrating over all age classes) and *λ* the female fecundity per unit of time.

The host population is divided into classes that are either infected or uninfected by *Wolbachia,* denoted by subscripts W and U, respectively. Mating is assumed to occur at random, and the frequency of incompatible matings at time *t* is thus As assumed in the reaction–diffusion model, the fecundity of infected mothers is reduced by a proportion *s*_{f} while incompatible matings cause a fraction *s*_{h} of the resulting (uninfected) offspring to die. We also consider the possibility that *Wolbachia* infection may increase adult mortality, particularly in older age-classes, and the proportional reduction in average adult longevity caused by *Wolbachia* is denoted *s*_{g} [26]. Finally, we assume that *Wolbachia* may not be perfectly maternally transmitted, so that a proportion of *ω* of the offspring of infected mothers are uninfected [28].

### 2.2. Spatial population dynamics in one dimension

Host insects are distributed among *N* patches evenly spaced on a straight line. Adult insects move between the patches, and the probability of insects moving *n* patches from where they were born by the time they are of age *a* is given by the dispersal kernal *θ*_{D}(*n*,*a*)*.* Using the superscript *i* to represent spatial location, the model is described by the following system of equations
2.1a
2.1b
2.1cand
2.1d

In any patch of the metapopulation, a threshold infection frequency, , must be exceeded in order for *Wolbachia* to spread. In the special case in which the effects of host movement can be ignored (the patch experiences no immigration or emigration), Hancock *et al.* [26] showed that the unstable equilibrium infection frequency above which *Wolbachia* spreads is
2.2

where *M* = *ω*(1 − *s*_{f}) and *J* = 1 − (1 − *ω*)(1 − *s*_{f})(1 − *s*_{g}). This is identical to the expression for the unstable equilibrium derived by Turelli [28] except that equation (2.2) includes a term for the effect of *Wolbachia* on adult mortality (*s*_{g}). Turelli & Hoffmann's [21] reaction–diffusion model assumes that *Wolbachia* is perfectly maternally transmitted (*ω* = 0) and does not affect host lifespan (*s*_{g} = 0), which gives . In this case, has the same interpretation as in the reaction–diffusion model, and represents the trade-off between the fitness costs of *Wolbachia* infection and the fitness advantage afforded by CI. In both models, the equilibrium condition for local spread of *Wolbachia* is perturbed by the immigration and emigration of hosts, and requires a more complex analysis [22,26,29].

### 2.3. Comparing the metapopulation model with the reaction–diffusion model

We compared the equilibrium speed of the travelling wave of *Wolbachia* infection predicted by the metapopulation equations (2.1*a*–*d*) and reaction–diffusion equation (1.1) models in the case where all patches have identical dynamics and hence equilibrium host abundance is spatially homogeneous. Host dispersal in the metapopulation is represented as a classical random walk which approaches a diffusion process as the number of patches becomes large [30].

We estimated the standard deviation of the distance between the birth place of parents and offspring (*σ* in equation (1.1)) by the average of the standard deviation of the net number of patches moved by individuals dying at age *a* weighted by the probability of death at this time (see the electronic supplementary material). Equation (1.1) assumes that *Wolbachia* does not affect host lifespan, so for the purpose of this comparison, we set *s*_{g} = 0. The average generation time of hosts in the metapopulation is *T*_{L} + *Θ*, where *Θ* is the average adult lifespan. For the parameters in table 1, *Θ* = 20 days and the average standard deviation of dispersal is *σ* = 2.7 patches per generation.

The reaction–diffusion model does not include density-dependent mortality in the host population. The effects of different forms of juvenile density-dependence on the wave speed predicted by the metapopulation model were explored by varying *β*^{i} across metapopulations. The second density-dependent mortality parameter, *α*^{i}, which scales population size, was adjusted so that the equilibrium abundance of adults in the absence of *Wolbachia* was the same for each metapopulation.

### 2.4. Periodic spatial heterogeneity

To explore the effect of spatial (but not temporal) heterogeneity in larval habitat quality on the dynamics of *Wolbachia* infection, we begin by allowing the environment to vary periodically in space [31]. Larval patch quality is determined by the parameter *α*^{i} in the function describing density-dependent mortality with lower values of this parameter corresponding to better-quality patches that allow higher larval survival for a given density. Patch quality is assumed to switch between poor (*α*_{p}) and good (*α*_{g}) values with runs of good patches of length *L*_{g} being followed by runs of poor patches of length *L*_{p}. This produces smooth periodic oscillations in the equilibrium adult abundance across space (see the electronic supplementary material).

We compare populations with different patterns of spatial heterogeneity in two ways. First, we keep the values of *α*_{p} and *α*_{g} the same across all metapopulations. Second, we adjust the values of these parameters so that the amplitude of the periodic spatial variation in the equilibrium abundance of adults is constant across metapopulations (details in the electronic supplementary material). In all simulations, the model is initialized by introducing *Wolbachia-*infected adults at a constant rate into the three adjacent patches at the far left of the linear array until the *Wolbachia* infection frequency in these patches reaches the upper stable equilibrium, after which the introduction ceases.

### 2.5. Spatial variability and correlation

Realistic spatial variation in the quality of larval breeding habitats will be less regular than the periodic variation described earlier. To explore this, we assume that the parameter that determines patch quality, *α*^{i}, in the function describing density-dependent mortality is a realization of a normal random variable with mean *α*_{m} and standard deviation *σ*_{m} (table 1). We allow for the quality of nearby patches to be similar by assuming the correlation coefficient of larval habitat quality in patches separated by *n* patches is e^{−rn}, where *r* is a measure of the spatial extent of the correlation. The spatial correlation range is defined as the maximum separation distance in number of patches for which the correlation coefficient is greater than 0.05. Spatially correlated values of *α*^{i} for all patches were generated by taking a Cholesky decomposition of the correlation matrix *C _{ij}* = e

^{−r|i−j|}[32]. We compare metapopulations in which the mean and standard deviation of

*α*

^{i}is the same but the range of spatial correlation varies.

### 2.6. Model parametrization

The default set of parameters used in the model runs are shown in table 1. They were chosen to represent as closely as possible a typical *Wolbachia* infection in an *A. aegypti* mosquito population (for further details, see the electronic supplementary material and Hancock *et al.* [33]). *Wolbachia* infection was assumed to cause a small reduction in average adult lifespan, with older insects experiencing enhanced mortality owing to endosymbiont carriage.

## 3. Results

### 3.1. Wave speed in a homogeneous environment

The equilibrium speed of the travelling wave of *Wolbachia* infection derived from the reaction–diffusion model equation (1.1) provides a good approximation to that of the more complex metapopulation model (figure 1). The agreement is best for low values of the unstable equilibrium infection frequency, and is also affected by the form of juvenile density-dependent mortality, especially for higher . The reaction–diffusion model predicts that the speed of spread declines to zero as approaches 0.5; the metapopulation model indicates that spread stops at lower values of for weak density-dependence (low *β*) but the limit is above 0.5 for strong density-dependence (high *β*). Other parameters being equal; it is thus more difficult for *Wolbachia* to spread when density-dependence is weak.

The reasons why the form of juvenile density-dependence affects the spatial dynamics of *Wolbachia* are complex. The immigration of infected individuals into a patch has opposing effects on juvenile density; the introduced infected females reproduce and increase larval density, but the resident uninfected females produce fewer offspring on average due to incompatible matings with the introduced infected males. The immigration of infected individuals may thus either increase or decrease larval density depending on the rate at which they arrive and the local abundance of infected and uninfected insects. The influence of the form of density-dependence on the speed of spatial spread depends on whether immigration causes larval density to increase or decrease. If the density of larvae decreases, their mortality will be comparatively lower if density-dependence is strong. More juveniles will survive to become adults and there will be a greater number of infected adults able to disperse to other patches and initiate *Wolbachia* spread. If density-dependence is weak then the larval mortality remains relatively high even though their density is reduced. This means that fewer infected adults disperse to other patches. The effect of the form of density-dependence on spatial spread is more marked when *Wolbachia* infection incurs a higher fitness cost because a greater rate of immigration of infected insects is then required to allow spread to occur.

### 3.2. Periodic variation in larval breeding habitat

Periodic spatial heterogeneity reduces the average wave speed at which the *Wolbachia* infection spreads. To explore this consider first compare metapopulations where different length runs of good- and poor-quality patches occur periodically (with the values of the patch quality *α*_{g} and *α*_{p} being the same across metapopulations). The abundance of adults approaches spatial homogeneity when the scale of environmental heterogeneity in larval habitat is finer than the standard deviation of adult dispersal per generation, *σ*. For the parameters in table 1, the equilibrium wave speed in a homogeneous environment is approximately one patch per generation.

In figure 2, wave speed is plotted as a function of the length of the runs of good- or poor-quality patches. Where adult densities do not vary spatially the wave speed remains approximately one patch per generation irrespective of the actual value of the density (as also predicted by reaction–diffusion models). This occurs when the environment is largely made up of good patches (towards the top left of figure 2) or bad patches (towards the bottom right) or where the periodicity of habitat change is short and adult movement homogenizes density (towards the bottom left).

When good- and poor-quality runs are more widely separated the wave speed markedly declines. This happens when short stretches of good patches are interspersed among longer stretches of poor patches, when the reverse occurs, and also when the length of runs of both types of patch increase together. Habitat heterogeneity thus reduces the speed with which *Wolbachia* infections can spread through a naive host population.

Increasing the spatial scale of variation in larval patch quality increases both the amplitude and period of the variation in adult density. To separate these effects, we again alternated the quality of larval patches periodically but adjusted *α*_{g} and *α*_{p} such that the amplitude of the variation in adult density remained constant across metapopulations (figure 3). This is not possible for fine-scale variation in patch quality hence the different axes ranges in figures 2 and 3. The presence of fluctuations in adult density always substantially reduces the wave speed below that observed in homogeneous environments. Spatial spread of *Wolbachia* only occurs when the length of the runs of poor-quality patches exceeds that of the high-quality patches, and longer sequences of poor-quality patches usually produce faster wave speeds. The finest scale of heterogeneity at which the constant amplitude could be obtained involved runs of five to seven high-quality patches and seven to 12 low-quality patches (figure 3). In this region, shorter sequences of poor-quality patches lead to faster wave speeds.

To understand why spatial variation in larval habitat quality slows the spread of infection consider the dynamics of a single patch in a metapopulation. *Wolbachia* will increase in frequency in the patch if the proportion of individuals infected exceeds the unstable equilibrium frequency (see also §4). As the wavefront nears the patch, the rate of immigration of infected insects (*I*_{W}) increases leading to this threshold being breached; however, the patch dynamics are also affected by immigration of uninfected individuals (*I*_{U}). Immigration of uninfected adults (that then reproduce) increases the proportion of larvae that are uninfected and reinforces the relative recruitment of uninfected adults. The threshold rate of immigration of infected adults required for spread (*I*_{WT}) and the unstable equilibrium infection frequency that is breached at this immigration rate can be calculated using a similar method to that developed by Hancock *et al.* [26] (see figure 4 and the electronic supplementary material). Both these quantities increase as the rate of immigration of uninfected adults increases (figure 4), therefore higher rates of uninfected insect immigration require a greater influx of infected individuals before *Wolbachia* can spread.

The wavefront can be slowed or stopped if it approaches a region of good-quality larval patches that are a source of high rates of immigration by uninfected individuals. This is illustrated in figure 5, which shows a wave of infection moving from left to right through an environment that consists of low- and high-quality patches in the ratio 4 : 1. The local speed of spread is slowest when the wavefront spans a region of poor-quality patches but is approaching a sequence of high-quality patches. At this point, the wave produces relatively low rates of immigration of infected individuals because they are recruited in low-quality patches. The immigration of uninfected insects from patches ahead of the wave is high as they are coming from good patches. However, once the infection begins to spread in the high-quality region, the reverse effect occurs. Now the high-quality patches are producing infected immigrants and this increases the speed of the wave, especially when there are low-quality patches ahead of the wave and fewer uninfected immigrants are being produced.

Longer regions of good-quality patches provide a greater barrier for spread as they increase the rate of immigration of uninfected individuals into low-quality patches. Greater differences in quality between the two patch types also impedes spread by increasing the relative magnitude of the rate of immigration of uninfecteds into low-quality patches.

An increase in the length of the sequence of low-quality patches has two effects on the speed of spread. First, and most obviously, the average speed increases because fewer barriers of regions of high immigration of uninfected insects are encountered. But there is a second more subtle effect. In a longer region of poor patches, the wave can accelerate to its maximum speed and assume its equilibrium shape which has a steep spatial gradient in infection frequency (a steep wavefront; figure 6). This means that at the wavefront a greater fraction of the immigrants moving ahead consist of infected individuals. This assists the infection to spread into regions of high-quality patches. Therefore, the wave is less likely to be halted by barriers of high-quality patches if it encounters them after travelling through a long region of relatively poor-quality patches (figure 3). Of course, when the length of runs of low-quality patches is within the insect's dispersal range, the travelling wave can ‘jump the gap’ and hence the wave speed can increase (figure 3).

### 3.3. Spatial variability and correlation

While the study of *Wolbachia* spread through a periodic environment provides clear insight into the processes involved, real environments are more irregularly structured. To explore this, we studied movement through a metapopulation where patch qualities were chosen from a random distribution with given mean, variance and spatial covariance.

An example of a set of simulations is shown in figure 7. Keeping the mean and variance constant, increasing spatial covariance in patch quality reduced the average speed of the wave of *Wolbachia* infection and also increased the variance in speed across replicate simulations. The probability that the wave comes to a stop is also higher when there is longer-range spatial covariance in habitat quality.

Longer-range spatial covariance increases the probability that runs of low-quality patches are followed by sequences of high-quality patches which reduces the speed of the travelling wave for the reasons given in the discussion of periodic environments. In this stochastic setting, environments will occur where the travelling wave is brought to a halt by an insuperable barrier of high-quality patches, and these occur more often where habitat quality is correlated over longer distance. The variance in the speed of the wave is higher in spatially correlated environments because the wave speed increases across long regions of similar patch quality but decreases when long sequences of higher patch quality are encountered. An interesting observation is that realizations of the metapopulation occur in which the speed of spread exceeds that in a homogeneous environment when patch quality is spatially correlated (figure 7). The reason for this is that by chance these environments contain regions across which habitat quality declines, leading to long regions where adult host abundance consistently decreases in the direction the wave is travelling. Because it is easier for *Wolbachia* carried by immigrants to spread from regions of high to low density the speed of spread can be higher in these regions than when adult densities are constant. However, this effect is much weaker than the deceleration in wave speed that occurs in the opposite situation where the immigrants travel from low- to high-density regions.

## 4. Discussion

Endosymbiotic bacteria such as *Wolbachia* that spread by manipulating their hosts' reproduction can have dramatic effects on the biology of their insect hosts which can be harnessed for the management of insect pests and vectors of disease [9,34,35]. A good understanding of *Wolbachia*'s spatial epidemiology is required both to explain their distribution in natural insect communities as well as how they might be manipulated [20,21,34]. The aim of this study was to advance our understanding of these processes by developing a model that included interactions between *Wolbachia* infection frequency, host population demography and spatial variation in host habitat quality.

For appropriate limiting cases of the metapopulation model, we show that the major insights from reaction–diffusion models of *Wolbachia* spread carry over to a more complex model that incorporates realistic aspects of the host population age-structure and density-dependent regulation. In particular, there is good agreement in the conditions for spatial propagation of infection and the equilibrium speed of the travelling wave [21,23]. However, for a wide range of parameter values, details of host population ecology such as the form of density-dependent regulation can affect whether spread occurs and at what speed.

We also show that spatial heterogeneity in larval habitat quality nearly always reduces the speed at which *Wolbachia* spreads through a landscape. The reason for this is that it is difficult for an infection to spread from regions of low- to high-population density. Our analysis has shown that this is due to two main factors. First, a region of poor-quality habitat produces comparatively few migrants which makes it hard for the threshold for *Wolbachia* spread to be exceeded in an adjacent better-quality region. Second, movement of uninfected individuals from high-quality regions into lower-quality areas increases the threshold frequency that must be overcome for local *Wolbachia* establishment in these areas. This means that spatial covariance in habitat quality relative to the typical range of host dispersal is critical in determining the rate of spread of the infection. It also means that the spread of the infection can be halted within a region of low host density. Our results extend those obtained using island–mainland models of *Wolbachia* spread, which show that immigration of hosts from an uninfected mainland increases the threshold frequency required for spread on an island [22,29]. Extension of these models to describe two patches connected by asymmetric two-way migration can mimic the effects of heterogeneous host density by specifying different rates of emigration from each patch [29]. These models assume that the dynamics of *Wolbachia* depend only on their infection frequency in a host population of infinite size with discrete generations.

The same processes underlying the spatial spread of *Wolbachia* have also been associated with the spread of genetic systems involving underdominant chromosomal arrangements [23,28]. At a single location, an advantageous chromosomal form can increase in frequency, but only above a threshold frequency. Several cases where underdominant chromosomal types coexist in stable narrow ‘hybrid zones’ that coincide with regions of low population density have been documented [36–38]. These have been interpreted as regions where the spatial spread of a superior chromosomal type has become arrested because too few migrants are produced to propagate further the wave [37]. In general, dynamics of this type are expected to occur whenever spread is through a bistable wave where at a single location a threshold frequency or density must be breached before spread occurs [22,23,28,39,40].

The potential for *Wolbachia* spread to be halted by movement of uninfected insects from ahead of the wavefront has implications for the design of strategies for deliberate *Wolbachia* introductions. Releasing infected insects may not succeed in triggering *Wolbachia* spread if the liberated insects need to compete with wild-type insects for the available breeding habitat. Spread may be assisted by suppressing wild-type populations in areas that are near the release site, within the insects' dispersal range, so providing breeding habitat for the introduced insects.

In this study, our analysis of the spatial spread of *Wolbachia* was restricted to one-dimensional scenarios in order to compare the results of our metapopulation approach to the analytic results obtained from reaction–diffusion models and develop a clear understanding of the models behaviour. Barton & Turelli [22] analyse *Wolbachia* dynamics using reaction–diffusion models in one and two dimensions and show that the asymptotic speed of the travelling wave in a homogeneous environment is the same in both cases. However, fewer analytical results are available for the two-dimensional case.

Further study of *Wolbachia* spread will require a better understanding of their hosts' biology to produce more realistic models. Where the spread of *Wolbachia* or other endosymbiont infections through insect populations has been observed in nature, discrepancies between observed and predicted movement rates have been attributed, at least in part, to the spatial complexity of field populations that was not captured by the model [21,34]. In *A. aegypti* mosquitoes, the current primary targets for *Wolbachia*-based vector control interventions, larval breeding habitats can be highly heterogeneous, with a small number of sites generating a large fraction of the total recruitment to the adult population [12]. However, during rainy periods, larval breeding habitats may become more widespread [10,41].

The metapopulation model, we have developed here is a step towards incorporating the complex spatio-temporal dynamics displayed by insect host populations into models of *Wolbachia* spread. Models that aim to predict *Wolbachia* dynamics in particular host species at specific geographical locations will need to incorporate further details of the hosts' demography and will be strongly reliant on field data in order to describe the complex ecology of *A. aegypti* and other target species [42]. Our metapopulation model is purely deterministic but stochastic effects can be important in regions of low habitat quality where they are likely to have a significant influence on *Wolbachia* spread due to the small size of the local host insect population [43]. In subdivided host populations where migration is limited, local fluctuations in infection frequency may increase the likelihood of *Wolbachia* invasion [44].

There are very few insect species for which we have a good understanding of movement patterns. In this study, we used a random walk model for adult dispersal to allow comparison between the metapopulation model and the reaction–diffusion approach. For *A. aegypti* mosquitoes, Russell *et al.* [45] estimated the average lifetime dispersal to be *σ* = 78 m, in which case the range of wave speeds reported in figure 7 scales to approximately 0–46 m per generation. However, Schofield [20] showed that the form of the dispersal kernel can have a large effect on the speed of *Wolbachia* spread. Moreover, movement patterns will depend on habitat heterogeneity, and will likely be oriented towards preferred breeding habitats [46], which will influence the likelihood and speed of spread.

A priority for further study of *Wolbachia* dynamics is to combine modelling and experimental approaches to examine the demographic and behavioural processes that drive spatio-temporal fluctuations in the abundance of the insect host population [47,48], and to obtain estimates of the critical biological and environmental parameters that determine the probability and speed of invasion. An understanding of the fitness effects of *Wolbachia* on its host will be essential to this process. Our model assumes that *Wolbachia* infections have fitness costs [9,11] but there is evidence that the bacteria may increase the fitness of its host when they are infected with pathogens [49], which could help to overcome the barriers to the spatial propagation of *Wolbachia* that we have discussed [50].

## Acknowledgments

The research leading to these results has received funding from the European Union's Seventh Framework Programme (FP7/2007–2013) under grant agreement no. 223241, ‘AnoPopAge’. We are also grateful for support from the Wellcome Trust, the Bill and Melinda Gates Foundation and the RAPIDD (Research and Policy in Infectious Disease Dynamics) programme. We thank Steve Sinkins, Michael Turelli, Scott O'Neill and the members of the Bill and Melinda Gates Foundation Grand Challenge Projects ‘Modifying mosquito population age structure to eliminate dengue transmission’ for helpful discussion. Finally, we thank three anonymous reviewers for helpful suggestions leading to improvement of the manuscript.

- Received March 27, 2012.
- Accepted May 15, 2012.

- This journal is © 2012 The Royal Society

This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.