- Email: [email protected]
Spatially Explicit Models of Turelli-Hoffmann Wolbachia Invasive Wave Fronts
J. theor. Biol. (2002) 215, 121–131 doi:10.1006/jtbi.2001.2493, available online at http://www.idealibrary.com on
Spatially Explicit Models of Turelli–Hoffmann Wolbachia Invasive Wave Fronts Peter Schofield Department of Mathematics, University of Dundee, DD1 4HN, U.K. (Received on 5 January 2001, Accepted in revised form on 2 December 2001)
This paper examines different mathematical models of insect dispersal and infection spread and compares these with field data. Reaction–diffusion and integro-difference equation models are used to model the spatio-temporal spread of Wolbachia in Drosophila simulans populations. The models include cytoplasmic incompatibility between infected females and uninfected males that creates a threshold density, similar to an Allee effect, preventing increase from low incidence of infection in the host population. The model builds on an earlier model (Turelli & Hoffmann, 1991) by incorporating imperfect maternal transmission. The results of simulations of the models using the same parameter values produce different dynamics for each model. These differences become very marked in the integro-difference equation models when insect dispersal patterns are assumed to be non-Gaussian. The success or failure of invasion by Wolbachia in the simulations may be attributed to the insect dispersal mechanism used in the model rather than the parameter values. As the models predict very different outcomes for the integro-difference models depending on the underlying assumptions of insect dispersal patterns, this emphasizes that good field data on real (rather than idealized) dispersal patterns need to be collected before models such as these can be used for predictive purposes. r 2002 Elsevier Science Ltd.
Introduction The intra-cellular organism Wolbachia is a much studied reproductive parasite. Primarily vertically transmitted in egg cytoplasm, it manipulates its host’s reproduction in ways that increase the proportion (or fitness) of infected females in the population and hence the proportion of infected offspring in the next generation. The details vary depending on host species and Wolbachia strain (Werren, 1997; O’Neill et al., 1997). Turelli and Hoffmann studied the sweep of Wolbachia through the Drosophila simulans population in California in a series of papers (Hoffmann et al., 1986; Turelli & Hoffmann, 1991, 1995; Turelli et al., 1992). The Wolbachia 0022-5193/02/050121 + 11 $35.00/0
strain studied caused a reproduction distorting effect referred to as cytoplasmic incompatibility (CI). First observed in mosquitoes (Hertig & Wolbach, 1924) and attributed to infection (Ghelelovitch, 1952), CI has the effect that crosses between infected males and uninfected females suffer an increased hatch rate failure. It is thought that Wolbachia imprints the sperm of infected males in a manner that requires infection in the egg to recover the sperm function (Breeuwer & Werren, 1990). Turelli showed that there is a threshold level in the dynamics of invasion by the infected phenotype into an uninfected population (Turelli & Hoffmann, 1995) which produces local dynamics similar to an Allee effect, where the unstable equilibrium r 2002 Elsevier Science Ltd.
122
P. SCHOFIELD
threshold must be exceeded for invasion to be successful. There are also two stable equilibria with either no infection or infection fixation. Imperfect maternal transmission has two effects: (1) it produces a stable equilibrium where uninfected individuals are maintained in the population; (2) it increases the level of the unstable equilibrium or threshold for invasion (Turelli & Hoffmann, 1995). In these models, we assume that the level of horizontal transmission of infection is negligible (Turelli et al., 1992). One of the main motivations for the current investigation is the observation of a wave-like advance of Wolbachia infection through the California D. simulans populations during the late 1980s (Turelli & Hoffmann, 1991). In their paper, Turelli and Hoffmann suggest that in order for the reaction–diffusion equation model’s results to fit the data, the D. simulans dispersal rate s must be of the order of 50 km gen1/2. They comment that this is almost an order of magnitude greater than the observed dispersal rates in mark-recapture experiments (Jones et al., 1981) and this might indicate that some dispersal process other than diffusion is acting. Simple diffusion is therefore unlikely to be the correct dispersal mechanism to use. However, it must be noted that that mark-recapture data in Jones et al. (1981) is for D. pseudoobscura rather than D. simulans, and therefore the lack of fit may be due to species differences rather than the diffusion asumption. It has also been shown that a convection process, such as prevailing wind, can increase the wave speed but will not otherwise alter the shape of the wave front (Murray, 1989; Lewis & Kareiva, 1983). Another approach that has been employed recently in modelling insect dispersal (in particular Drosophila sp. movements) is that of integro-difference equations (Kot et al., 1996). The work of Kot et al. (1996) reanalyses Drosophila pseudoobscura dispersal data from Dobzhansky & Wright (1947), showing that the data were not sufficient to distinguish between exponentially bounded and non-exponentially bounded dispersal functions. In an exponentially bounded distribution, the probability of movement from the initial location to a new location (within a fixed time) approached zero exponentially or faster with the increase in distance
between the two locations. With the nonexponentially bounded distribution, this is not the case. Long-range redistribution events have a finite probability, see Fig. 1. The non-exponentially bounded case could represent a situation where most movement is extremely local, for example, insects with weak flight, but there is a significant probability of much longer range movement, perhaps by human activity such as fruit transportation. Kot et al. (1996) also investigated the Allee effect showing that this slows down the effective wave speed. This has motivated the present examination of the California Wolbachia sweep using integro-difference models with contrasting dispersal kernels. Mathematical Models Models of invasion using Fisher-type reaction diffusion equations (Fisher, 1937; Murray, 1989), integro-difference equations (Mollison, 1977) and stepping stone models (Kimura, 1953; Renshaw, 1991), that produce a travelling wave of advance have been well studied. Analytical results have been obtained for diffusion equation models of the advance of advantageous genes (Fisher, 1937; Murray, 1989), movement in hybrid zones (Nagylaki, 1975; Barton, 1979) and invasion in predator–prey systems (Dunbar, 1983). The dynamics of travelling waves from a Fisher-type equation that included an Allee effect have also been investigated (Lewis & Kareiva, 1993). The reproduction distortion mechanism of CI, studied in Hoffmann and Turelli’s work on D. simulans, also produces an effect similar to an Allee effect. This means that if infected males are at a high frequency in the population, then uninfected females are at a disadvantage. If the infection is above a critical frequency then infection will increase to fixation, that is, all members of the population are infected. If there is imperfect maternal transmission fixation is not reached, but rather an equilibrium balance between fitness increase due to infection and infection loss is achieved. Turelli & Hoffmann (1991) assumed perfect maternal transmission and a fecundity penalty for infected females. This produced the following equation for the change Dp in the proportion of infected
SPATIALLY EXPLICIT MODELS
123
Fig. 1. Subplots (a) and (b) show the graph of the Gaussian distribution function with dispersal rate s=50 km gen1/2. Subplots (c) and (d) show the graph of the leptokurtic distribution function with dispersal rate s=5 km gen1/2. Comparison of the graphs (b) and (d) shows that the Gaussian function tends to zero faster than the leptokurtic function.
individuals p in the population between each generation: Dp ¼
# sh pð1 pÞðp pÞ ; 1 sf p sh pq
ð1Þ
where p# ¼ sf =sh ; sf is the reduction in fecundity caused by Wolbachia infection and sh is the reduction in hatch success in incompatible crosses (infected male cross uninfected female). Turelli and Hoffmann make the approximation that the terms s f p and sh pq are negligible and hence mean fitness can be assumed to be 1. The addition of spatial movement by a diffusion mechanism gives the following equation for the spatial spread of infected individuals: @p @2 p # ¼ D 2 þ sh pð1 pÞðp pÞ; @t @x
ð2Þ
where D is the diffusion coefficient. It is assumed that D=s2/2 where s is the variance of the individual movement probability distribution (Okubo, 1980; Turchin, 1998).
Turelli and Hoffmann use the analytical result that in the limit t-N (i.e. time tends to infinity) eqn (2) produces travelling wave fronts that can be approximated by a hyperbolic tangent function (Barton, 1979). This models the asymptotic constant velocity wave of advance of infection through the population. Here I examine a modified model that extends eqn (2) to include an imperfect maternal transmission rate m producing the equation @p @2 p # mð1 sf Þp: ð3Þ ¼ D 2 þ sh pð1 pÞðp pÞ @t @x If m=0 then eqn (3) reduces to eqn (2). In the next section, I compare numerical solutions of the reaction–diffusion eqn (3) with the solution using an asymptotic approximation of the equation, similar to Barton (1979), when m=0, given by " pffiffiffiffi !# # sh t=2 x sð1 2pÞ 1 pffiffiffiffi : pðx; tÞ ¼ 1 tanh 2 2s= sh ð4Þ
124
P. SCHOFIELD
These (continuous) reaction–diffusion models can be compared to other models that use a different explicit spatial modelling approach, that of discrete-time, integro-difference equation models (Kot & Schaffer, 1986; Kot, 1992), where movement is modelled by dispersal at the end of each discrete generation using a probability density function. Integro-difference models are based on equations of the generic form ptþ1 ðxÞ ¼
Z
b
kðjx yjÞ f ð pt ð yÞÞ dy;
ð5Þ
a
where k(|xy|) is the dispersal kernel and f( pt(y)) a function representing the local kinetics at spatial location y. The particular form I derived using the kinetics of eqn (3) is as follows: ptþ1 ðxÞ ¼
Z
b
kðjx yjÞ½sh pt ð yÞð1 pt ð yÞÞ a
# mð1 sf Þpt ð yÞ dy: ð6Þ ð pt ð yÞ pÞ It has been shown that the shape of the dispersal kernel can have significant effects on the resulting wave of invasion (Kot et al., 1996). Using eqn (6), I compare results for two probability density functions: (I) The Gaussian (normal ) distribution function: 2 1 2 kðjx yjÞ ¼ pffiffiffiffiffiffi eðjxyjÞ =2s : s 2p
ð7Þ
This is the integro-difference equivalent to diffusion and represents the situation where long-range dispersal in a short time period is extremely rare. See Fig. 1 (a) and (b) for a graph of this function. (II) A leptokurtic distribution (Zwillinger & Kokoska, 2000): pffiffiffiffiffiffiffiffi l2 l jxyj ; ð8Þ kðjx yjÞ ¼ e 4 where l=1/s. The above distribution is sometimes referred to as ‘‘fat-tailed’’ and long-range dispersal is relatively much more common than with the Gaussian distribution. See Fig. (1c) and (d). Kot et al. (1996) showed that a leptokurtic distribution function can produce waves of
increasing speed when there is no Allee effect. Lewis and Kareiva investigated the diffusion model with an Allee effect (Lewis & Kareiva, 1993), showing that dispersal can reduce local concentrations below the threshold required for invasion, even when the initial inoculum (number of infected individuals) is above the threshold. The initial size of the inoculum can therefore affect the ultimate success of invasion and I investigate the dynamic with various initial inoculum areas. I will present results for both increasing loss of maternal transmission and for the effect of initial incoculum size on the ultimate success of the invading wave. Model Results In this section and the remainder of the paper, I follow the conventions of Turelli and Hoffmann and formulate following definitions: K The asymptotic travelling wave speed is defined to be the rate at which the point with 50% of infected individuals moves through space once the wave has achieved its final shape. K The wave width w=Dx is defined to be the distance between the location on the wave front where the proportion of infected individuals in the population is 5% and the location where it has risen to 95%.
The numerical simulations were performed on a spatial grid divided into 200 grid points, the grid point spacing representing 10 km, thus covering a total domain of 2000 km. Numerical results were also produced on a grid of 1000 points with the same 10 km spacing to eliminate possible boundary effects, but the results of the simulation on the two grids were not different in the time frame examined in this paper. The reaction– diffusion equation was solved numerically with a standard central difference solver (Mitchell & Griffiths, 1980) using a finite difference discretization with zero flux boundary conditions and a time increment of 102 generations. Numerical simulations were performed for two different initial set conditions (see Fig. 2 for details): (I) a semi-infinite infection inoculum of density 90% for 0pxp500 km and zero elsewhere;
SPATIALLY EXPLICIT MODELS
125
Fig. 2. (a) Semi-infinite conditions (I) used in numerical simulations. (b) Finite width initial conditions (II) used in numerical simulations.
(II) an infection inoculum of density 90% centred at x=500 km with width L=60, 120 or 200 km and zero elsewhere. For all the following results, values for the effects of infection on hatch rate and fecundity were taken from Turelli and Hoffmann’s paper i.e. sh=0.45, sf=0.05 (Turelli & Hoffmann, 1991). The three models, reaction–diffusion, Gaussian and leptokurtic integro-difference models, were solved numerically using the same parameters and initial conditions. The results were also compared to the asymptotic (hyperbolic tangent) solution. The results using the semiinfinite initial condition (I) and perfect maternal transmission m=0 are shown in Fig. 3. Wave speeds and widths w agree with previous results for Gaussian and diffusion models using the same value of s=50 km gen1/2 (Turelli & Hoffmann, 1991). The wave speed was calculated for the asymptotic model as 12.5 km gen1 [Fig. 3(a–c)], the diffusion model as 12.5 km gen1 [Fig. 3(d–f )], the Gaussian kernel model as 12.5 km gen1 [Fig. 3(g–i)] and the leptokurtic kernel model as 40 km gen1 (not shown in the figure). Thus, the wave speed predicted by the leptokurtic model is very different from that predicted by the other models using the same dispersal parameter s. However, using a leptokurtic dispersal kernel with a smaller value of s=5, it is possible to obtain
the observed wave speed of 12.5 km gen1 [Fig. 3( j)–(l)]. The wave widths were calculated for the asymptotic model as 430 km, the diffusion model as 450 km and the Gaussian model as 440 km. Next the effects of infidelity of maternal transmission were investigated, again using semi-infinite initial conditions (I), but this time with a fixed s and with varying infidelity of maternal transmission rate m>0. The results are shown in Table 1. Columns 3 and 4 of the table show results for the diffusion and Gaussian kernel models with s=30, 60 and 90 km gen1/2 and column 5 shows comparable results for the leptokurtic kernel model achieved with s=3, 6 and 9 km gen1/2. The results for the diffusion model and the Gaussian kernel model are of the same order of magnitude but similar results are only obtained from the leptokurtic kernel model by using a dispersal rate s which is an order of magnitude smaller. As can be seen from the results, a loss rate m of 9% is sufficient to prevent the wave from advancing, and this is independent of the model dispersal mechanism used. It is not possible to derive a wave speed in this case. Travelling wave profiles for the reaction– diffusion (top), Gaussian (middle) and leptokurtic (bottom) models are shown in Fig. 4, with s=50 km gen1/2 in the reaction–diffusion and Gaussian models, and s=5 km gen1/2 in the
126
P. SCHOFIELD
Fig. 3. Plots of travelling wave solutions for all models with varying dispersal rate s at times t=15, 30, 45, 60 generations. (a)–(c) show results from the analytical hyperbolic tangent solution; (d)–(f ) show results from the reaction– diffusion model; (g)–(i) show the results from the Gaussian kernel model; ( j)–(l) show the results from the leptokurtic kernel model. Parameter values are s=30 km gen1/2 in (a), (d), (g); s=60 km gen1/2 in (b), (e), (h); s=90 km gen1/2 in (c), (f ), (i); s=3 km gen1/2 in ( j); s=6 km gen1/2 in (k); s=9 km gen1/2 in (1). ( ) 15; ( ) 30; ( ) 45; ( ) 60.
leptokurtic model. As can be seen from the results in Fig. 4, even a relatively small loss in maternal transmission can have a big effect on the ultimate success of the invasion. Fig. 4(b), (e) and (h) all have m=5% and invasion occurs in each case whereas Fig. 4(c), (f) and (i) have m=9% and invasion fails in each case. Where invasion does succeed the infection level does not
reach fixation, but rather a coexistent equilibrium is reached; that is, a balance between the increase due to the advantage of infection in a population, where infection is above the threshold level, and the production of uninfected individuals due to imperfect transmission. To complement the numerical results, a local stability analysis was performed on the spatially
SPATIALLY EXPLICIT MODELS
Table 1 Wave speeds in km gen1 for various dispersal rates s and transmission loss rate m s m 0.00 0.05 0.09 0.00 60 0.05 0.09 0.00 90 0.05 0.09 30
Diffusion Gaussian Leptokurtic m wave speed wave speed wave speed 7.8 3.6 * 15.0 7.7 * 22.5 12.9 *
7.5 3.1 * 15.0 6.7 * 22.5 10 *
6.0 2.7 * 14.0 8.3 * 22.5 14.0 *
s
0.00 0.05 3 0.09 0.00 0.05 6 0.09 0.00 0.05 9 0.09
Note: Column 3 shows wave speeds calculated from the reaction–diffusion model, column 4 shows wave speeds calculated from the Gaussian kernel model, column 5 shows wave speeds calculated from the leptokurtic kernel model. *indicates where it is not possible to derive a +ve wave speed as the threshold for m is exceeded.
127
homogeneous solution of eqn (3). This produces three equilibria: p1 ¼ 0;
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2f 2sf sh þ s2h 4mðsh þ sh sf Þ p2 ¼ ; 2sh qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðsf þ sh Þ þ s2f 2sf sh þ s2h 4mðsh þ sh sf Þ p3 ¼ ; 2sh ðsf þ sh Þ
ð9Þ where p1 is the stable no infection equilibrium, p3 is the stable coexistence equilibrium and p2 is the unstable threshold. From this it can be seen that if mXðs2f 2sf sh þ s2h Þ=ð4sh ð1 sf ÞÞ then the threshold disappears and there is only a single stable equilibrium at p=0 and infection cannot
Fig. 4. Plots of travelling wave solutions for the three model systems with varying transmission rates m at times t=15, 30, 45, 60 generations. (a)–(c) show the results from the reaction–diffusion model with s=50 km gen1/2; (d)–(f ) show the results from the Gaussian kernel model with s=50 km gen1/2; (g)–(i) show the results from the leptokurtic kernel model with s=5 km gen1/2. Parameter values are m=0 (a), (d) and (g); m=0.05 (b), (e) and (h); m=0.09 (c), (f) and (i). ( ) 15; ( ) 30; ( ) 45; ( ) 60.
128
P. SCHOFIELD
invade. For the values sh=0.45, sf=0.05 this requires mo9% for invasion, which corresponds to the result found from numerical simulations. Figure 5 shows the bifurcation diagram of change of equilibrium and threshold with changing m. As m approaches 9%, the threshold rises and the stable equilibrium value decreases, the coexistence equilibrium vanishes and there is only a single equilibrium where infection is removed from the population. Turelli & Hoffmann (1995) estimate mX4.7% and this could still be within the range of m where infection can invade the population. However, as can be seen in Table 1, a value of m=5% reduces the wave speed to less than half that in a population with perfect transmission. Finally, the effect of the width L of a local, spatially finite inoculum, was investigated. Infected individuals dispersing from a local inoculum and being replaced by uninfected individuals from the area surrounding the inoculum will reduce the proportion of infected individuals in the inoculum. If the dispersal is rapid enough, or the advantage of infection is not large enough, the proportion infected will fall below the threshold level before infected numbers can build up in the surrounding areas and invasion will not be possible. To investigate this, the diffusion, Gaussian kernel and leptokurtic kernel models were solved numerically using initial condition (II). A dispersal rate s=50 km gen1/2 was used in the diffusion and
Fig. 5. Bifurcation diagram showing the change in local stable coexistence equilibrium S and unstable threshold U with increasing infidelity of maternal transmission m. ( )S; ( )U.
Gaussian kernel models, and a dispersal rate s=5 km gen1/2 was used in the leptokurtic model. The width L of the initial inoculum and the rate transmission loss m were varied for each simulation run. The transient dynamics produced by these simulations are shown in Fig. 6. The three rows (top, middle and bottom) correspond to the results from the three different model types. Row one [Fig. 6(a–c)] gives the results of the reaction– diffusion equation (RDE) model, row two [Fig. 6(e–f)] the results of the Gaussian kernel integro-difference (GKI) model and row three [Fig. 6(g–i)] the results from the leptokurtic kernel integro-difference (LKI) model. Each column of Fig. 6 represents a different initial inoculum width L. The results shown in column 1 [Fig. 6(a),(d) and (g)] are for L=60 km, column 2 [Fig. 6(b),(e) and (h)] for L=120 km and column 3 [Fig. 6(c), (f ) and (i)] for L=200 km. With an initial inoculum width of L=60 km it can be seen that the RDE model predicts that infection cannot invade even with perfect transmission m=0 [Fig. 6(a)]. The GKI model predicts invasion only with perfect transmission [Fig. 6(d)], but the LKI model predicts invasion with loss rates of up to 4% (m=0.04) [Fig. 6(e)]. Increasing the initial inoculum width to L=120 km, the RDE model predicts invasion with mp0.02 (2% loss) [Fig. 6(b)]. The GKI model predicts invasion with perfect transmission [Fig. 6(e)] while the LKI model predicts invasion with mp0.06 [Fig. 6(h)]. Finally, increasing the width to L=200 km the RDE model predicts invasion with m=0.04 [Fig. 6(c)]. The GKI model predicts invasion mp0.06 (Fig. 6(f ) and the LKI model still predicts invasion with mp0.006 [Fig. 6(i)]. This shows that increasing the size of the inoculum increases the ultimate success of infection spread in all models. In all the figures, it can be seen that increasing the transmission loss rate m reduces the equilibrium level of infection when invasion is successful, although as can be seen, this is not affected by either L or the model chosen but depends only on m. In summary, the predicted long-term spread of infection in the population is dependent not only the parameters m and s but also on the assumptions
129
SPATIALLY EXPLICIT MODELS
Fig. 6. Plots of the level of infection I at x=500 km for times 0ptp150 generations using different initial inoculum widths L. (a)–(c) show results from the reaction–diffusion model with s=40 km gen1/2; (d)–(f) show results from the Gaussian kernel model with s=40 km gen1/2; (g)–(i) show results from the leptokurtic kernel model with s=4 km gen1/2. (a), (d) and (g) have an initial inoculum width of L=60 km; (b), (e) and (h) have L=120 km; (c), (f) and (i) and L=200 km. ( ) 0.00; ( ) 0.02; ( ) 0.04; ( ) 0.06; ( ) 0.08.
made about initial conditions, and the initial (spatial) size of the infected population. Thus we can say: K An increase in the initial width of the inoculum L increases the predicted ultimate success of the invasion of infection. K An increase in the loss rate of maternal transmission m decreases the predicted ultimate success of the invasion of infection. K An increase in the loss rate of maternal transmission m reduces the predicated eventual equilibrium level of infection after successful invasion.
Additionally, unlike the simulation results when semi-infinite initial condition (I) was used, significant variation is evident between the results of the three models when the finite initial condition (II) was used. This indicates that
assumptions about the dispersal behaviour (probability distribution functions) of the insects and the precise modelling technique chosen also affect the prediction success of invasion. The reaction–diffusion model requires a larger initial inoculum width L than the Gaussian kernel model to achieve invasion, and in turn, the Gaussian kernel model requires a larger initial inoculum width L than the leptokurtic kernel model. Discussion This paper has examined the variation in results obtained from simulations of different mathematical models of infection spread that use varying explicit assumptions about insect dispersal patterns. Some comparison is made between these results and earlier results and field
130
P. SCHOFIELD
data. The model of Turelli & Hoffmann (1991) is extended with the inclusion of a maternal loss rate which raises the threshold for invasion. This is shown to slow down and eventually reverse the invasion wave with increasing rate of loss. This would confirm further the unsatisfactory nature of the diffusion model. Some analytical results for various classes of dispersal kernel with birth death kinetics the include a simple Allee effect have been investigated (Kot et al., 1996). They have also shown that modelling the movement of insects with the diffusion model or the equivalent Gaussian kernel, integro-difference equation model may not be appropriate, leading to significant underestimation of wave speeds. The results presented in this paper on the whole agree with this conclusion and with that of Turelli and Hoffmann, that the actual spread of Wolbachia in the California D. simulans population is not best modelled with a diffusion-based model. More generally, the results show that the transient dynamics of invasion subject to a threshold similar to an Allee effect are sensitive to the modelling scheme chosen, and it is possible that asymptotic solutions for the reaction–diffusion and Gaussian dispersal models may not provide robust results. The results also tend to confirm the views of other research (Minogue, 1989; Shaw, 1994) that for some populations a better model of redistribution is achieved using ‘‘fat-tailed’’ non-exponentially bounded kernels. It is possible that this might be the case for the California D. simulans data (Turelli & Hoffmann, 1991)Fthe leptokurtic kernel giving faster wave speeds for the same value of s than the Gaussian. With the model chosen here, this was enough to account for the discrepancy between the observed movement of Drosophila and the estimated dispersal rate s (Turelli & Hoffmann, 1991). However, it is possible that this faster wave speed could be produced by an additional convection term in the reaction–diffusion model that reflected some external pressure, such as human activity transporting fruit or abiotic factors such as prevailing winds. This study is also in keeping with other studies of insect dispersal systems (Hasting, 2000) which suggest that non-Gaussian dispersal kernels may better fit the empirical data.
In showing that the wave speed is sensitive to the dispersal kernel chosen in the integrodifference model, this would suggest that arbitrary choices of dispersal kernel based purely on ease of mathematical analysis may miss important aspects of the system being modelled. This may be particularly important when it comes to the success of invasion from a finite inoculum into a population where rapid dispersal may be disadvantageous to the ultimate success of the invasion Kot et al. (1996) reanalysis of Drosophila pseudoobscura dispersal data from Dobzhansky & Wright (1947) showed that the data were not sufficient to distinguish between an exponentially bounded and a non-exponentially bounded dispersal function. The conclusion was made that better data concerning the probability of rare extremely long distance dispersal are needed to distinguish between the two functions. Attempts to discover more accurate dispersal probability distributions and parameters from field data have been made. A series of papers on mark recapture experiments with Drosophila sp. (Jones et al., 1981; Coyne et al., 1982, 1987) indicated that individuals released in desert conditions can be found to have travelled much further from the release site in a short period of time than diffusion estimates would predict. These were studies of organisms dispersing in non-native environments, and, as pointed out in the Introduction with the possible problems of using Drosophila pseudoobscura parameter estimates in a model of Drosophila simulans, estimates derived for species in non-native environment may not be adequate to model the species in its natural habitat (Lewis, 1997). I am very grateful to Michael Turelli for his guidance and suggestions during this work and to the University of California Davis for hosting me. I am also grateful for the advice of Professor Chaplain and the comments of two anonymous referees in helping revise the manuscript. This work has been funded in part by the BBSRC and by a special travel grant from the Scottish International Education Trust. REFERENCES Barton, N. H. (1979). The dynamics of hybrid zones. Heredity 43, 341–359. Breeuwer, J. & Werren, J. H. (1990). Microorganisms associated with chromosome destruction and
SPATIALLY EXPLICIT MODELS
reproductive isolation between two insect species. Nature 346, 558–560. Coyne, J. A., Boussy, I. A., Prout, T., Bryant, S. H., Jones, J. & Moore, J. A. (1982). Long distance migration in Drosophila. Am. Nat. 119, 589–595. Coyne, J. A., Bryant, S. H. & Turelli, M. (1987). Long distance migration of Drosophila 2. Presence in desolate sites and dispersal near desert oasis. Am. Nat. 129, 847–861. Dobzhansky, T. & Wright, S. (1947). Genetics in natural population. XV. Rate of difussion of a mutant gene through a population of Drosophila pseudoobscura. Genetics 32, 303–324. Dunbar, S. R. (1983). Travelling wave solution of diffusive Lotka–Volterra equations. J. Math. Biol. 17, 11–32. Fisher, R. (1937). The wave of advance of an advantageous gene. Ann. Eugenics 7, 353–369. Ghelelovitch, S. (1952). Sur le determinime genetique de sterilite dans le croisement entre differentes souch de Culex autogenicus Roubaud. C. R. Acad. Sci. Paris 24, 2386–2388. Hasting, A. (2000). Parasitoid spread: lessons for and from invasion biology. In: Parasitoid Population Biology (Hochberg, M. E. & Ives, A. R., eds), pp. 70–82. Princeton, NJ: Princeton Unviersity Press. Hertig, M. & Wolbach, S. (1924). Studies in rickettsialike microorganism in insects. J. Med. Res. 44, 329–374. Hoffmann, A. A., Turelli, M. & Simmons, G. (1986). Unidirectional incompatibility between populations of Drosophila simulans. Evolution 40, 692–701. Jones, J., Bryant, S. H., Lewontin, R., Moore, J. A. & Prout, T. (1981). Gene flow and the geographicaldistribution of a molecular polymorphism in Drosophilapseudoobscura. Genetics 98, 157–178. Kimura, M. (1953). ‘Stepping stone’ model of population. Annual Report of the National Institute of Genetics, Japan, Vol. 3, 62–63. Kot, M. (1992). Discrete-time travelling waves: ecological examples. J. Math. Biol. 30, 413–436. Kot, M. & Schaffer, W. M. (1986). Discrete-time growth-dispersal models. Math. Biosci. 80, 109–136. Kot, M., Lewis, M. A. & van den Driessche, P. (1996). Dispersal data and the spread of invading organisms. Ecology 77, 2027–2042. Lewis, M. A. (1997). Spread if an invading population. In: Spatial Ecology (Tilman, D. & Karevia, P. M., eds),
131
pp. 46–69. Monographs in Population Biology, No. 30. Princeton, NJ: Princeton University Press. Lewis, M. A. & Kareiva, P. M. (1993). Allee dynamics and the spread of invading organisms. Theoret. Popul. Biol. 43, 141–158. Minogue, K. (1989). Diffusion and spatial probability models for disease. In: Spatial Components of Plat Disease Epidemics (Jeger, M., ed), pp. 127–143. Englewood Cliffs, NJ: Prentice-Hall. Mitchell, A. & Griffiths, D. (1980). The Finite Difference Method in Partial Differential Equations. New York: John Wiley & Sons. Mollison, D. (1977). Spatial contact models for ecological and epidemic spread. J. R. Stat. Soc. B 39, 283–326. Murray, J. (1989). Mathematical Biology. Biomathematics Texts, No. 19. Berlin: Springer-Verlag. Nagylaki, T. (1975). Conditions for the existence of clines. Genetics 80, 595–615. Okubo, A. (1980). Diffusion and Ecological Problems: Mathematical Models. Berlin: Springer-Verlag. O’neill, S. L., Hoffmann, A. A. & Werren, J. H. (eds) (1997). Influential Passangers Inherited Microorganisms and Arthropod Reproduction. Oxford: Oxford University Press. Renshaw, E. (1991). Modelling Biological Populations in Space and Time, Cambridge Studies in Mathematical Biology. Cambridge: Cambridge University Press. Shaw, M. (1994). Modeling stochastic processes in plant pathology. Annu. Rev. Phytopathol. 32, 523 –544. Turchin, P. (1998). Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants. Sinauer Associates. Turelli, M. & Hoffmann, A. A. (1991). Rapid spread of an inherited incompatability factor on California Drosophila. Nature 353, 440 – 442. Turelli, M. & Hoffmann, A. A. (1995). Cytoplasmic incompatibility in Drosophila simulans: dynamics and parameter estimates from natural populations. Genetics 140, 1319–1338. Turelli, M., Hoffmann, A. A. & Mckechnie, S. W. (1992). Dynamics of cytoplasmic incompatability and mtDNA variation in natural Drosophila simulans populations. Genetics 132, 713–723. Werren, J. H. (1997). Biology of Wolbachia. Annu. Rev. Entomol. 42, 587–609. Zwillinger, D. & Kokoska, S. (2000). CRC Standard Probability and Statistics Tables and Formulae. London/ Boca Raton, FL: Chapman & Hall/CRC.
Spatially Explicit Models of Turelli–Hoffmann Wolbachia Invasive Wave Fronts Peter Schofield Department of Mathematics, University of Dundee, DD1 4HN, U.K. (Received on 5 January 2001, Accepted in revised form on 2 December 2001)
This paper examines different mathematical models of insect dispersal and infection spread and compares these with field data. Reaction–diffusion and integro-difference equation models are used to model the spatio-temporal spread of Wolbachia in Drosophila simulans populations. The models include cytoplasmic incompatibility between infected females and uninfected males that creates a threshold density, similar to an Allee effect, preventing increase from low incidence of infection in the host population. The model builds on an earlier model (Turelli & Hoffmann, 1991) by incorporating imperfect maternal transmission. The results of simulations of the models using the same parameter values produce different dynamics for each model. These differences become very marked in the integro-difference equation models when insect dispersal patterns are assumed to be non-Gaussian. The success or failure of invasion by Wolbachia in the simulations may be attributed to the insect dispersal mechanism used in the model rather than the parameter values. As the models predict very different outcomes for the integro-difference models depending on the underlying assumptions of insect dispersal patterns, this emphasizes that good field data on real (rather than idealized) dispersal patterns need to be collected before models such as these can be used for predictive purposes. r 2002 Elsevier Science Ltd.
Introduction The intra-cellular organism Wolbachia is a much studied reproductive parasite. Primarily vertically transmitted in egg cytoplasm, it manipulates its host’s reproduction in ways that increase the proportion (or fitness) of infected females in the population and hence the proportion of infected offspring in the next generation. The details vary depending on host species and Wolbachia strain (Werren, 1997; O’Neill et al., 1997). Turelli and Hoffmann studied the sweep of Wolbachia through the Drosophila simulans population in California in a series of papers (Hoffmann et al., 1986; Turelli & Hoffmann, 1991, 1995; Turelli et al., 1992). The Wolbachia 0022-5193/02/050121 + 11 $35.00/0
strain studied caused a reproduction distorting effect referred to as cytoplasmic incompatibility (CI). First observed in mosquitoes (Hertig & Wolbach, 1924) and attributed to infection (Ghelelovitch, 1952), CI has the effect that crosses between infected males and uninfected females suffer an increased hatch rate failure. It is thought that Wolbachia imprints the sperm of infected males in a manner that requires infection in the egg to recover the sperm function (Breeuwer & Werren, 1990). Turelli showed that there is a threshold level in the dynamics of invasion by the infected phenotype into an uninfected population (Turelli & Hoffmann, 1995) which produces local dynamics similar to an Allee effect, where the unstable equilibrium r 2002 Elsevier Science Ltd.
122
P. SCHOFIELD
threshold must be exceeded for invasion to be successful. There are also two stable equilibria with either no infection or infection fixation. Imperfect maternal transmission has two effects: (1) it produces a stable equilibrium where uninfected individuals are maintained in the population; (2) it increases the level of the unstable equilibrium or threshold for invasion (Turelli & Hoffmann, 1995). In these models, we assume that the level of horizontal transmission of infection is negligible (Turelli et al., 1992). One of the main motivations for the current investigation is the observation of a wave-like advance of Wolbachia infection through the California D. simulans populations during the late 1980s (Turelli & Hoffmann, 1991). In their paper, Turelli and Hoffmann suggest that in order for the reaction–diffusion equation model’s results to fit the data, the D. simulans dispersal rate s must be of the order of 50 km gen1/2. They comment that this is almost an order of magnitude greater than the observed dispersal rates in mark-recapture experiments (Jones et al., 1981) and this might indicate that some dispersal process other than diffusion is acting. Simple diffusion is therefore unlikely to be the correct dispersal mechanism to use. However, it must be noted that that mark-recapture data in Jones et al. (1981) is for D. pseudoobscura rather than D. simulans, and therefore the lack of fit may be due to species differences rather than the diffusion asumption. It has also been shown that a convection process, such as prevailing wind, can increase the wave speed but will not otherwise alter the shape of the wave front (Murray, 1989; Lewis & Kareiva, 1983). Another approach that has been employed recently in modelling insect dispersal (in particular Drosophila sp. movements) is that of integro-difference equations (Kot et al., 1996). The work of Kot et al. (1996) reanalyses Drosophila pseudoobscura dispersal data from Dobzhansky & Wright (1947), showing that the data were not sufficient to distinguish between exponentially bounded and non-exponentially bounded dispersal functions. In an exponentially bounded distribution, the probability of movement from the initial location to a new location (within a fixed time) approached zero exponentially or faster with the increase in distance
between the two locations. With the nonexponentially bounded distribution, this is not the case. Long-range redistribution events have a finite probability, see Fig. 1. The non-exponentially bounded case could represent a situation where most movement is extremely local, for example, insects with weak flight, but there is a significant probability of much longer range movement, perhaps by human activity such as fruit transportation. Kot et al. (1996) also investigated the Allee effect showing that this slows down the effective wave speed. This has motivated the present examination of the California Wolbachia sweep using integro-difference models with contrasting dispersal kernels. Mathematical Models Models of invasion using Fisher-type reaction diffusion equations (Fisher, 1937; Murray, 1989), integro-difference equations (Mollison, 1977) and stepping stone models (Kimura, 1953; Renshaw, 1991), that produce a travelling wave of advance have been well studied. Analytical results have been obtained for diffusion equation models of the advance of advantageous genes (Fisher, 1937; Murray, 1989), movement in hybrid zones (Nagylaki, 1975; Barton, 1979) and invasion in predator–prey systems (Dunbar, 1983). The dynamics of travelling waves from a Fisher-type equation that included an Allee effect have also been investigated (Lewis & Kareiva, 1993). The reproduction distortion mechanism of CI, studied in Hoffmann and Turelli’s work on D. simulans, also produces an effect similar to an Allee effect. This means that if infected males are at a high frequency in the population, then uninfected females are at a disadvantage. If the infection is above a critical frequency then infection will increase to fixation, that is, all members of the population are infected. If there is imperfect maternal transmission fixation is not reached, but rather an equilibrium balance between fitness increase due to infection and infection loss is achieved. Turelli & Hoffmann (1991) assumed perfect maternal transmission and a fecundity penalty for infected females. This produced the following equation for the change Dp in the proportion of infected
SPATIALLY EXPLICIT MODELS
123
Fig. 1. Subplots (a) and (b) show the graph of the Gaussian distribution function with dispersal rate s=50 km gen1/2. Subplots (c) and (d) show the graph of the leptokurtic distribution function with dispersal rate s=5 km gen1/2. Comparison of the graphs (b) and (d) shows that the Gaussian function tends to zero faster than the leptokurtic function.
individuals p in the population between each generation: Dp ¼
# sh pð1 pÞðp pÞ ; 1 sf p sh pq
ð1Þ
where p# ¼ sf =sh ; sf is the reduction in fecundity caused by Wolbachia infection and sh is the reduction in hatch success in incompatible crosses (infected male cross uninfected female). Turelli and Hoffmann make the approximation that the terms s f p and sh pq are negligible and hence mean fitness can be assumed to be 1. The addition of spatial movement by a diffusion mechanism gives the following equation for the spatial spread of infected individuals: @p @2 p # ¼ D 2 þ sh pð1 pÞðp pÞ; @t @x
ð2Þ
where D is the diffusion coefficient. It is assumed that D=s2/2 where s is the variance of the individual movement probability distribution (Okubo, 1980; Turchin, 1998).
Turelli and Hoffmann use the analytical result that in the limit t-N (i.e. time tends to infinity) eqn (2) produces travelling wave fronts that can be approximated by a hyperbolic tangent function (Barton, 1979). This models the asymptotic constant velocity wave of advance of infection through the population. Here I examine a modified model that extends eqn (2) to include an imperfect maternal transmission rate m producing the equation @p @2 p # mð1 sf Þp: ð3Þ ¼ D 2 þ sh pð1 pÞðp pÞ @t @x If m=0 then eqn (3) reduces to eqn (2). In the next section, I compare numerical solutions of the reaction–diffusion eqn (3) with the solution using an asymptotic approximation of the equation, similar to Barton (1979), when m=0, given by " pffiffiffiffi !# # sh t=2 x sð1 2pÞ 1 pffiffiffiffi : pðx; tÞ ¼ 1 tanh 2 2s= sh ð4Þ
124
P. SCHOFIELD
These (continuous) reaction–diffusion models can be compared to other models that use a different explicit spatial modelling approach, that of discrete-time, integro-difference equation models (Kot & Schaffer, 1986; Kot, 1992), where movement is modelled by dispersal at the end of each discrete generation using a probability density function. Integro-difference models are based on equations of the generic form ptþ1 ðxÞ ¼
Z
b
kðjx yjÞ f ð pt ð yÞÞ dy;
ð5Þ
a
where k(|xy|) is the dispersal kernel and f( pt(y)) a function representing the local kinetics at spatial location y. The particular form I derived using the kinetics of eqn (3) is as follows: ptþ1 ðxÞ ¼
Z
b
kðjx yjÞ½sh pt ð yÞð1 pt ð yÞÞ a
# mð1 sf Þpt ð yÞ dy: ð6Þ ð pt ð yÞ pÞ It has been shown that the shape of the dispersal kernel can have significant effects on the resulting wave of invasion (Kot et al., 1996). Using eqn (6), I compare results for two probability density functions: (I) The Gaussian (normal ) distribution function: 2 1 2 kðjx yjÞ ¼ pffiffiffiffiffiffi eðjxyjÞ =2s : s 2p
ð7Þ
This is the integro-difference equivalent to diffusion and represents the situation where long-range dispersal in a short time period is extremely rare. See Fig. 1 (a) and (b) for a graph of this function. (II) A leptokurtic distribution (Zwillinger & Kokoska, 2000): pffiffiffiffiffiffiffiffi l2 l jxyj ; ð8Þ kðjx yjÞ ¼ e 4 where l=1/s. The above distribution is sometimes referred to as ‘‘fat-tailed’’ and long-range dispersal is relatively much more common than with the Gaussian distribution. See Fig. (1c) and (d). Kot et al. (1996) showed that a leptokurtic distribution function can produce waves of
increasing speed when there is no Allee effect. Lewis and Kareiva investigated the diffusion model with an Allee effect (Lewis & Kareiva, 1993), showing that dispersal can reduce local concentrations below the threshold required for invasion, even when the initial inoculum (number of infected individuals) is above the threshold. The initial size of the inoculum can therefore affect the ultimate success of invasion and I investigate the dynamic with various initial inoculum areas. I will present results for both increasing loss of maternal transmission and for the effect of initial incoculum size on the ultimate success of the invading wave. Model Results In this section and the remainder of the paper, I follow the conventions of Turelli and Hoffmann and formulate following definitions: K The asymptotic travelling wave speed is defined to be the rate at which the point with 50% of infected individuals moves through space once the wave has achieved its final shape. K The wave width w=Dx is defined to be the distance between the location on the wave front where the proportion of infected individuals in the population is 5% and the location where it has risen to 95%.
The numerical simulations were performed on a spatial grid divided into 200 grid points, the grid point spacing representing 10 km, thus covering a total domain of 2000 km. Numerical results were also produced on a grid of 1000 points with the same 10 km spacing to eliminate possible boundary effects, but the results of the simulation on the two grids were not different in the time frame examined in this paper. The reaction– diffusion equation was solved numerically with a standard central difference solver (Mitchell & Griffiths, 1980) using a finite difference discretization with zero flux boundary conditions and a time increment of 102 generations. Numerical simulations were performed for two different initial set conditions (see Fig. 2 for details): (I) a semi-infinite infection inoculum of density 90% for 0pxp500 km and zero elsewhere;
SPATIALLY EXPLICIT MODELS
125
Fig. 2. (a) Semi-infinite conditions (I) used in numerical simulations. (b) Finite width initial conditions (II) used in numerical simulations.
(II) an infection inoculum of density 90% centred at x=500 km with width L=60, 120 or 200 km and zero elsewhere. For all the following results, values for the effects of infection on hatch rate and fecundity were taken from Turelli and Hoffmann’s paper i.e. sh=0.45, sf=0.05 (Turelli & Hoffmann, 1991). The three models, reaction–diffusion, Gaussian and leptokurtic integro-difference models, were solved numerically using the same parameters and initial conditions. The results were also compared to the asymptotic (hyperbolic tangent) solution. The results using the semiinfinite initial condition (I) and perfect maternal transmission m=0 are shown in Fig. 3. Wave speeds and widths w agree with previous results for Gaussian and diffusion models using the same value of s=50 km gen1/2 (Turelli & Hoffmann, 1991). The wave speed was calculated for the asymptotic model as 12.5 km gen1 [Fig. 3(a–c)], the diffusion model as 12.5 km gen1 [Fig. 3(d–f )], the Gaussian kernel model as 12.5 km gen1 [Fig. 3(g–i)] and the leptokurtic kernel model as 40 km gen1 (not shown in the figure). Thus, the wave speed predicted by the leptokurtic model is very different from that predicted by the other models using the same dispersal parameter s. However, using a leptokurtic dispersal kernel with a smaller value of s=5, it is possible to obtain
the observed wave speed of 12.5 km gen1 [Fig. 3( j)–(l)]. The wave widths were calculated for the asymptotic model as 430 km, the diffusion model as 450 km and the Gaussian model as 440 km. Next the effects of infidelity of maternal transmission were investigated, again using semi-infinite initial conditions (I), but this time with a fixed s and with varying infidelity of maternal transmission rate m>0. The results are shown in Table 1. Columns 3 and 4 of the table show results for the diffusion and Gaussian kernel models with s=30, 60 and 90 km gen1/2 and column 5 shows comparable results for the leptokurtic kernel model achieved with s=3, 6 and 9 km gen1/2. The results for the diffusion model and the Gaussian kernel model are of the same order of magnitude but similar results are only obtained from the leptokurtic kernel model by using a dispersal rate s which is an order of magnitude smaller. As can be seen from the results, a loss rate m of 9% is sufficient to prevent the wave from advancing, and this is independent of the model dispersal mechanism used. It is not possible to derive a wave speed in this case. Travelling wave profiles for the reaction– diffusion (top), Gaussian (middle) and leptokurtic (bottom) models are shown in Fig. 4, with s=50 km gen1/2 in the reaction–diffusion and Gaussian models, and s=5 km gen1/2 in the
126
P. SCHOFIELD
Fig. 3. Plots of travelling wave solutions for all models with varying dispersal rate s at times t=15, 30, 45, 60 generations. (a)–(c) show results from the analytical hyperbolic tangent solution; (d)–(f ) show results from the reaction– diffusion model; (g)–(i) show the results from the Gaussian kernel model; ( j)–(l) show the results from the leptokurtic kernel model. Parameter values are s=30 km gen1/2 in (a), (d), (g); s=60 km gen1/2 in (b), (e), (h); s=90 km gen1/2 in (c), (f ), (i); s=3 km gen1/2 in ( j); s=6 km gen1/2 in (k); s=9 km gen1/2 in (1). ( ) 15; ( ) 30; ( ) 45; ( ) 60.
leptokurtic model. As can be seen from the results in Fig. 4, even a relatively small loss in maternal transmission can have a big effect on the ultimate success of the invasion. Fig. 4(b), (e) and (h) all have m=5% and invasion occurs in each case whereas Fig. 4(c), (f) and (i) have m=9% and invasion fails in each case. Where invasion does succeed the infection level does not
reach fixation, but rather a coexistent equilibrium is reached; that is, a balance between the increase due to the advantage of infection in a population, where infection is above the threshold level, and the production of uninfected individuals due to imperfect transmission. To complement the numerical results, a local stability analysis was performed on the spatially
SPATIALLY EXPLICIT MODELS
Table 1 Wave speeds in km gen1 for various dispersal rates s and transmission loss rate m s m 0.00 0.05 0.09 0.00 60 0.05 0.09 0.00 90 0.05 0.09 30
Diffusion Gaussian Leptokurtic m wave speed wave speed wave speed 7.8 3.6 * 15.0 7.7 * 22.5 12.9 *
7.5 3.1 * 15.0 6.7 * 22.5 10 *
6.0 2.7 * 14.0 8.3 * 22.5 14.0 *
s
0.00 0.05 3 0.09 0.00 0.05 6 0.09 0.00 0.05 9 0.09
Note: Column 3 shows wave speeds calculated from the reaction–diffusion model, column 4 shows wave speeds calculated from the Gaussian kernel model, column 5 shows wave speeds calculated from the leptokurtic kernel model. *indicates where it is not possible to derive a +ve wave speed as the threshold for m is exceeded.
127
homogeneous solution of eqn (3). This produces three equilibria: p1 ¼ 0;
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2f 2sf sh þ s2h 4mðsh þ sh sf Þ p2 ¼ ; 2sh qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðsf þ sh Þ þ s2f 2sf sh þ s2h 4mðsh þ sh sf Þ p3 ¼ ; 2sh ðsf þ sh Þ
ð9Þ where p1 is the stable no infection equilibrium, p3 is the stable coexistence equilibrium and p2 is the unstable threshold. From this it can be seen that if mXðs2f 2sf sh þ s2h Þ=ð4sh ð1 sf ÞÞ then the threshold disappears and there is only a single stable equilibrium at p=0 and infection cannot
Fig. 4. Plots of travelling wave solutions for the three model systems with varying transmission rates m at times t=15, 30, 45, 60 generations. (a)–(c) show the results from the reaction–diffusion model with s=50 km gen1/2; (d)–(f ) show the results from the Gaussian kernel model with s=50 km gen1/2; (g)–(i) show the results from the leptokurtic kernel model with s=5 km gen1/2. Parameter values are m=0 (a), (d) and (g); m=0.05 (b), (e) and (h); m=0.09 (c), (f) and (i). ( ) 15; ( ) 30; ( ) 45; ( ) 60.
128
P. SCHOFIELD
invade. For the values sh=0.45, sf=0.05 this requires mo9% for invasion, which corresponds to the result found from numerical simulations. Figure 5 shows the bifurcation diagram of change of equilibrium and threshold with changing m. As m approaches 9%, the threshold rises and the stable equilibrium value decreases, the coexistence equilibrium vanishes and there is only a single equilibrium where infection is removed from the population. Turelli & Hoffmann (1995) estimate mX4.7% and this could still be within the range of m where infection can invade the population. However, as can be seen in Table 1, a value of m=5% reduces the wave speed to less than half that in a population with perfect transmission. Finally, the effect of the width L of a local, spatially finite inoculum, was investigated. Infected individuals dispersing from a local inoculum and being replaced by uninfected individuals from the area surrounding the inoculum will reduce the proportion of infected individuals in the inoculum. If the dispersal is rapid enough, or the advantage of infection is not large enough, the proportion infected will fall below the threshold level before infected numbers can build up in the surrounding areas and invasion will not be possible. To investigate this, the diffusion, Gaussian kernel and leptokurtic kernel models were solved numerically using initial condition (II). A dispersal rate s=50 km gen1/2 was used in the diffusion and
Fig. 5. Bifurcation diagram showing the change in local stable coexistence equilibrium S and unstable threshold U with increasing infidelity of maternal transmission m. ( )S; ( )U.
Gaussian kernel models, and a dispersal rate s=5 km gen1/2 was used in the leptokurtic model. The width L of the initial inoculum and the rate transmission loss m were varied for each simulation run. The transient dynamics produced by these simulations are shown in Fig. 6. The three rows (top, middle and bottom) correspond to the results from the three different model types. Row one [Fig. 6(a–c)] gives the results of the reaction– diffusion equation (RDE) model, row two [Fig. 6(e–f)] the results of the Gaussian kernel integro-difference (GKI) model and row three [Fig. 6(g–i)] the results from the leptokurtic kernel integro-difference (LKI) model. Each column of Fig. 6 represents a different initial inoculum width L. The results shown in column 1 [Fig. 6(a),(d) and (g)] are for L=60 km, column 2 [Fig. 6(b),(e) and (h)] for L=120 km and column 3 [Fig. 6(c), (f ) and (i)] for L=200 km. With an initial inoculum width of L=60 km it can be seen that the RDE model predicts that infection cannot invade even with perfect transmission m=0 [Fig. 6(a)]. The GKI model predicts invasion only with perfect transmission [Fig. 6(d)], but the LKI model predicts invasion with loss rates of up to 4% (m=0.04) [Fig. 6(e)]. Increasing the initial inoculum width to L=120 km, the RDE model predicts invasion with mp0.02 (2% loss) [Fig. 6(b)]. The GKI model predicts invasion with perfect transmission [Fig. 6(e)] while the LKI model predicts invasion with mp0.06 [Fig. 6(h)]. Finally, increasing the width to L=200 km the RDE model predicts invasion with m=0.04 [Fig. 6(c)]. The GKI model predicts invasion mp0.06 (Fig. 6(f ) and the LKI model still predicts invasion with mp0.006 [Fig. 6(i)]. This shows that increasing the size of the inoculum increases the ultimate success of infection spread in all models. In all the figures, it can be seen that increasing the transmission loss rate m reduces the equilibrium level of infection when invasion is successful, although as can be seen, this is not affected by either L or the model chosen but depends only on m. In summary, the predicted long-term spread of infection in the population is dependent not only the parameters m and s but also on the assumptions
129
SPATIALLY EXPLICIT MODELS
Fig. 6. Plots of the level of infection I at x=500 km for times 0ptp150 generations using different initial inoculum widths L. (a)–(c) show results from the reaction–diffusion model with s=40 km gen1/2; (d)–(f) show results from the Gaussian kernel model with s=40 km gen1/2; (g)–(i) show results from the leptokurtic kernel model with s=4 km gen1/2. (a), (d) and (g) have an initial inoculum width of L=60 km; (b), (e) and (h) have L=120 km; (c), (f) and (i) and L=200 km. ( ) 0.00; ( ) 0.02; ( ) 0.04; ( ) 0.06; ( ) 0.08.
made about initial conditions, and the initial (spatial) size of the infected population. Thus we can say: K An increase in the initial width of the inoculum L increases the predicted ultimate success of the invasion of infection. K An increase in the loss rate of maternal transmission m decreases the predicted ultimate success of the invasion of infection. K An increase in the loss rate of maternal transmission m reduces the predicated eventual equilibrium level of infection after successful invasion.
Additionally, unlike the simulation results when semi-infinite initial condition (I) was used, significant variation is evident between the results of the three models when the finite initial condition (II) was used. This indicates that
assumptions about the dispersal behaviour (probability distribution functions) of the insects and the precise modelling technique chosen also affect the prediction success of invasion. The reaction–diffusion model requires a larger initial inoculum width L than the Gaussian kernel model to achieve invasion, and in turn, the Gaussian kernel model requires a larger initial inoculum width L than the leptokurtic kernel model. Discussion This paper has examined the variation in results obtained from simulations of different mathematical models of infection spread that use varying explicit assumptions about insect dispersal patterns. Some comparison is made between these results and earlier results and field
130
P. SCHOFIELD
data. The model of Turelli & Hoffmann (1991) is extended with the inclusion of a maternal loss rate which raises the threshold for invasion. This is shown to slow down and eventually reverse the invasion wave with increasing rate of loss. This would confirm further the unsatisfactory nature of the diffusion model. Some analytical results for various classes of dispersal kernel with birth death kinetics the include a simple Allee effect have been investigated (Kot et al., 1996). They have also shown that modelling the movement of insects with the diffusion model or the equivalent Gaussian kernel, integro-difference equation model may not be appropriate, leading to significant underestimation of wave speeds. The results presented in this paper on the whole agree with this conclusion and with that of Turelli and Hoffmann, that the actual spread of Wolbachia in the California D. simulans population is not best modelled with a diffusion-based model. More generally, the results show that the transient dynamics of invasion subject to a threshold similar to an Allee effect are sensitive to the modelling scheme chosen, and it is possible that asymptotic solutions for the reaction–diffusion and Gaussian dispersal models may not provide robust results. The results also tend to confirm the views of other research (Minogue, 1989; Shaw, 1994) that for some populations a better model of redistribution is achieved using ‘‘fat-tailed’’ non-exponentially bounded kernels. It is possible that this might be the case for the California D. simulans data (Turelli & Hoffmann, 1991)Fthe leptokurtic kernel giving faster wave speeds for the same value of s than the Gaussian. With the model chosen here, this was enough to account for the discrepancy between the observed movement of Drosophila and the estimated dispersal rate s (Turelli & Hoffmann, 1991). However, it is possible that this faster wave speed could be produced by an additional convection term in the reaction–diffusion model that reflected some external pressure, such as human activity transporting fruit or abiotic factors such as prevailing winds. This study is also in keeping with other studies of insect dispersal systems (Hasting, 2000) which suggest that non-Gaussian dispersal kernels may better fit the empirical data.
In showing that the wave speed is sensitive to the dispersal kernel chosen in the integrodifference model, this would suggest that arbitrary choices of dispersal kernel based purely on ease of mathematical analysis may miss important aspects of the system being modelled. This may be particularly important when it comes to the success of invasion from a finite inoculum into a population where rapid dispersal may be disadvantageous to the ultimate success of the invasion Kot et al. (1996) reanalysis of Drosophila pseudoobscura dispersal data from Dobzhansky & Wright (1947) showed that the data were not sufficient to distinguish between an exponentially bounded and a non-exponentially bounded dispersal function. The conclusion was made that better data concerning the probability of rare extremely long distance dispersal are needed to distinguish between the two functions. Attempts to discover more accurate dispersal probability distributions and parameters from field data have been made. A series of papers on mark recapture experiments with Drosophila sp. (Jones et al., 1981; Coyne et al., 1982, 1987) indicated that individuals released in desert conditions can be found to have travelled much further from the release site in a short period of time than diffusion estimates would predict. These were studies of organisms dispersing in non-native environments, and, as pointed out in the Introduction with the possible problems of using Drosophila pseudoobscura parameter estimates in a model of Drosophila simulans, estimates derived for species in non-native environment may not be adequate to model the species in its natural habitat (Lewis, 1997). I am very grateful to Michael Turelli for his guidance and suggestions during this work and to the University of California Davis for hosting me. I am also grateful for the advice of Professor Chaplain and the comments of two anonymous referees in helping revise the manuscript. This work has been funded in part by the BBSRC and by a special travel grant from the Scottish International Education Trust. REFERENCES Barton, N. H. (1979). The dynamics of hybrid zones. Heredity 43, 341–359. Breeuwer, J. & Werren, J. H. (1990). Microorganisms associated with chromosome destruction and
SPATIALLY EXPLICIT MODELS
reproductive isolation between two insect species. Nature 346, 558–560. Coyne, J. A., Boussy, I. A., Prout, T., Bryant, S. H., Jones, J. & Moore, J. A. (1982). Long distance migration in Drosophila. Am. Nat. 119, 589–595. Coyne, J. A., Bryant, S. H. & Turelli, M. (1987). Long distance migration of Drosophila 2. Presence in desolate sites and dispersal near desert oasis. Am. Nat. 129, 847–861. Dobzhansky, T. & Wright, S. (1947). Genetics in natural population. XV. Rate of difussion of a mutant gene through a population of Drosophila pseudoobscura. Genetics 32, 303–324. Dunbar, S. R. (1983). Travelling wave solution of diffusive Lotka–Volterra equations. J. Math. Biol. 17, 11–32. Fisher, R. (1937). The wave of advance of an advantageous gene. Ann. Eugenics 7, 353–369. Ghelelovitch, S. (1952). Sur le determinime genetique de sterilite dans le croisement entre differentes souch de Culex autogenicus Roubaud. C. R. Acad. Sci. Paris 24, 2386–2388. Hasting, A. (2000). Parasitoid spread: lessons for and from invasion biology. In: Parasitoid Population Biology (Hochberg, M. E. & Ives, A. R., eds), pp. 70–82. Princeton, NJ: Princeton Unviersity Press. Hertig, M. & Wolbach, S. (1924). Studies in rickettsialike microorganism in insects. J. Med. Res. 44, 329–374. Hoffmann, A. A., Turelli, M. & Simmons, G. (1986). Unidirectional incompatibility between populations of Drosophila simulans. Evolution 40, 692–701. Jones, J., Bryant, S. H., Lewontin, R., Moore, J. A. & Prout, T. (1981). Gene flow and the geographicaldistribution of a molecular polymorphism in Drosophilapseudoobscura. Genetics 98, 157–178. Kimura, M. (1953). ‘Stepping stone’ model of population. Annual Report of the National Institute of Genetics, Japan, Vol. 3, 62–63. Kot, M. (1992). Discrete-time travelling waves: ecological examples. J. Math. Biol. 30, 413–436. Kot, M. & Schaffer, W. M. (1986). Discrete-time growth-dispersal models. Math. Biosci. 80, 109–136. Kot, M., Lewis, M. A. & van den Driessche, P. (1996). Dispersal data and the spread of invading organisms. Ecology 77, 2027–2042. Lewis, M. A. (1997). Spread if an invading population. In: Spatial Ecology (Tilman, D. & Karevia, P. M., eds),
131
pp. 46–69. Monographs in Population Biology, No. 30. Princeton, NJ: Princeton University Press. Lewis, M. A. & Kareiva, P. M. (1993). Allee dynamics and the spread of invading organisms. Theoret. Popul. Biol. 43, 141–158. Minogue, K. (1989). Diffusion and spatial probability models for disease. In: Spatial Components of Plat Disease Epidemics (Jeger, M., ed), pp. 127–143. Englewood Cliffs, NJ: Prentice-Hall. Mitchell, A. & Griffiths, D. (1980). The Finite Difference Method in Partial Differential Equations. New York: John Wiley & Sons. Mollison, D. (1977). Spatial contact models for ecological and epidemic spread. J. R. Stat. Soc. B 39, 283–326. Murray, J. (1989). Mathematical Biology. Biomathematics Texts, No. 19. Berlin: Springer-Verlag. Nagylaki, T. (1975). Conditions for the existence of clines. Genetics 80, 595–615. Okubo, A. (1980). Diffusion and Ecological Problems: Mathematical Models. Berlin: Springer-Verlag. O’neill, S. L., Hoffmann, A. A. & Werren, J. H. (eds) (1997). Influential Passangers Inherited Microorganisms and Arthropod Reproduction. Oxford: Oxford University Press. Renshaw, E. (1991). Modelling Biological Populations in Space and Time, Cambridge Studies in Mathematical Biology. Cambridge: Cambridge University Press. Shaw, M. (1994). Modeling stochastic processes in plant pathology. Annu. Rev. Phytopathol. 32, 523 –544. Turchin, P. (1998). Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants. Sinauer Associates. Turelli, M. & Hoffmann, A. A. (1991). Rapid spread of an inherited incompatability factor on California Drosophila. Nature 353, 440 – 442. Turelli, M. & Hoffmann, A. A. (1995). Cytoplasmic incompatibility in Drosophila simulans: dynamics and parameter estimates from natural populations. Genetics 140, 1319–1338. Turelli, M., Hoffmann, A. A. & Mckechnie, S. W. (1992). Dynamics of cytoplasmic incompatability and mtDNA variation in natural Drosophila simulans populations. Genetics 132, 713–723. Werren, J. H. (1997). Biology of Wolbachia. Annu. Rev. Entomol. 42, 587–609. Zwillinger, D. & Kokoska, S. (2000). CRC Standard Probability and Statistics Tables and Formulae. London/ Boca Raton, FL: Chapman & Hall/CRC.