Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Q1 16 17 Q5 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Mathematical modelling of spatial sorting and evolution in a host–parasite system Matthew H. Chan a, Richard Shine b, Gregory P. Brown b, Peter S. Kim a a b

School of Mathematics and Statistics, University of Sydney, NSW 2006, Australia School of Biological Sciences, University of Sydney, NSW 2006, Australia

H I G H L I G H T S

    

We model the evolution and spatial distribution of traits in cane toad and lungworm populations. We incorporate life-history trade-offs affected by these traits and examine the resulting spatial distribution of traits. The model predicts prepatent period and larval size of parasites evolve according to host densities, which matches previous empirical studies. The host wavespeed generated matches data on the invasion history of cane toads in Australia. Reproductive or survival disadvantage is likely to be small for strong spatial sorting to occur.

art ic l e i nf o

a b s t r a c t

Article history: Received 29 October 2014 Received in revised form 13 April 2015 Accepted 3 June 2015

There have been numerous empirical and agent-based modelling studies on the spatial self-structuring of traits, particularly in regard to dispersal ability (termed spatial sorting) of cane toads in northern Australia, but few mathematical modelling studies. In this study, we formulate a reaction–diffusion based partial-integro-differential equation model based on an earlier model by Bouin et al. (2012) to examine this spatial self-structuring of traits in both a cane toad population and lungworm parasite population, which evolves with the cane toad population. In particular, the traits we focus on are dispersal ability for the cane toad population and both prepatent period and larval size for the lungworm parasite population. Apart from the spatial self-structuring of these traits, our results confirm a number of observations made in empirical and agent-based studies; particularly, that there is a noticeable lag between the host and parasite population which is critically dependent on the parasite functional response to host densities, that older populations regress back to lower dispersal speeds and that spatial sorting can still occur with a disadvantage in reproductivity and/or survival in more motile individuals. Moreover, we find that such a disadvantage in reproductivity and/or survival is unlikely to be large if spatial sorting is to have a noticeable effect on the rate of range expansion, as it has been observed to Q2 have over the last 60 years in northern Australia. & 2015 Published by Elsevier Ltd.

Keywords: Population dynamics Evolutionary dynamics Spatial structure

1. Introduction A phenomenon that is often observed in biological invasions is an accelerating rate of range expansion. The current literature suggests that in many cases this is due to increased dispersal ability at the leading front of the invasion caused by “spatial sorting”, whereby heritable traits that enhance rates of dispersal accumulate at the leading front. This leads to increased rates of dispersal in successive generations due to interbreeding and produces offspring at the front with higher mean dispersal rate (Shine et al., 2011). Spatial sorting has

E-mail addresses: [email protected] (M.H. Chan), [email protected] (R. Shine), [email protected] (G.P. Brown), [email protected] (P.S. Kim).

been observed in many populations, such as cane toads (Phillips et al., 2006), butterflies (Hughes et al., 2007) and bush crickets (Simmons and Thomas, 2004). Shine et al. (2011) further asserted that such a phenomenon allows individuals with higher rates of dispersal to invade even if they do not necessarily exhibit a greater lifetime reproductive success than their competitors, in contrast to the theory of classical natural selection. Although this has not been theoretically confirmed, both mathematical and agent-based models have provided results that support this hypothesis (Shine et al., 2011; Phillips et al., 2008; Bouin et al., 2012). While the majority of studies have concentrated on examining the spatial sorting of traits of free-living organisms, there is likely to be an influence on selective pressures on life-history traits in parasites infecting hosts undergoing spatial sorting. This is due to rapidly evolving parasite life-history strategies according to host densities; for example, at low host densities, there

http://dx.doi.org/10.1016/j.jtbi.2015.06.027 0022-5193/& 2015 Published by Elsevier Ltd.

Please cite this article as: Chan, M.H., et al., Mathematical modelling of spatial sorting and evolution in a host–parasite system. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.06.027i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85

2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

M.H. Chan et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

should be a strong selection pressure for traits which maximize the probability of transmission, although trade-offs with other parasite traits imply that the optimal life-history strategy varies with host densities. The cane toad (Rhinella marina) in northern Australia and lungworm nematode Rhabdias pseudosphaerocephala, a common parasite of cane toads, are a well-studied example of such a host– parasite system. Cane toads are a highly successful invasive species in the various areas it has been introduced (e.g. Fiji, New Guinea, Philippines and various island countries in the Caribbean), usually for the purpose of controlling damage to crops by pests. Their invasion in northern Australia has been no exception; since its introduction in 1935, as an attempt to control damage in sugarcane crops in Queensland, the cane toad invasion front has accelerated with time due to spatial sorting. Current estimates place the speed of the moving front at 55–60 km/year with the front having passed the border into Western Australia in 2009 (Shine, 2010; Brown et al., 2014). Rhabdias pseudosphaerocephala is the most common lung parasite in Australian cane toads, with over 80% prevalence in a sample of 580 toads in a study by Barton (1998). The toads at the invading front lack the parasite due to low transmission rates at low host densities; this has resulted in the parasite invading front lagging behind its host's by approximately 2 years. The parasite exhibits a direct life-cycle; infective larvae invade a cane toad through the skin or alimentary tract and enter the lungs, where they feed on the host and mature. Eggs released in the lungs are passed onto the digestive tract where they hatch and enter the soil with the host faeces which completes the lifecycle. Rhabdias pseudosphaerocephala also exhibits spatial sorting in the form of life-history traits rather than dispersal ability. This is due to the rapid evolution of life-history traits associated with transmission and survival in response to changing host densities. In particular, low host densities at the invading front should favour the selection of variations in life-history strategies which enhance the parasite's ability to successfully infect hosts and reproduce. This was evidenced by Kelehear et al. (2012) who compared life-history traits among different lungworm populations, suggesting that those at the range edge exhibited larger eggs, larger free-living adults and larger infective larvae, due to increased parasite juvenile mortality at the range edge, and reduced prepatent period (reduced age at maturity). This is presumably a form of life-history trade-off; large larvae have higher establishment success but maternal fitness at high host densities is maximised by producing many small offspring even though their establishment success is lower. Similarly, parasites with shorter prepatent period exhibit reduced longevity and lifetime fecundity. Mathematical modelling of accelerating dispersal, such as in the case of the cane toads, has received increasing attention over the recent years. It is well known from mark-recapture studies that many populations disperse according to a leptokurtic (fat-tailed) dispersal kernel as opposed to a Gaussian dispersal kernel, as is commonly assumed in modelling studies. Such fat-tailed dispersal kernels have been attributed to population heterogeneity and can generate an increased probability of long-distance dispersal, which can in turn generate accelerating dispersal (Stover et al., 2014; Hastings, 2005). They have been the main ingredient of various models, such as integro-difference equations and fractional reaction-diffusion models (Kot et al., 1996; Baeumer et al., 2008). However, these models intrinsically include a fat-tailed dispersal kernel, without modelling the population heterogeneities which generate the kernel. Moreover, they assume a fixed dispersal kernel, whereas in reality it is likely to be changing over time. It is only until recently that there has been focus on spatial sorting as a source of the leptokurtosis (see Phillips et al., 2008; Perkins et al., 2013; Stover et al., 2014). In this study, we examine the spatial variation of fitness traits of both a host and parasite, where we use the cane toad and lungworm

R. pseudosphaerocephala as a specific example of such a host–parasite system. We formulate a partial integro-differential equation (PIDE) model, accounting for traits influencing dispersal ability for the cane toads and prepatent period and larvae size for the lungworms. Firstly, we qualitatively and quantitatively examine the acceleration of the cane toad invading front due to spatial sorting, and how this varies with trade-offs in survival and reproduction associated with increased dispersal ability. Secondly, we qualitatively examine the variation in life-history traits of the lungworm population emerging from the host undergoing spatial sorting and determine how this variation is influenced by host densities. Although there have been previous agent-based modelling studies done on this problem, to our knowledge this is the first PIDE model in the literature which describes the evolution of life-history traits in this specific host–parasite system. Note that we do not examine whether the lungworm nematode Rhabdias pseudosphaerocephala is a suitable biological control for the cane toad invasion; instead, we focus on the interplay between the evolutionary and spatial dynamics of the host–parasite system, and compare with observations made in previous empirical studies.

2. Model formulation Denote the host and parasite population density at position x and time t by uðx; θ; tÞ and vðx; a; γ ; tÞ, where θ A ½0; 1 is a measure of dispersal ability, a A ½0; 1 is a measure of parasite larval size and γ A ½0; 1 is a measure of the parasite prepatent period. We assume that parasite larval size is positively correlated with survivability in the environmental stage of the parasite lifecycle and negatively correlated with number of offspring produced. We model spatial movement using a reaction–diffusion system where the spatial diffusion coefficient is dependent on the dispersal ability of the cane toad. The mutation of the life-history traits θ, a and γ are modelled with diffusion terms, where the diffusion coefficient for these diffusion terms is assumed to be constant. We assume that if there is an absence of hosts, then there will be an absence of R1 parasites at a particular x position, that is, if 0 u dθ ¼ 0 then R1 R1 0 0 v da dγ ¼ 0 for a particular x. This is equivalent to assuming that the dispersal rate of the parasites is directly tied with the dispersal rate of the cane toads. We denote the spatial diffusion coefficient for the cane toads by f ðθÞ (defined in Section 3) and define the spatial diffusion coefficient for the parasite population R1 f ðθ Þu dθ by the normalised average of f ðθÞ, that is, 0R 1 . Note that we u dθ

0 R1 f ðθÞu dθ 0 have set limu-0 R 1 ¼ 0 as this is the biologically relevant 0

u dθ

choice for the limit. For the cane toad population, we assume logistic growth with the assumption of homogeneous mixing. We denote the reproductive rate rh and mortality rate μh and let them be functions of θ. Throughout this study, we only consider r h ðθÞ and μh ðθÞ to be nonincreasing and non-decreasing respectively with a linear relationship in θ, representing a cost in enhanced dispersal. They are given by r h ðθÞ ¼ r h ð1  cθÞ;

ð2:1Þ

μh ðθÞ ¼ μ h ð1 þ dθÞ:

ð2:2Þ

Similar to the equation for cane toads, we denote the reproductive rate and mortality rate of the parasites by rp and μp respectively. These are defined such that the following trade-offs, mentioned in Section 1, are captured in the dynamics: 1. Low host density selects for low vice versa.

γ (prepatent period) and

Please cite this article as: Chan, M.H., et al., Mathematical modelling of spatial sorting and evolution in a host–parasite system. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.06.027i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M.H. Chan et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

uðx; θ; 0Þ ¼ δðxÞδðθÞ and vðx; a; γ ; 0Þ ¼ δðxÞδða  0:5Þδðγ  0:5Þ. The system was solved using a finite-difference scheme, see the Appendix for details.

2. Low host density selects for high a (larval size) and vice versa. Because virus and parasite transmission are density-dependent and also because of trade-offs (1) and (2), we let the reproductive rate depend on the availability of hosts at a particular position x, γ and a. We let the reproductive rate decrease linearly with a and have a Holling-type   functional response in γ and define it to be of ^ γ ð1 kaÞ, where the form r p 〈u〉;     ^ γ ¼ r p 〈u〉 ^ b 2ð1  γ Þð1 〈u〉 ^ q Þ þ γ 〈u〉 ^ q ; ð2:3Þ r p 〈u〉; R1 u 〈f 〉 ¼ 0 f ðθÞ dθ, u^ ¼ rh ð0Þ  μ ð0Þ and r p , b, m, q, k are constants. h ^ is a rescaling of the host density so that it falls in the Essentially, 〈u〉 interval ½0; 1 for convenience, as the equilibrium host density is not equal to 1. Thus, the infection rate of toads is then given by vr p and we let the growth be controlled by density dependence in v. Parameter b controls the functional response in host availability, k controls the proportion of reproductive loss due to a higher a value and r p denotes the base reproductive rate for individuals with γ ¼ 1 and a¼0 when ^ ¼ 1. Parameters q and m control the intensity of trade-off (1); m 〈u〉 determines how many times more successful per-capita an individual with γ ¼ 0 is, compared to an individual with γ ¼ 1 at low host densities and q controls the shape of the functional response. This form for rp ensures selection for a long prepatent period at high host densities, which declines as host density declines. See Fig. 1a for a plot     ^ γ ; for example, when 〈u〉 ^ γ is ^ is less than 0.5, r p 〈u〉; of r p 〈u〉; decreasing with respect to γ. The function for parasite mortality is required to be a decreasing function of the availability of hosts as per trade-off (2); we define it to be of the form   ð1  aÞ ^ aÞ ¼ μ p 1 þ μp ð〈u〉; ; ð2:4Þ ^ 7 0:05 þ 5〈u〉

3. Parameter estimates Throughout the study, we use estimated parameter values for the cane toad population and hypothetical parameter values for the parasite population. The reproductive rate r h of toads was calculated by assuming an average clutch size of 15,000 eggs, a maturation period of 18 months and that 0.5% of these eggs hatch and survive to sexual maturity. We assume that the significant proportion of mortality within a given period of time is due to predation and estimate the base mortality rate μ h from data on radio-tracked toads using a least-squares fit of the model y ¼ e  μt (i.e. assuming a proportional mortality rate), see Fig. 2. The cost in reproductive and mortality rate associated with enhanced dispersal ability, c and d in Eq. (2.5a), was set to an arbitrary value of 0 and 0.1, although we vary both of these parameters in Section 5 (see Fig. 9). The diffusion rate f ðθÞ was calculated based on the minimum wave speedpformula for Fisher's equation with a linear ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi mortality term v ¼ 2 f ðθÞðr h  μ h Þ while also noting that the speed of the invasion front of the cane toad population in 1950 was approximately 10 km/year. Without loss of generality, we letffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi this wavespeed correspond to θ ¼ 0. Note that v ¼ 2 f ðθÞðr h  μ h Þ is merely an approximation(an underestimate actually, see Bouin et al., 2012 for a proof) of the wavespeed, although we expect this approximation to be accurate as Dθ is expected to be small. 2 Assuming a quadratic form v ¼ v0 þ v1 θ for the wavespeed, based on data on the invasion history from (Phillips et al., 2006), yields

where μ p , f, g and h are constants. We have let μp scale linearly with a for simplicity, similar to rp. Parameter μ p is the base mortality rate associated with a¼ 1, f is a measure of the additional morality suffered with decreasing a; thus μ p =f is the additional mortality suffered for ^ ¼ 0. Parameters g and h control the individuals with a¼0 when 〈u〉 shape of the functional response to host density. See Fig. 1b for a plot ^ aÞ. A concise list of parameter definitions and values is given of μp ð〈u〉; in Tables 1 and 2. The model is then given by 8 ∂u ∂2 u ∂2 u > > ¼ f ðθÞ 2 þ Dθ 2 þ r h ðθÞuð1  〈u〉Þ  μh ðθ Þu; > > > ∂t ∂x > ∂θ > < ∂v ∂2 v ∂2 v 〈fu〉 ∂2 v ¼ Da 2 þ Dγ 2 þ > > ∂t 〈u〉 ∂x2 ∂a ∂γ > >  >   R1 R1 > > ^ ^ aÞv: : þ vr p 〈u〉; γ ð1  kaÞ 1  0 0 v da dγ  μp ð〈u〉;

3

f ðθ Þ ¼

 1 2 2 v0 þv1 θ ; 4ðr h  μ h Þ

ð3:1Þ

where v0 and v1 represent the wavespeeds corresponding to θ ¼ 0 and θ ¼ 1 respectively. Throughout this study we set v0 ¼ 10 km/ year and v1 ¼ 70 km/year.

4. Stability analysis In this section, we examine the equilibrium solution of the model, assuming a finite spatial domain, and show that at equilibrium, slow diffusers dominate in the host population, while the parasite population responds by reverting to a long prepatent period (γ ¼ 1) and small larval size (a¼ 0) to maximize reproductivity. We first consider the equation for the cane toad population, Eq. (2.5a). For rh and μh non-constant, we discretise the system with respect to θ into n compartments, indexed by subscript i, to avoid working with a non-autonomous system. Thus, we obtain

ðaÞ

ð2:5Þ ðbÞ

For simplicity, we set no-flux boundary conditions, i.e. all partial derivatives of u and v are equal to 0, and the initial condition to be

0.1

20 0.05

10 0 0

0 1

1 0.5 0 0

0.5 γ

1

0.5

0.5

a 1 0

Fig. 1. Plots of rp and μp using parameter values from Table 2.

Please cite this article as: Chan, M.H., et al., Mathematical modelling of spatial sorting and evolution in a host–parasite system. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.06.027i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M.H. Chan et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

solution tends to Uð0Þ when the solution reaches the boundary and spatial sorting ceases. Since Jð0Þ clearly has full rank and is thus invertible, by the Implicit Function Theorem there exists a smooth local mapping gðεÞ such that ðgðεÞ; εÞ is an equilibrium of the system, where ε is contained within an open neighbourhood around 0 and gðεÞ is contained within a local neighbourhood of Uð0Þ. This shows that Uð0Þ perturbs smoothly under ε, and so we expect that UðεÞ to retain the hyperbolic and locally asymptotically stable attributes of Uð0Þ. In a similar way, it is possible to deduce the stability for the parasite equation. Due to system (2.5) being decoupled, we assume that u is at equilibrium. We discretise a and γ, indexed by i and j respectively, to obtain

1

0.8 Proportion surviving

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 Q4 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

0.6

0.4

0.2

0

0

50

100 Days

0 1 n n n X n X X X ∂vi;j ∂2 vi;j M 1i;k vk;j þ ε2 M 2j;k vi;k þ C 1 2 þ vi;j r j ð1  kai Þ@1  vi;j A  μi vi;j ; ¼ ε1 ∂t ∂x j¼1i¼1 k¼1 k¼1

150

Fig. 2. The proportion of toads that were radiotracked and killed by predators (blue markers) and the proportion of toads that were radiotracked but not killed by predators (red markers). The black line represents a least-squares fit using the model y ¼ e  μt with μ ¼ μ h . (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

the system

ð4:6Þ R1 f ðθÞu dθ and M 1i;j and M 2i;j denote the mutation where C 1 ¼ 0R 1 0

matrices for a and are given by ∂v0i;j

0

n n X X ∂ui ∂2 u ¼ Di 2i þ r i ui @1  uj A  μi ui þ ε M i;j uj ; ∂t ∂x j¼1 j¼1

¼ε

∂vi;j

1

i ¼ 1; 2; …; n; ð4:1Þ

where M i;j is the matrix of mutation rates, with M i;i 4 M i;j for j a i. This is similar to a system examined in a study by Dockery et al. (1998), who showed that the equilibrium U 1 ð0Þ corresponding to ε ¼ 0 perturbs smoothly under ε and thus U 1 ðεÞ retains the stability properties of U 1 ð0Þ, for ε sufficiently small. We follow this approach for System (4.1). We first consider the case where ε ¼ 0. The Jacobian of the system is given by 8 ! n X > > < εM i;i þ r i 1  ui  uk  μi if i ¼ j; ð4:2Þ J i;j ðεÞ ¼ k¼1 > > : εM i;j  r i ui if i aj: Let 1 rk rn and consider steady states of the form 8 < uk ¼ 1  μk if i ¼ k rk : u ¼ 0 if i a k: i Evaluating the Jacobian at such a steady state gives 8 εMi;j if i aj a k; > > > < μk if i ¼ ja k; J i;j ðεÞ ¼ εM i;i þ r i r  μi k > > > : εMi;j þ μi  r i if i ¼ k:

u dθ

γ respectively. The derivatives for the Jacobian 0

1 2 @ 1 M i;i vi;j þ 2 M j;j vi;j þ hj ð1  kai Þ 1 

ε

n X n X

1 vi;j A  vi;j r j ð1  kai Þ  μi ;

j¼1i¼1

∂v0i;j ∂vi;l ∂v0i;j ∂vl;j ∂v0i;j ∂vl;k

¼ ε2 M 2j;l  vi;j r l ð1  kai Þ;

j a l;

¼ ε1 M 1i;l  vi;j r j ð1  kal Þ;

i a l;

¼  vl;k r k ð1  kal Þ;

where v0i;j ¼

∂vi;j ∂t .

ð4:7Þ

i a l and ja k;

We define V ðε1 ; ε2 Þ to be the equilibrium of the

parasite equation. Similarly to above, consider ε1 ¼ ε2 ¼ 0 and the equilibrium solution Vð0; 0Þ of the form 8 μk < vk;l ¼ 1  hl ð1  kak Þ ð4:8Þ : v ¼ 0 if ia k and j al; i;j where we assume that

ð4:3Þ

ð4:4Þ

Letting ε ¼ 0, the eigenvalues are given by the diagonal elements, that is 8 if i ¼ k; < μk r k λi ¼ r μk  μ ð4:5Þ if ia k: : i i rk We assume that r i 4 μi to avoid the trivial steady state ui ¼0 for all i being locally asymptotically stable. Assuming that ri is monotonically decreasing and μi is monotonically increasing, such that u1 and un corresponds to the slowest and fastest disperser respectively, the only stable equilibrium of the form given by Eq. (4.3) is when k ¼1 as r i μr11  μi o 0 for i Z2. We denote this equilibrium as Uð0Þ. Thus, when there is a reproductive disadvantage and increased mortality for increasing dispersal ability θ, the

maxhj ð1  kai Þ=μi ¼

hl

μk

ð1 kak Þ;

ð4:9Þ

where hl ð1  kak Þ 4 μk . The Jacobian evaluated at this equilibrium has eigenvalues given by its diagonal. That is, ( hj ð1  kai Þð1  2vk;l Þ  μi if i ¼ k and j ¼ l λ ¼ h ð1  ka Þð1  v Þ  μ ð4:10Þ if ia k or j al: j i k;l i  h ð1  ka Þμ Since hj ð1  kai Þð1  vk;l Þ  μi ¼  μi 1  hj ð1  kai Þμk , by assumpl

k

i

tion (4.9) all eigenvalues are negative and thus the equilibrium Vð0; 0Þ is locally asymptotically stable. The Jacobian clearly has full rank and hence by the above reasoning for the host equation Eq. (2.5a), we expect Vðε1 ; ε2 Þ to retain the hyperbolic and locally asymptotically stable attributes of Vð0; 0Þ for sufficiently small and

ε2.

ε1

5. Results Plots of the total population densities varying over time and space are shown in Fig. 3; a more in-depth view of the solution at certain t values are given in Figs. 4–6. In Fig. 4 we show plots of u

Please cite this article as: Chan, M.H., et al., Mathematical modelling of spatial sorting and evolution in a host–parasite system. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.06.027i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M.H. Chan et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

4

x 10

x 10 2.5

2.5

0.8

0.8 0.7

0.7 2

2

0.6 Time t (days)

Time t (days)

0.6 0.5

1.5

0.4 1

0.3 0.2

0.5

0.5

1.5

0.4 1

0.3 0.2

0.5

0.1

0.1 0

0

5

10

15

20

0

0

0

5

10

15

20

0

Distance x (100km)

Distance x (100km) Fig. 3. Plot a and b represent

R1 0

u dθ and

R1 R1 0

0

v da dγ respectively, varying over t and x.

1

0.8

2

0.8

0.6

Dispersal ability θ

Population density

0.7

0.5 0.4 0.3 0.2

1.5 0.6 1 0.4 0.5

0.2

0.1 0

0

5

10

15

0

20

0

5

Distance x (100km)

10

15

20

0

Distance x (100km)

1

1

16

2

14 0.8

1.5 0.6 1

0.4

0.5

0.2

Prepatent period γ

0.8

Size a

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

5

12 10

0.6

8 0.4

6 4

0.2

2 0

0

5

10

15

20

0

0

5

Distance x (100km)

10

15

20

0

Distance x (100km)

R1 R1 Fig. 4. Plots of uðx; θ; tÞ and vðx; a; γ; tÞ for t ¼ 15; 000 days. The solid and dashed lines in (a) represent 0 u dθ and 0 0 v da dγ respectively. The colours in plots (b), (c) and R1 R1 (d) are an indicator of the magnitude of u, 0 v dγ and 0 v da, where blue and red represent a low and high value respectively. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.) R1

and v with respect to θ, γ and a at t¼ 15,000, with the initial host population at x¼ 0 and θ ¼ 0 and initial parasite population at a ¼ γ ¼ 1 and x ¼0. In Fig. 4a we show the cane toad and parasite population densities along x, which exhibits a clear lag between the host and parasite population resulting from a declining host density over x. This declining host density is due to the front being occupied by individuals with fast dispersal ability (i.e. spatial

sorting), for which there is an extra mortality cost, as observable in Fig. 4b, where dispersal ability θ is plotted against x and where the colours represent population density levels (blue and red correspond to low and high receptively, according to a logarithmic colour mapping). Note that the cost in enhanced dispersal is the reason for a high population density along θ ¼ 0. Fig. 4c and d shows the parasite population corresponding to the host

Please cite this article as: Chan, M.H., et al., Mathematical modelling of spatial sorting and evolution in a host–parasite system. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.06.027i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M.H. Chan et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

6

population (Fig. 4c). From Fig. 4c, there is an evident shift to lower larval sizes a in the parasite population near x ¼0, since the optimal strategy for the parasites is to increase offspring numbers when host densities are high. Conversely, Fig. 4c shows strong selection for larger larval sizes at the front  due to declining host ^ γ in Eq. (2.3), the densities. However, due to the form of r p 〈u〉; selection pressure is not as strong in the prepatent period γ towards the parasite wave front, as shown by the relatively low population density deviating from γ ¼ 1 in Fig. 4d. Fig. 5 shows the solution of u and v at t ¼20,000. As observable in Fig. 5a, the lag between the host and parasite population is still present, but is greatly reduced. This is due to the parasites having adequate time to adapt to the maximum values of a and γ, shown in Fig. 5c and d. The curve between approximately x¼5 to x ¼10 in the parasite population density is due to a combination of decreasing host densities and adaptation to the extremes of a and γ by the parasites along these x values. Moreover, the range of the host population with dispersal ability θ ¼ 0 is increasing over time, evident from comparison of Figs. 4b and 5b. This follows the results in Section 4 and matches observations from empirical studies that suggest leg lengths of individuals decline to shorter lengths over time (Phillips et al., 2006). These qualitative behaviours observed in Fig. 5 are maintained for larger t, as shown in Fig. 6. Note that in Section 4 we have assumed a finite domain for x, which is precisely why it is possible to obtain an equilibrium solution along x. This would not hold in general for an infinite domain as there will be spatial sorting, however we expect that the equilibrium solution is achieved behind the wavefront in older

populations, as there is little interaction between new populations at the front and older populations behind the front. This indicates that a cost of dispersal must be included in a model to capture this decline over time in the mean θ value of the host population. As expected from the density-dependent nature of parasite transmission, and indeed evident from Figs. 4–6, the lag is sensitive to the parameters q and b controlling the functional response of parasite growth to host density; this can be observed in Figs. 7 and 8. As observable in Fig. 8, it is possible for the lag to increase by large amounts during the process of spatial sorting in the host due to a combination of the parasite adapting to lower host densities and lower overall reproductive rates at such host densities, which is heavily dependent on the value for q. However, note that the lag distance always reaches an equilibrium when the mean θ of the host population stabilizes at its maximum value 1. From Fig. 7, the wavespeed produced by the model with parameters values given by Table 1 matches data of the invasion history of cane toads in northern Australia for earlier years (t¼0 to t ¼ 32 years), but seemingly underestimates the invasion speed for later years (at t ¼ 55 years) based on the available data (Phillips et al., 2006). In Fig. 9 we show plots of dispersal ability θ over x, similar to Figs. 4b, 5b, and 6b, given a certain measure of reproductive disadvantage (the magnitudes of c and d), where each simulation was stopped when the wavefront reached x¼20. Corresponding plots showing the position of the wavefront over time is shown in Fig. 10. This provides a simple sensitivity analysis to determine the degree of spatial selection in the host population. Clearly, the

1

0.8

2

0.7 0.6

Dispersal ability θ

Population density

0.8

0.5 0.4 0.3 0.2

1.5 0.6 1 0.4 0.5

0.2

0.1 0

0

5

10

15

0

20

0

5

10

15

20

16

1

1

16

14 12 10

0.6

8 0.4

6 4

0.2

14

0.8 Prepatent period γ

0.8

12 0.6

10 8

0.4

6 4

0.2

2 0

0

5

10

15

Distance x (100km)

0

Distance x (100km)

Distance x (100km)

Size a

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

20

0

2 0

0

5

10

15

20

0

Distance x (100km)

R1 R1 R1 Fig. 5. Plots of uðx; θ; tÞ and vðx; a; γ; tÞ for t ¼ 20; 000 days. The solid and dashed lines in (a) represent 0 u dθ and 0 0 v da dγ respectively. The colours in plot (b), (c) and R1 R1 (d) are an indicator of the magnitude of u, 0 v dγ and 0 v da, where blue and red represent a low and high value respectively. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

Please cite this article as: Chan, M.H., et al., Mathematical modelling of spatial sorting and evolution in a host–parasite system. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.06.027i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M.H. Chan et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1

0.8

2

0.7 0.6

Dispersal ability θ

Population density

0.8

0.5 0.4 0.3 0.2

1.5 0.6 1 0.4 0.5

0.2

0.1 0

5

10

15

0

20

0

5

10

15

0

20

Distance x (100km)

Distance x (100km)

1

1

16

16

14 0.8 10

0.6

8 0.4

6 4

0.2

Prepatent period γ

Size a

14

0.8

12

12 0.6

10 8

0.4

6 4

0.2

2 0

0

5 10 15 Distance x (100km)

2

0

20

0

0

5 10 15 Distance x (100km)

0

20

R1 R1 R1 Fig. 6. Plots of uðx; θ; tÞ and vðx; a; γ; tÞ for t ¼ 25; 000 days. The solid and dashed line in (a) represent 0 u dθ and 0 0 v da dγ respectively. The colours in plot (b), (c) and R1 R1 (d) are an indicator of the magnitude of u, 0 v dγ and 0 v da, where blue and red represents a low and high value respectively. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

50

15

Velocity (km/year)

Distance x (100km)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

7

10

5

0 0

40

30

20

10

10

20

30 40 Time t (years)

50

60

10

20

30 40 Time t (years)

50

60

Fig. 7. Plots of the distance and the velocity of the host and parasite wavefront over time using parameter values as listed in Tables 1 and 2. The solid and dashed lines represent the wavefront of the host and parasite population respectively. The red circles in plot (b) mark the radial increase in range of that particular year, with t¼ 0 corresponding to 1950. Plot (b) was obtained by numerically differentiating the data points in plot (a).

degree of spatial sorting is sensitive to c and d; for no disadvantage in either reproduction or mortality, there is a high degree of spatial sorting, however this quickly disappears as c and d are increased as θ ¼ 0 becomes increasingly favourable. The maximum reproductive disadvantage before spatial sorting is severely limited occurs at approximately c ¼0.12 and d ¼0.15, with a complete lack

of spatial sorting at c Z 0:15 and d Z 0 and similarly for d Z 0:20 and c Z0. Although spatial sorting is guaranteed to be present at the wave front, strong or weak depending on the reproductive and survival costs associated with enhanced dispersal ability, a relatively small extent of spatial sorting is unlikely to have a significant impact on the invasion speed. Given the substantial

Please cite this article as: Chan, M.H., et al., Mathematical modelling of spatial sorting and evolution in a host–parasite system. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.06.027i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M.H. Chan et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

8

30

50 45 Velocity (km/year)

25 Distance x (100km)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

20 15 10 5 0

40 35 30 25 20 15 10

0

20

40

60

80

20

Time t (years)

40

60

80

Time t (years)

Fig. 8. Plots of the distance and the velocity of the host and parasite wavefront over time using parameter values as listed in Tables 1 and 2 except with b¼ 1.5 and q¼ 0.7. The solid and dashed lined represent the wavefront of the host and parasite population respectively.

Table 1 Host parameter descriptions and estimates. Symbol

Definition

Value

rh

Rate of juveniles reaching sexual maturity (including survival)

8:22  10  2 day

μh

Mortality rate

6:95  10



Mutation rate of dispersal ability

c d

Proportion of decrease in reproductive rate due to θ Proportion of increase in mortality rate due to θ

5  10  3 mm day 0 0.1

increase in speed over the cane toad invasion history, as shown in Fig. 7b, Figs. 9 and 10 show that faster dispersal speeds in cane toads are unlikely to incur significant costs in survival or reproductivity. Comparing with Fig. 7a, it is clear that the only viable (c,d) values are ðc; dÞ ¼ ð0:10; 0:00Þ, ð0:05; 0:05Þ or ð0:00; 0:10Þ.

6. Discussion In this study, we have provided a general mathematical framework to study the spatial self-structuring and evolution of traits in a host–parasite system, where we have used cane toads and the lungworm parasite Rhabdias pseudosphaerocephala as an explicit example. Numerical solutions, using rough estimates for the host and parasite population, show that the tradeoffs in reproductivity and mortality occurring with these traits are sufficient to generate spatial structuring of traits in both the host and parasite population. In addition to spatial self-structuring of traits, the model captures various other aspects described in previous empirical studies, such as spatial sorting occurring without a dependence on classical natural selection, a lag between parasite and host populations and the decline in dispersal ability over time in older cane toad populations (Kelehear et al., 2012; Phillips et al., 2006). The parasite lag occurs even though the parasites are of no detriment to the hosts; this was also found by Phillips et al. (2010) using an agent-based model. Moreover, by observing the solution under variation of both c and d, we show that a significant cost in reproductivity and survival associated with enhanced dispersal is unlikely as such a cost prevents noticeable spatial sorting, which in turn prevents a significant increase in invasion speed over time, as observed in Phillips et al. (2006). Throughout the study we have made a few simplifying assumptions to ease analysis and computation. Firstly, we have assumed that the parasite infection does not induce an increased

2

day

1 1 1

mortality rate in the cane toad population. Although the lungworm parasite does not induce a significant mortality rate in adult cane toads, it has been shown to have a significant impact on juvenile toads (Kelehear et al., 2009). We leave a detailed study of this for future work, however we note that this can be modelled by adding an extra v-dependent mortality term to the cane toad equation Eq. (2.5a), such as μ^ vu=ðε þ uÞ. Qualitatively, we expect the addition of this extra mortality term would cause the invasion front of the cane toad population to peak due to the absence of parasites, but decline as the parasites arrive, as observed in a study by Phillips et al. (2010). Secondly, the trade-off for a short prepatent period γ is reduced longevity and lifetime fecundity. We have represented the result of this trade-off (individuals with low γ are selected for at low host densities and vice versa) directly in the function rp. This is due to the lack of an age-structure component in the model; the inclusion of such a component would represent the trade-off more naturally, but would be computationally expensive and expected to achieve the same result after summing over all age compartments. A major hurdle of the study was obtaining parameter estimates for both the cane toad and the parasite, particularly for the latter. The parameter values for the parasites in Table 2 are mainly hypothetical and were chosen such that the parasite lag generated by the model matches with the lag of 1–3 years observed by Phillips et al. (2010). Due to the number of parameters, we have not included a sensitivity analysis for the parameters of the parasite as we have done for the parameters c and d for the cane toads. As expected, the parasite lag period is sensitive to the parameters, however, we have observed that the qualitative behaviours of the life-history traits a and γ across space remain the same as the parameters vary within a reasonable range. We give a simple example of this in Fig. 8. A discrepancy in the invasion speed during the later years between empirical data and our model is shown in Fig. 7b. This

Please cite this article as: Chan, M.H., et al., Mathematical modelling of spatial sorting and evolution in a host–parasite system. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.06.027i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M.H. Chan et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1

2 0.6 1.5 0.4

1

0.2

0.6 0.4

0.5 0.2 0

20

0

10

1

0.4 0.5 0.2 20

Dispersal ability θ

Dispersal ability θ

0.6

0.8 1.5 0.6 1 0.4

0

0

0.5

0.2 0

Distance x (100km)

0

1

0.4

0.5

0.2 0 20

10

20

0

20

0

0.8

2

0.6

1.5

0.4

1

0.2

0.5

0

0

10

20

0

Distance x (100km) 1 2.5

Dispersal ability θ

Dispersal ability θ

1.5 0.6

0.8

2

0.6

1.5

0.4

1

0.2

0.5

0

0

0 0

10 Distance x (100km)

0.8 2 0.6

1.5

0.4

1

0.2

0.5 0

0

20

0

10 Distance x (100km)

20

1

1

1

3

2

0.6

1.5

0.4

1

0.2

0.5 0

10

20

Dispersal ability θ

2.5

0.8

0

10

1

2

0.8

10 Distance x (100km)

0.5

0.2

Distance x (100km)

1

0

1 0.4

1

2

1.5 0.8

10

1.5 0.6

Distance x (100km)

1

0

0.8

Distance x (100km)

Distance x (100km) 1

0

2

0

20

Dispersal ability θ

10

1

Dispersal ability θ

0

0.8

0.8

2

0.6

1.5

0.4

1

0.2

0.5

0

0 0

Distance x (100km)

10 Distance x (100km)

20

Dispersal ability θ

0

0.5

1

1.5 Dispersal ability θ

2.5

0.8

Dispersal ability θ

Dispersal ability θ

1

Dispersal ability θ

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

9

0.8 2

0.6 0.4

1

0.2 0

0

10

20

0

Distance x (100km)

Fig. 9. Plots of u, where blue and red represent a low and high value respectively, for varying c and d (the cost in reproduction and survival associated with enhanced dispersal, see Eqs. (2.1) and (2.2)). The simulation was stopped when the wavefront reached 2000 km.

could be due to many reasons; we have selected a relatively simple form for the spatial diffusion coefficient f ðθÞ for the cane toad population Eq. (2.5a) based on data of the invasion history. A more accurate method of estimating of f ðθÞ is to relate θ with tibia length, as it has been observed by Phillips et al. (2006) that tibia length is positively correlated with distance moved. However, data on tibia length versus distance moved, to ascertain a clear functional relationship between the two, is still lacking due to logistical difficulties. Moreover, it is possible we have underestimated the mutation rate of dispersal ability Dθ or the rate of juveniles reaching sexual maturity r h . Lastly, it is likely that there are other traits which affect the invasion speed that are not included in our model; for example, a recent empirical study by Brown et al. (2014) showed that there is strong selection for “path straightness” at the range edge. A difficult aspect of examining spatial sorting is quantifying the extent or “strength” of spatial sorting under the trade-off between increased spatial motility and reduced reproductive and survival probability. We show in Section 5 that the system can still exhibit spatial sorting when individuals of high motility have a reproductive disadvantage, but that the degree of spatial sorting is sensitive

to the magnitude of this disadvantage. However, these results have included the advantage of reduced competition from conspecifics (due to the density dependence in Eq. (2.5a)) that motile individuals have when undergoing range expansion, and so it is unclear whether spatial sorting is occurring due solely to spatial selfstructuring of dispersal ability or a mix with classical natural selection. It is likely that a different model formulation needs to be used to examine this question, as it is difficult to quantify the benefit of reduced competition from conspecifics at range edge compared to the range core from the model in this study. Currently, we are unaware of any mathematical modelling studies which have addressed the effect of spatial sorting alone, however an agentbased model by Shine et al. (2011) demonstrated that this is true when the lifetime reproductive success is equal across θ by using a fixed reproductive rate for all individuals and global density dependence to remove correlation from classical natural selection. A crude way of showing that spatial sorting can occur with a reproductive disadvantage and without the help of reduced competition from conspecifics at the range edge is by assuming a one-dimensional spatially and temporally discrete system where at

Please cite this article as: Chan, M.H., et al., Mathematical modelling of spatial sorting and evolution in a host–parasite system. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.06.027i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M.H. Chan et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

10

15 10 5

0

10

20

30

20 Distance x (100km)

20 Distance x (100km)

Distance x (100km)

20

15 10 5

40

0

20

Time t (years)

10 5

10 5

0

20

10 5

40

0

10 5

0

20

40

60

10 5

60

10 5

80

0

50

80

100

Time t (years) 20 Distance x (100km)

Distance x (100km)

15

100

15

Time t (years) 20

40

50 Time t (years)

15

60

20

20

5

20

Time t (years)

0

10

60

Distance x (100km)

Distance x (100km)

Distance x (100km)

40

20

15

60

15

Time t (years)

20

40 Time t (years)

15

40

20

20

20

Time t (years)

0

5

0

Distance x (100km)

Distance x (100km)

Distance x (100km)

15

20

10

40

20

0

15

Time t (years)

20

Distance x (100km)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

15 10 5

0

Time t (years)

50

100

15 10 5

0

50

Time t (years)

100

Time t (years)

Fig. 10. Plots of the wavefront of u, for varying c and d. These correspond to the plots in Fig. 9. The simulation was stopped when the wavefront reached 2000 km.

Table 2 Parasite parameter descriptions and hypothetical values. Symbol

Definition

Value

rp

Rate of juveniles reaching sexual maturity (including survival)

0:1 day

μp

Mortality rate

3  10

Da

Mutation rate of lifespan in environmental stage

5  10  3 day



Mutation rate of prepatent period

b q k

Intensity of infection functional response to host density Parameter to control functional response curve Proportion of reproductive loss at high a

each time step an individual can either disperse to the left or right without reproducing, or stay at their current position and reproduce. In addition, we assume that the individuals vary in their dispersal ability θ, but we do not include mutation for simplicity. We denote this system at time t by the vector xt ðθÞ, whose entries represent the population numbers at a certain spatial position, for a particular dispersal ability θ. The solution for xt ðθÞ (see the Appendix for a

5  10 2 2 0.3

1

2

3

day day

1 1 1

derivation) is given by 0R1

1 ðð1 þ γ ðθÞÞð1  2vðθÞÞ þ 2vðθÞ cos ðπ xÞÞt sin π x sin π x dx B R0 C B 1 ðð1 þ γ ðθÞÞð1  2vðθÞÞ þ 2vðθÞ cos ðπ xÞÞt sin π x sin 2π x dx C B 0 C xt ðθÞ ¼ 2B C B⋮ C @R A 1 t ðð1 þ γ ð θ ÞÞð1  2vð θ ÞÞ þ 2vð θ Þ cos ð π xÞÞ sin π x sin n π x dx 0

Please cite this article as: Chan, M.H., et al., Mathematical modelling of spatial sorting and evolution in a host–parasite system. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.06.027i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M.H. Chan et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

15

Scaled population size

x 10 Population size

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

11

10 5 0 20 n

4 2 0 20

1

10 0 0

1

10

0.5 θ

n

0.5 θ

0 0

Fig. 11. Plot (a) shows the solution of x50 ðθÞ with γ ¼ 1, v ¼ 0:25, C 1 ¼ 0:1 and C2 ¼ 0.5. Plot (b) is equivalent to (a) except rescaled according to x50 ð0Þ.

where γ ðθÞ and vðθÞ denote the proportional increase in population numbers after one time step and the proportion of the population dispersing, both functions of dispersal ability θ. We define γ ðθÞ ¼ γ ð1  C 1 θÞ and vðθÞ ¼ vð1  C 2 θÞ. Thus they are decreasing linearly, similar to the growth term in Eq. (2.5a). Plots of xt , shown in Fig. 11, suggest that spatial sorting can still occur without help of reduced competition from conspecifics and with a small reproductive disadvantage, although the strength of spatial sorting is dampened as expected. We leave the formulation and analysis of a more sophisticated model targeting this question to a future work. Acknowledgements The work of MHC was supported by the Australian Postgraduate Award. P.S.K. was supported by the ARC Discovery Early Career Research Award. The empirical studies were funded by ARC Grants to R.S. The authors would also like to thank Mary Myerscough for very helpful suggestions on the manuscript.

n þ 1=2

where U^ is the half step approximation obtained by running n þ 1=2 (A.2) with time step Δ2t and Un replacing U^ (Douglas, 1962; Thomas, 1995).

A.2. Discrete model without density-dependent growth We consider a one-dimensional spatially and temporally discrete system where at each time step an individual can either disperse to the left or right without reproducing, or stay at their current position and reproduce. Let vðθÞ be the probability of dispersal and γ ðθÞ be the proportion of growth after one time step. Both are functions of θ, denoting dispersal ability. Then, the model of such a system is given by the matrix equation xt þ 1 ðθÞ ¼ MðθÞxt ðθÞ;

where MðθÞ is a symmetric tridiagonal toeplitz matrix defined by 0 B B M¼B B @

Appendix A

ðA:3Þ

ð1 þ γ ðθÞÞð1  2vðθÞÞ

vðθÞ



0

vðθÞ ⋮

ð1 þ γ ðθÞÞð1  2vðθÞÞ ⋮

⋯ ⋱

0 ⋮

0

0



ð1 þ γ ðθ ÞÞð1  2vðθ ÞÞ

We use the Peaceman–Rachford and Douglas–Gunn scheme for the diffusion terms of Eq. (2.5a) and (2.5b) respectively and solve the reaction terms explicitly. This results in solving a tridiagonal system at each time step, which can be efficiently performed using Δt Thomas’ algorithm. Note that below we have denoted μx ¼ 2Δ x2 2 and δx U i;j ¼ U i  1;j  2U i;j þ U i þ 1;j . A.1.1. Peaceman–Rachford 

ðA:1Þ

   n þ 1=2 2 2 2 2 1  μx δx U ni;j;k ¼ 1 þ μx δx þ 2μy δy þ 2μz δz U ni;j;k þ f U^ i;j;k  2 2 n n 1  μy δy U nn ðA:2Þ i;j;k ¼ U i;j;k  μy δy U i;j;k 2 2 n þ1 ¼ U nn 1  μz δz U ni;j;k i;j;k  μz δz U i;j;k

It is well known that the eigenvalues of M are given by (see Noschese et al., 2013)

λk ðθÞ ¼ ð1 þ γ ðθÞÞð1  2vðθÞÞ þ 2vðθÞ cos

kπ ; nþ1

k ¼ 1; …; n

ðA:5Þ

and that the matrix of eigenvectors is given by the discrete sine transform matrix rffiffiffiffiffiffiffiffiffiffiffi 2 jkπ sin : ðA:6Þ Sj;k ¼ n þ1 nþ1 Since Sj;k is symmetric and orthogonal, we can solve (A.3) to obtain

A.1.2. Douglas–Gunn



C C C: C A

ðA:4Þ

A.1. Numerical scheme

Δt n 2 2 1  μx δx U ni;j ¼ U ni;j þ μy δy U ni;j þ f i;j 2  Δt n 2 2 1  μy δy U ni;jþ 1 ¼ U ni;j þ μx δx U ni;j þ f i;j 2

1

xt ðθÞ ¼ S diagðλðθÞÞSx0 ðθÞ: Assuming that x0 ðθÞ ¼ ð1; 0; 0; …; 0ÞT for all θ, then 0 n 1 X t πj πj sin λ ðθÞ sin B C nþ1 nþ1 C Bj¼1 j B C B n C BX t C π j 2 π j B C sin λ ð θ Þ sin 2 B j nþ1 nþ1 C xt ðθÞ ¼ Bj¼1 C: C n þ 1B B C ⋮ B C BX C πj nπ j C B n t @ A sin λj ðθÞ sin nþ1 nþ1 j¼1

ðA:7Þ

ðA:8Þ

Please cite this article as: Chan, M.H., et al., Mathematical modelling of spatial sorting and evolution in a host–parasite system. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.06.027i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

M.H. Chan et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Each entry is a Riemann sum and so for large t (since we require n ¼t for M to be a square matrix) the solution converges to 0 R 1 1 t 0 ðð1 þ γ ðθ ÞÞð1  2vðθ ÞÞ þ 2vðθ Þ cos ðπ xÞÞ sin π x sin π x dx BR C B 1 ðð1 þ γ ðθÞÞð1 2vðθ ÞÞ þ 2vðθÞ cos ðπ xÞÞt sin π x sin 2π x dx C B C xt ðθÞ ¼ 2B 0 C: B C ⋮ @R A 1 t 0 ðð1 þ γ ðθ ÞÞð1 2vðθ ÞÞ þ 2vðθ Þ cos ðπ xÞÞ sin π x sin nπ x dx ðA:9Þ References Baeumer, B., Kovács, M., Meerschaert, M.M., 2008. Numerical solutions for fractional reaction diffusion equations. Comput. Math. Appl. 55 (10), 2212–2226. Barton, D.P., 1998. Dynamics of natural infections of Rhabdias cf. Hylae (nematoda) in bufo marinus (amphibia) in Australia. Parasitology 117, 505–513. Bouin, E., Calvez, V., Meunier, N., Mirrahimi, S., Perthame, B., Raoul, G., Voituriez, R., 2012. Invasion fronts with variable motility: phenotype selection, spatial sorting and wave acceleration. C. R. Math. 350 (1516), 761–766. Brown, G.P., Phillips, B.L., Shine, R., 2014. The straight and narrow path: the evolution of straight-line dispersal at a cane toad invasion front. Proc. R. Soc. B: Biol. Sci., 281(1795). Dockery, J., Hutson, V., Mischaikow, K., Pernarowski, M., 1998. The evolution of slow dispersal rates: a reaction diffusion model. J. Math. Biol. 37 (1), 61–83. Douglas, J., 1962. Alternating direction methods for three space variables. Numer. Math. 4 (1), 41–63. Hastings, A., Cuddington, K., Davies, K., Dugaw, C., Elmendorf, S., Freestone, A., Harrison, S., Holland, M., Lambrinos, J., Malvadkar, U., Melbourne, B., Moore, K., Taylor, C., Thomson, D., 2005. The spatial spread of invasions: new developments in theory and evidence. Ecol. Lett. 8 (1), 91–101.

29 Hughes, C.L., Dythm, C., Hill, J.K., 2007. Modelling and analysing evolution of dispersal in populations at expanding range boundaries. Ecol. Entomol. 32 (5), 30 437–445. 31 Kelehear, C., Brown, G.P., Shine, R., 2012. Rapid evolution of parasite life history 32 traits on an expanding range-edge. Ecol. Lett. 15 (4), 329–337. Kelehear, C., Webb, J.K., Shine, R., 2009. Rhabdias pseudosphaerocephala infection 33 in bufo marinus: lung nematodes reduce viability of metamorph cane toads. 34 Parasitology 136, 919–927. 35 Kot, M., Lewis, M.A., van den Driessche, P., 1996. Dispersal data and the spread of invading organisms. Ecology 77, 2027–2042. 36 Noschese, S., Pasquini, L., Reichel, L., 2013. Tridiagonal toeplitz matrices: properties 37 and novel applications. Numer. Linear Algebra Appl. 20 (2), 302–326. 38 Perkins, T.A., Phillips, B.L., Baskett, M.L., Hastings, A., 2013. Evolution of dispersal and life history interact to drive accelerating spread of an invasive species. Ecol. 39 Lett. 16, 1079–1087. 40 Phillips, B.L., Brown, G.P., Travis, J.M.J., Shine, R., 2008. Reids paradox revisited: the 41 evolution of dispersal kernels during range expansion. Am. Nat. 172 (S1), S34–S48. 42 Phillips, B.L., Brown, G.P., Webb, J.K., Shine, R., 2006. Invasion and the evolution of 43 speed in toads. Nature 439 (7078), 803. 44 Phillips, B.L., Kelehear, C., Pizzatto, L., Brown, G.P., Barton, D., Shine, R., 2010. 45 Parasites and pathogens lag behind their host during periods of host range advance. Ecology 91 (3), 872–881. 46 Shine, R., 2010. The ecological impact of invasive cane Toads (Bufo Marinus) in 47 Australia. Q. Rev. Biol. 85 (3), 253–291. 48 Shine, R., Brown, G.P., Phillips, B.L., 2011. An evolutionary process that assembles phenotypes through space rather than through time. Proc. Natl. Acad. Sci. 108 49 (14), 5708–5711. 50 Simmons, A.D., Thomas, C.D., 2004. Changes in dispersal during species’ range 51 expansions. Am. Nat. 164, 378–395. Stover, J.P., Kendall, B.E., Nisbet, R.M., 2014. Consequences of dispersal hetero52 geneity for population spread and persistence. Bull. Math. Biol. 76 (11), 53 2681–2710. 54 Thomas, J., 1995. Numerical Partial Differential Equations: Finite Difference Methods. Number v.1 in Graduate Texts in Mathematics. Springer. Q3 55

Please cite this article as: Chan, M.H., et al., Mathematical modelling of spatial sorting and evolution in a host–parasite system. J. Theor. Biol. (2015), http://dx.doi.org/10.1016/j.jtbi.2015.06.027i

Mathematical modelling of spatial sorting and evolution in a host-parasite system.

There have been numerous empirical and agent-based modelling studies on the spatial self-structuring of traits, particularly in regard to dispersal ab...
2MB Sizes 0 Downloads 9 Views