Host--Parasitoid Dynamics and Climate-Driven Range Shifts


 Climate change has created new and evolving environmental conditions, impacting all species, including hosts and parasitoids. I therefore present integrodifference equation (IDE) models of host--parasitoid systems to model population dynamics in the context of climate-driven shifts in habitats. I describe and analyze two IDE models of host--parasitoid systems to determine criteria for coexistence of the host and parasitoid. Specifically, I determine the critical habitat speed, beyond which the parasitoid cannot survive. By comparing the results from two IDE models, I investigate the impacts of assumptions that reduce the system to a single-species model. I also compare critical speeds predicted by a spatially-implicit difference-equation model with critical speeds determined from numerical simulations of the IDE system. The spatially-implicit model uses approximations for the dominant eigenvalue of an integral operator. The classic methods to approximate the dominant eigenvalue for IDE systems do not perform well for asymmetric kernels, including those that are present in shifting-habitat IDE models. Therefore, I compare several methods for approximating dominant eigenvalues and ultimately conclude that geometric symmetrization and iterated geometric symmetrization give the best estimates of the parasitoid critical speed.

In the second stage, individuals disperse. For plants, dispersal occurs for seeds, while dispersers of other species like insects are adults. For each source, y, function k(x, y) is a nonnegative probability density function in x, describing the movement of the dispersers from y. All locations y ∈ [−L/2 + ct, L/2 + ct] contribute individuals that move to location x, according to probabilities determined by the dispersal kernel. The density at location x after dispersal is determined by integrating over source contributions. Individuals may be inside place.
In the first stage, which is sedentary, individuals within the habitat patch, [−L/2 + ct, L/2 + ct] undergo local processes, such as reproduction, death, birth, and growth. The growth or recruitment function f is chosen for a specific species or category of species to describe biological processes that happen in This model includes two discrete stages to map the population density at generation t, n t (x), to the population density at generation t + 1.
(1) n t+1 (x) = model species spread. Zhou and Kot's (2011;) IDE model for a shifting habitat is L/2+ct Since the early 2000s, researchers have been incorporating climate change velocity into mechanistic mathematical models that describe the growth and dispersal of biological populations. Potapov and Lewis (2004), Berestycki et al (2008), and Leroux et al (2013) each incorporated a moving habitat into continuous-time reaction-diffusion models. These models inherently assume that growth and dispersal occur simultaneously. To better model distributions for species that have distinct growth and dispersal stages, Kot (2011, 2013) introduced climate-driven range shifts in a discrete-time integrodifference-equation (IDE) model. IDE models also have the advantage of being able to incorporate a vast range of dispersal kernels to more accurately ecologists are using quantitative tools and mathematical models.
As a response to climate changes, suitable habitat ranges are shifting latitudinally and upwards in elevation (Brondizio et al, 2019;Chen et al, 2011;Parmesan, 2006;Parmesan and Yohe, 2003). Local extirpation and eventual global extinction of a species will result if populations of a species are unable to move with their shifting habitat. To better understand and predict the conditions that determine whether or not a species can persist as its range shifts,

Introduction
The earth is currently experiencing an unprecedented global rate of species extinction (Brondizio et al, 2019). Extinction is driven by numerous factors, including pollution, spread of invasive species, natural resource extraction, habitat fragmentation, and other changes due to anthropogenic land use. One of main drivers of extinction is climate change, and it has been shown that extinction will continue to accelerate with increased global warming (Brondizio et al, 2019;Urban, 2015). Ωn (2a) k n (x, y)F n [N t (y), P t (y)]dy, N t+1 (x) =

Model formulation
Consider interacting host and parasitoid populations on a one-dimensional habitat where there are discrete growth and dispersal stages, described by the IDE system geometric symmetrization and iterated geometric symmetrization.
Methods of analysis for the Case 1 system, which reduces to a single-species IDE, are explained in Section 3. Section 4 gives a brief summary of the impact of biological parameters on the critical speed for the Case 1 model. I then explain methods of analysis for the Case 2 model in Section 5. In Section 6, I compare a number of approximation methods and use geometric symmetrization to explore how biological parameters impact critical speed. I conclude with a discussion of the benefits and limitations of using spatially-implicit difference equations to analyze IDE systems and the potentially surprising accuracy of and I model population density of both species on the patch.
In Section 2, I describe the two models considered in this paper, including biological assumptions. The assumptions in Case 1 describe an ideal scenario for the parasitoid, such that the host density is at carrying capacity everywhere. In Case 2, the host and parasitoid have the same suitable habitat-patch, host-parasitoid system subject to climate change. For individual species and competing species on a moving habitat, there is a critical speed, c * , beyond which a population cannot keep up with the movement of its habitat (Harsch et al, 2014;Phillips and Kot, 2015;Rinnan, 2018b;Kot, 2011, 2013). However, the mathematical tools appropriate for single-species models are not sufficient for an IDE system with multiple species. I therefore explain a reduction from the spatially-explicit IDE model to a spatially-implicit system of difference equations. I then compare approximation methods for determining critical speed for persistence and conclude that an approximation known as geometric symmetrization is the best method for estimating the critical speed for parasitoid persistence. Using geometric symmetrization, I explore how different biological parameters impact host-parasitoid dynamics and persistence criteria.
need for further consideration of biotic interactions in IDE models.
In this paper, I describe and analyze an IDE model for the dynamics of a Biotic interactions are also important to consider when modeling species distributions in the context of climate change (Urban et al, 2013). Rinnan (2018b) recently used an IDE model to consider the possible outcomes for two competing species on a moving habitat, and Harsch et al (2017) stated the or outside the patch at the end of the second stage, but only individuals within the shifted patch will reproduce at the start of the next cycle. Figure 1 shows host and parasitoid local population dynamics occurring prior to dispersal with density dependence preceding parasitism in the life cycle of the host, in accordance with system (3). For a system with a forest tent-caterpillar host and a generic parasitoid, Cobbold et al (2005) describe how host adults disperse, reproduce, and lay eggs in midsummer, such that the host overwinters in the egg stage. On the other hand, parasitoids overwinter as pupae, and the adults disperse in spring to reproduce and lay eggs when they encounter and parasitize hosts. Cobbold et al (2005) also explain that parasitism of the forest tent caterpillar can occur "before, during, or after density-dependent host mortality, depending on the parasitoid species." The mathematical model analyzed by Cobbold et al (2005) assumes that parasitism precedes density dependence in the life cycle of the host and that parasitized hosts subsequently compete for resources. This differs from system where G[N t (y)] is the host per-capita-recruitment and H[P t (y)] is the fraction of hosts that succumb to parasitism per adult female parasitoid. The clutch size, c, is the average number of female parasitoid eggs laid on or in a single host that emerge and successfully become reproducing adults.
Ωn parasitism in the life cycle of the host. The model for this situation is A number of researchers have shown that population dynamics of hosts and parasitoids are impacted by the timing of density-dependent competition and parasitism in the life cycle of hosts (Cobbold et al, 2009;May et al, 1981;Wang and Gutierrez, 1980). To provide biological specificity for this analysis, I have chosen to focus on the scenario in which density dependence precedes is the probability that an adult from location y disperses to location x.
The sedentary stages for hosts and parasitoids are captured by the functions F n and F p , which describe intraspecific competition and interspecific parasitism between individuals at location y. Dispersal occurs after competition and parasitism, when the new adults of each species move from location y to location x. The density of adults of each species after dispersal at location x is the sum of adults from all possible contributing locations y that are at x at the end of the dispersal stage. In this model, the dispersal kernel, k i (x, y), For this model, let N t (x) be the density of reproducing host adults in generation t at location x, and let P t (x) be the density of adult female parasitoids in generation t at location x. Further, let Ω i be the patch of habitat that is suitable for species i. If an individual is outside the patch, that individual does not reproduce.
zero term of the negative binomial distribution, Common expressions for host density-dependent recruitment include the Beverton-Holt (Beverton and Holt, 1957) growth function, which is compensatory, and the Ricker (Ricker, 1954) growth function, which is overcompensatory. There is also a range of appropriate functions for parasitism, and it is common to represent the fraction of hosts that escape parasitism with the (3) for which parasitism follows host density-dependent competition. See also Marcinko and Kot (2020) for a more in-depth explanation of the local host and parasitoid interactions assumed in this paper.

5
Host-Parasitoid Dynamics and Climate-Driven Range Shifts Springer Nature 2021 L A T E X template where a is parasitoid searching efficiency (Cobbold et al, 2009;Livadiotis et al, 2016;May et al, 1981 As shown in Marcinko and Kot (2020), the simplest dynamics for the nonspatial host-parasitoid system arise from the model that uses Beverton-Holt H(P t ) = a and a fractional form of parasitism resulting from using κ = 1 in equation (4). For κ = 1, the fraction of hosts that succumb to parasitism per adult female parasitoid is To nondimensionalize, let n t (x) = N t (x)/K and let p t (x) = aP t (x). Define the nondimensionalized "clutch-size" as b = acK. Then per-capita-recruitment which pattern formation may occur.
For the purposes of this analysis, I use equation (5) for host densitydependent recruitment and equation (6) for parasitism. This choice will lead to an understanding of the basic persistence question and avoid systems in 1 + aP t In these expressions, R 0 is the host net reproductive rate, K is the host carrying capacity, and a is the parasitoid searching efficiency.
(6) . h[p t (y)] = 1 and the fraction of hosts that succumb to parasitism, p t (y)h[p t (y)], can be written with 1 + (R 0 − 1)n t (x) (7) , Following the same notation choices as in Marcinko and Kot (2020), system (3) can be rewritten as

Dispersal
Dispersal is an important area of study for both mathematical modeling and experimental understanding of species spread. Within the realm of possible dispersal kernels, the Laplace kernel has the benefit of being both mathematically elegant (Kot et al, 1996;Kot and Phillips, 2015;Kot and Schaffer, 1986;Lutscher, 2019;Reimer et al, 2016;Van Kirk and Lewis, 1997) and ecologically relevant as a leptokurtic distribution (Cobbold et al, 2005;Kot et al, 1996;Kot and Schaffer, 1986). Examples of analysis of insect dispersal data that exhibits leptokurtosis can be found in Aikman and Hewitt (1972); Bianchi et al (2009) Taylor (1978); Wallace (1966); Wolfenbarger (1946Wolfenbarger ( , 1959Wolfenbarger ( , 1975. While experimental dispersal kernels tend to exhibit fatter tails than those of the Laplace kernel, kernels with exponentially decaying tails are still widely used (Nathan et al, 2012).
The Laplace kernel is rightwards at speed c.
To investigate the host-parasitoid system with a shifting habitat, I consider two versions of system (9). The first case reduces the problem to a single IDE such that traditional IDE methods can be used. In the second case, I assume that the host and the parasitoid exist on the same habitat patch that is moving Writing this system with functions u and v will simplify the analysis needed in Section 5.1.
(11) v[n t (y), p t (y)] = bn t (y)g[n t (y)]h[p t (y)]. and Host-Parasitoid Dynamics and Climate-Driven Range Shifts 8 Springer Nature 2021 L A T E X template

Parameters
This paper seeks to consider parameter ranges that are biologically reasonable, but also broad enough to explore mathematical implications of the model. where δ is the mean dispersal distance and variance is σ 2 = 2δ 2 . The Laplace kernel only depends on the distance between x and y, and is an example of a kernel that describes homogeneous and isotropic dispersal. For examples and specifics in this paper, I assume a Laplace kernel. However, most of the following work can be easily numerically replicated with other kernels, including those with fatter tails.
(12) e −x-y/δ , where n t (x) is the nondimensionalized host density at location x and p t (x) is the nondimensionalized parasitoid density at location x. With Beverton-Holt 2.3 Case 1: Best possible scenario for parasitoid First consider a scenario in which a parasitoid population exists on a finite onedimensional patch of habitat with length L that is moving at speed c. Suppose also that the corresponding host population exists on a patch that is large enough to be approximated as infinite in comparison to the parasitoid habitatpatch. Additionally, assume that the host has a large enough net reproductive rate, R 0 , that its density is at carrying capacity throughout the infinite habitat. In some sense, this scenario represents an ideal set of circumstances for the parasitoid species since the host species exists at carrying capacity wherever the parasitoid encounters the host. As is common with IDE population models, I also assume that dispersal is homogeneous and isotropic so that dispersal depends on the distance between locations x and y. The corresponding IDE system is thus nondimensionalized by the use of L = 1.
As described by Cobbold et al (2005), space in an IDE model with a finite patch length can be scaled with respect to the patch length, L 0 , such that the scaled habitat-patch length is L = 1 and the scaled average dispersal-distance is given by δ = δ 0 /L 0 . I will primarily use L = 1 in the following sections. If the system has a true patch-length of L 0 = 1 km, scaled average dispersaldistances δ n and δ p are in km. Otherwise, the spatial parameters are effectively and σ p 2 = 0.241 km 2 . These parameters are based on data from Taylor (1995, 1997). For a slightly wider range of values, I will consider average host dispersal-distances from δ n = 1/30 km to δ n = 1/2 km and average parasitoid dispersal-distances as small δ p = 1/70 km. When it is relevant to display the full range of behavior of the results, I will consider δ p as large as δ p = 3.5 km.
Data on dispersal distances for insect hosts and parasitoids is limited, making it challenging to accurately estimate dispersal distance (Cobbold et al, 2005;Hughes et al, 2015). For a study on population dynamics of a forest tent caterpillar and a generic parasitoid, Cobbold et al (2005) estimated that host dispersal variance is between σ n 2 = 0.241 km 2 and σ n 2 = 0.377 km 2 for a Laplace distribution and that parasitoid dispersal variance is between σ p 2 = 0.060 km 2 rate, b, Hughes et al ( Then substitute this into the convolution equation (14) so that with h[p t (y)] from equation (8).
Since the habitat is assumed to be moving with constant speed, c, I look for a traveling pulse, Now, change variables such that x = x − ct and ŷ = y − ct, Then shift x by c, to get After dropping the hats for notational convenience, this equation is 2.4 Case 2: Host and parasitoid share a patch The second model assumes that both the parasitoid and host populations exist on a finite one-dimensional patch of habitat with length L that is moving at speed c. The IDE system for this case is L/2+ct Because this case reduces to a single IDE, I will analyze the persistence of the parasitoid population using techniques from Zhou and Kot (2011) in Section 3.
1 + p t (y) (19) dy.  (20) with c = 0 for 120 generations starting from an initial condition with both populations having a uniform distribution in the patch. The density for the host's uniform initial condition was 50, and the density for the parasitoid's uniform initial condition was 6. The t = 0 distributions in the figure represent the steady-state population distributions with no habitat movement. The population distributions for the moving habitat are shown every 10 generations. Note that the distributions shown are the adults after dispersal, but only those within the habitat patch Ω = [−L/2 + ct, L/2 + ct] reproduce at the next step. The parasitoid population decreases to a new steady-state distribution on the moving patch, while the lower population of parasitoids allows the host population to increase to a new steady-state distribution on the moving patch. If speed is sufficiently increased, the parasitoid population will die out. The convolutions at each iteration were computed using MATLAB's built-in 'conv' function with dx = 2 −10 .
where I again assume homogeneous and isotropic dispersal. Similarly to Case 1, look for solutions of the form Zhou and Kot (2011Kot ( , 2013 first developed methods for determining the critical speed from an IDE model for a single species on a moving habitat. Recall from Section 2.3 that the Case 1 model, system (13), reduces to the single IDE, equation (19). Therefore, the critical speed for parasitoid survival in optimal circumstances can be determined using methods from Kot (2011, 2013). Numerical simulations of this host-parasitoid system can either use system (20) or system (24). The first shows the population distributions from the observer's perspective, while the second shows changes in the population distribution from the moving-habitat frame of reference. In either case, a shifting habitat leads to a decrease in the total parasitoid population and an increase in the total host population as it experiences an ecological release from the parasitoid species. The qualitative shape of the population distribution on the habitat patch changes for both species, becoming narrower. The parasitoid population distribution skews towards the back end of the patch, as seen in Again dropping the hats for notational convenience, this system becomes Now use the change of variables x = x − ct − c and ŷ = y − ct, so the system is Substituting these into system (20) gives

Numerical approach
To find the critical speed for Case 1, I follow the numerical approach used by Kot and Phillips (2015); Kot (2011, 2013) in which I discretize equation (27), which yields a matrix eigenvalue problem. Specifically, I use Nystrom's method (Delves and Walsh, 1974;Lutscher, 2019;Press et al, 1992) with the trapezoid rule to write assumptions.
Furthermore, a positive kernel guarantees a simple, positive dominant eigenvalue, using Jentzch's theorem (Jentzsch, 1912). In this case, stability of the trivial solution will be lost as the dominant eigenvalue passes through λ = 1. The Laplace kernel satisfies the criteria for Jentzch's theorem. Thus, the critical speed for parasitoid persistence is the speed c * for which the spectral radius of the integral operator is one. Kot and Phillips (2015) give a nice overview of methods available for finding the critical speed under the previous discrete set.
While the general eigenvalue problem expressed in equation (27) is nasty, certain classes of kernels provide tractable solution methods. If we restrict the domain to x ∈ [−L/2, L/2] and y ∈ [−L/2, L/2] with finite L and use a continuous dispersal kernel with infinite support, then the linear integral operator is compact. In this case, the eigenvalues of equation (27) form a where λ is an eigenvalue of the integral operator with corresponding eigenfunction ξ(x). Equation (27) is a homogeneous Fredholm integral equation of the second kind.
Perturbations off the trivial solution are thus governed by Next, I look for solutions of the form to persistence of the parasitoid population. The linearization is Specifically, after assuming a traveling pulse and shifting to the movinghabitat frame, as was done in Section 2.3, I now linearize equation (19) around the trivial solution, p t (x) = 0. Instability of the trivial solution will correspond Springer Nature 2021 L A T E X template where ξ is an N ×1 vector of the densities and K is an N ×N matrix defined by an unrealistic ideal for the parasitoid species.
Figure 3 also shows that for larger average dispersal distances, increases in b have a larger impact on the critical speed for persistence. In all cases, the critical speed predicted by Case 1 will be an upper bound for the critical speed predicted by Case 2, since Case 1 represents a set of host assumptions that are dispersal, and the population cannot survive at all. Figure 3 shows the dependence of critical speed on average parasitoid dispersal-distance, δ p , for eight different values of b ranging from b = 2.2 to b = 3.6. Initially, the critical speed increases as dispersal distance increases, as dispersers are better able to make use of newly available habitat in front of the previous generation's patch. However, if average dispersal-distance increases too much, the population suffers larger losses from both dispersers left behind as the habitat patch advances and dispersers that overdispersed in front of the new habitat. Eventually, too many dispersers are lost from longer distance (2011,2013).

Case 1 results
Unsurprisingly, the critical speed for parasitoid persistence in Case 1 increases monotonically with increasing b. Note that this is a result previously shown by Zhou and Kot (2013) for the Laplace dispersal kernel, where b in equation (25) is effectively the net reproduction rate for the parasitoid species. Biologically, higher rates of reproduction for the parasitoid increase its ability to survive as the habitat shifts.
One can either use integral equation (27) or matrix system (29) to obtain the dominant eigenvalue, which depends continuously on the model parameters. For the results in Section 4, I used the method of bisection to find the critical speed, c * , which is the speed such that the spectral radius of matrix K, ρ(K), equals 1. At each step of the algorithm, I used MATLAB's built-in 'eigs' function to find the dominant eigenvalue, and then adjusted c accordingly. Other variations on these numerical techniques are discussed in Harsch et al (2014); Kot and Phillips (2015); Phillips and Kot (2015); Zhou and Kot 2 with k(·) as the Laplace kernel, equation (12).
A spatially-implicit approximation to IDE system (24) requires an approximation for the probability that a disperser remains inside the patch during the dispersal phase. Alternately, one can think of this same expression as the proportion of the population that remains in the patch after the dispersal (24). Cobbold et al (2005) analyzed a host-parasitoid IDE system with finite, stationary habitat by linearizing around the exclusion equilibrium, (n * (x), 0). They then used Van Kirk and Lewis' 1997 average dispersal success approximation for the exclusion equilibrium to determine a critical patch size. In the same study, Cobbold et al (2005) also used a spatially-implicit system of difference equations to approximate the spatially-explicit host-parasitoid IDE model. This second approach is applicable for a broader range of kernels, and informs my methods of analysis for the shifting-habitat IDE model, system

Case 2 methods
For studies of a single-species IDE model with moving habitat, the standard approach for determining critical speed involves linearizing around the trivial equilibrium and using an ansatz to generate an eigenvalue problem (Harsch et al, 2014;Kot and Phillips, 2015;Phillips and Kot, 2015;Rinnan, 2018a;Kot, 2011, 2013). Persistence of the population then corresponds to instability of the trivial equilibrium, under various assumptions. However, instability of the trivial or extinction equilibrium of an IDE system with two species does not correspond to persistence of both species, due to the existence of a state or states where one of the species persists and the other dies out. One can also derive two difference equations with the same form as system (31). Begin by denoting the averages of the two populations over the patch as where n t and p t are average host and parasitoid populations at generation t. Φ n is the proportion of host dispersers that remain in the patch after dispersal, and Φ p is the proportion of parasitoid dispersers that remain in the patch after dispersal.

Approximation methods
Heuristically, the average host and parasitoid populations on the patch [−L/2, L/2] in the moving-habitat frame of reference can be expressed as from simulations in Section 5.3.
In Section 5.1, I apply a spatially implicit approximation approach to derive criteria for parasitoid persistence. Van Kirk and Lewis' 1997 average dispersal success approximation has been shown to perform poorly in approximating IDEs with asymmetric kernels arising from a moving habitat Rinnan, 2018a). Therefore, I describe a number of possible approximations for the proportion of the population remaining in the patch after dispersal in Section 5.2. In order to compare the usefulness of these approximations, I then explain how the parasitoid critical speed can be determined stage. For IDEs with stationary patches, this concept has been used in a number of ecological contexts with great success (Cobbold et al, 2005;Fagan and Lutscher, 2006;Hughes et al, 2015;Lockwood et al, 2002;Lutscher and Lewis, 2004;Reimer et al, 2016;Van Kirk and Lewis, 1997).
Springer Nature 2021 L A T E X template L/2 following the approach in Fagan and Lutscher (2006); Lutscher and Lewis (2004); Rinnan (2018a,b). Note that for this paper, these averages refer to the average population for the patch in the moving-habitat frame of reference.
Next, compute the average of both sides of the equations for Case 2, system (24), Host-Parasitoid Dynamics and Climate-Driven Range Shifts Linearizing around (n t , p t ) and dropping higher-order terms yields the system, This system can also be expressed as In the context of a single-species IDE, average dispersal success, equation (37), is the Rayleigh-quotient approximation of the dominant eigenvalue of the critical speeds .
The difference equation approximation, system (35), is mathematically valid when the equilibrium population distributions for the host and parasitoid are close to the spatial average for all x. Even when the equilibrium distribution is not close to constant, the approximation tends to be relatively accurate, for cases with symmetric kernels (Lutscher, 2019;Van Kirk and Lewis, 1997).
The complicating aspect of moving-habitat IDE models is that the kernel is not symmetric. Reimer et al (2016) showed that the average dispersal success approximation, equation (37), decreases in accuracy as the kernel becomes more asymmetric. Equation (37) has also been derived from analytic approximations to the single-species moving-habitat model IDE using Legendre polynomials Zhou and Kot, 2013). Use of this approximation provides decent estimates of critical speed for some parameter regimes, but yields increasingly poor approximations for scenarios with higher S i is average dispersal success, as defined by Van Kirk and Lewis (1997).

Fixed points and stability
Recall that Φ n and Φ p represent the proportion of the population that remains in the patch after dispersal. For any approximation of Φ n and Φ p , we can determine persistence criteria from difference equation system (31). Assume compensatory density-dependence and that parasitism has the fractional form of equation (8). Then host and parasitoid persistence can be determined from the fixed points and stability of the work in this paper.) Analytic approximations include the average dispersal success approximation, S, Reimer et al.'s 2016 modification to the average dispersal success approximation, Ŝ, and the geometric success function, G, from Rinnan (2018a,b) and Kot and Phillips (2015). Numerical approximations begin with discretizing integral operator (38) to obtain a matrix K i , where I use the discretization given in Section 3.1. The dominant eigenvalues of K n and K p will be approximated with maximum and minimum row-sums, iterated row-sums from Yeh (1997), and iterated geometric symmetrization from Kot and Phillips (2015). These approximations are described in detail in Section 5.2 and are used to approximate Φ n and Φ p to obtain results in Section 6. (Note that the power method is an alternative numerical approach to approximating the dominant eigenvalue and is worth future consideration, but is beyond the scope of in the context of a two-species IDE moving-habitat model. Because Case 2 IDE system (24) has asymmetric kernels, I compare a number of approximations to the dominant eigenvalue of the linear integral operator (38), specifically using the Laplace kernel evaluated at x + c − y to account for moving habitat. Spectral radii can either be determined using the integral operator in equation (38) or a matrix discretization of the integral operator, previously described in Section 3.1. The following approximation methods were all compared by Kot and Phillips (2015) or Rinnan (2018a) for a single-species moving-habitat IDE, and have not previously been compared The validity of this approximation depends on the kernel being symmetric, as discussed by Kot and Phillips (2015).
Host-Parasitoid Dynamics and Climate-Driven Range Shifts 18 Springer Nature 2021 L A T E X template The three fixed points of system (39) are the extinction equilibrium, (0, 0), the exclusion equilibrium, (n, p) = Φ n and the coexistence equilibrium, The dynamics of system (39) are determined by the parameters R 0 , b, and the dispersal-survival proportions, Φ n and Φ p . The extinction equilibrium, (0,0), is asymptotically stable if Φ n R 0 < 1. The exclusion equilibrium, point The coexistence equilibrium, point (41), is asymptotically stable if Φ n R 0 > 1 and bΦ p = stability. Note that Φ n and Φ p will depend on c, L, δ n , and δ p . To find the critical speed, c * , for the parasitoid population, one can numerically solve Φ n (R 0 − 1) Parasitoid extinction occurs when the coexistence equilibrium loses stability and the exclusion equilibrium gains stability. The host dies out when the exclusion equilibrium loses stability and the extinction equilibrium gains Details for the stability analysis are given in Appendix A. Based on the nonspatial difference equation previously studied in Marcinko and Kot (2020), I assume R 0 > 1 and b > 1. Additionally, note that 0 < Φ n < 1 and 0 < Φ p < 1 by definition.

Dominant eigenvalue approximations
Each of the following describes an approximation to the dominant eigenvalue of the integral operator (38).
(44) Φ n R 0 − 1 for fixed R 0 , b, L, δ n , and δ p . For each approximation method for Φ n and Φ p , I compare the resulting parasitoid critical speed to the critical speed found from simulations in Section 6.1.
Springer Nature 2021 L A T E X template |Ω| Ω Ω

Geometric symmetrization
Following Kot and Phillips (2015) and Rinnan (2018a,b), I define geometric symmetrization of the average dispersal success as G = 1 is For short, I will refer to this simply as "geometric symmetrization." For the Laplace kernel on a habitat interval of length L centered at 0, this The geometric symmetrization with the Laplace kernel evaluates as (55) e − 2δ x+c-y e −y+c-x/(2δ) dxdy. The primary benefit of this approximation is that the kernel in equation (54) is symmetric and more robust to increased climate velocity. Kot and Phillips (2015) and Rinnan (2018a) showed that geometric symmetrization provides a relatively tight lower bound for the critical speed for a single-species moving-habitat IDE model. Lutscher (2019) gives a heuristic explanation for why geometric symmetrization provides a good approximation to the dominant eigenvalue of the original IDE, but also notes that finding a rigorous connection between the original IDE problem and the geometric symmetrization is an

Iterated maximum and minimum row-sums
One can also find increasingly tight bounds on the spectral radius of nonnegative K beginning with a result from Minc (1988). For a nonnegative matrix with nonzero row-sums, where r(x) is the redistribution function, equation (48). The analytic version of this method has a computational advantage over the numeric version if the integrals can be computed analytically. However, the next method builds on the matrix version of maximum and minimum row-sums, so I use the numeric version described by (58) to build towards the iterated method discussed next. is The analytic equivalent of these bounds for the spectral radius of the original integral operator K from equation (38) with the habitat-shifted kernel where ρ(K) is the spectral radius of nonnegative matrix K. This is a classic result from Frobenius (1912).

Maximum and minimum row-sums
The simplest bounds on the spectral radius of nonnegative matrix K are based on the row sums of the matrix. Let r i be the ith row sum of K, defined as N Host-Parasitoid Dynamics and Climate-Driven Range Shifts 22 Springer Nature 2021 L A T E X template where D (m) is the diagonal matrix with the ith row sum of matrix K (m−1) as the ith entry on the diagonal of D (m) . From this sequence of matrices, define the sequence Springer Nature 2021 L A T E X template g ij = where • represents the elementwise-product of the two matrices (the Hadamard product) and K T is the transpose of K. Equivalently, the elements of G are

Iterated geometric symmetrization
The analytic expression for geometric symmetrization, equation (54) was originally based on geometric symmetrization of an asymmetric matrix . For an asymmetric nonnegative matrix K, Schwenk (1986) defined the matrix G as the unique matrix such that Yeh's method were needed to get reasonable bounds for the critical speed. Kot and Phillips (2015) used Yeh's iterated row-sums with a numerical root finder to attain upper and lower bounds for the critical speed for the singlespecies moving-habitat IDE. For certain parameter values, several iterates of such that the sequences {min i r i (m) } and {max i r i (m) } converge, providing tighter bounds on ρ(K) as m increases. Yeh (1997) proved that for each m, computationally faster to compute for the Laplace kernel than ρ(G).
While ρ(G) provides a lower bound on the spectral radius of K, the analytic geometric symmetrization expression from equation (54) is significantly (66) ρ(G) ≤ ρ(K).
I will refer to G as the matrix geometric symmetrization. Importantly, G is symmetric and (65) k ij k ji . e T Ge vector e = (1, 1, ..., 1) T . This sequence is Kolotilina (1993) built on Schwenk's result and provided a sequence of increasingly tight lower bounds on ρ(K), based on the Rayleigh quotient with Springer Nature 2021 L A T E X template N ≤ e T G 2 e N ≤ · · · ≤ e T G 2 m e 1/2 Figure 5 zooms in on a range of c values near the sharp drop-off. Based on simulations like this one for various sets of parameters, one can reasonably conclude that if the parasitoid's maximum density is below 10 −3 after 1000 iterations, then the population will die out asymptotically. Technically, this definition of critical speed is an upper bound on the critical speed from is very different from the maximum density after 4000 iterations. Figure 4 shows that as c increases, the parasitoid's maximum density initially decreases slowly. For sufficiently small c, the parasitoid's maximum density after 200 iterations is indiscernible from the parasitoid's maximum density after 4000 iterations. As c continues to increase, the parasitoid's maximum density drops off sharply, and the maximum density after 200 iterations a set of biological parameters are shown in Figures 4 and 5.
To explore how many generations or iterations are necessary to determine whether the parasitoid population survives for a given c, I ran simulations for a range of c values and calculated the parasitoid's maximum density after 200 iterations up to 4000 iterations, in intervals of 200 iterations. The results for

Simulations
The Case 2 IDE system can be numerically simulated using either system (20) from the observer's perspective or system (24) from the moving-habitat frame of reference. In the moving-habitat frame of reference, the parasitoid population reaches a steady state distribution. As c increases for a given set of biological parameters, the total population of the parasitoid's steady state decreases. For some speeds, the steady state of the parasitoid distribution is very small. For slightly larger speeds, the parasitoid population slowly declines over time and asymptotically dies out. To determine the critical speed from simulations, it is necessary to distinguish between asymptotic survival at small values and asymptotic extinction.
previously mentioned approximations.
For the purposes of this paper, I will compare the use of iterated geometric symmetrization, (67), in difference-equation system (39) with the use of N where G 2 m is the geometric symmetrization of K 2 m . As a note of caution, this expression is misstated in Kot and Phillips (2015) as being based on powers of matrix G. The correct expressions are in terms of the geometric symmetrization matrices of powers of K, as written in Kolotilina (1993); Rinnan (2018a). Rinnan (2018a) also provides a sequence of analytic geometric symmetrization expressions based on recursively-defined kernels with increasingly many nested integrals. In this section, I compare the parasitoid critical speeds from each of the methods discussed in Section 5. I also compare the critical speed predicted by Case 1 with the critical speed for Case 2. Finally, I discuss the impact of biological parameters on the parasitoid critical speed.
For all critical speeds from simulations in Section 6, I used a bisection method to find the smallest c for which the parasitoid's maximum density was below 10 −3 after 1000 iterations. Numerical simulations for these results used the moving-habitat frame, system (24), for simulations, which is much more computationally efficient than simulating with the observer's frame of reference, system (20). Numerically simulating with system (24) restricts the values of c to integer multiples of the grid spacing, dx. For the simulations in Section 6, I used dx = 2 −11 ≈ 4.88 × 10 −4 . The integrals at each iteration were computed using MATLAB's built-in 'conv' function. For additional discussion of simulations, but it provides a significantly more precise metric than has been previously used (Rinnan, 2018b).  For all parameter ranges, geometric symmetrization resulted in a lower bound for the parasitoid critical speed that qualitatively matched the critical speed from the simulation method as parameters were varied. The lower bound from geometric symmetrization is not particularly tight, but it is tight enough to be useful as a conservative estimate on the maximum climate velocity that computational advantage as the other analytic methods compared here.
For complete disclosure, the integrals for the analytic version of the rowsum method from the inequalities in (59) can be evaluated for the Laplace kernel. Using those expressions to compute critical speed would have the same (1997) iterated row-sums, and iterated geometric symmetrization could not be computed on a laptop due to insufficient memory. Each of these methods nevertheless ran faster than simulations on an external server.  Note that geometric symmetrization gives a lower bound for the critical speed. The critical speed from the maximum rowsum method has huge error, but is an upper bound. Yeh's iterated maximum and minimum row-sums method gives an upper and lower bound, though the geometric-symmetrization method gives a tighter lower bound than two iterations of Yeh's iterated minimum row-sums.  Fig. 8: A comparison of the approximated critical speed from geometric symmetrization (solid) and one iteration of iterated geometric symmetrization (dotted). The critical speed from simulations is shown with connected markers. Biological parameters were L = 1, R 0 = 2, and δ p = 1/15 for both figures, δ n = 1/5 for the left figure, and b = 3 for the right figure. Both geometric symmetrization and iterated geometric symmetrization give lower bounds for the critical speed. Iterated geometric symmetrization is a tighter lower bound, but is more computationally and memory intensive.
6.2 Critical speed for parasitoid survival: Compare Case 1 as an upper bound for Case 2 While geometric symmetrization and iterated geometric symmetrization provide nice lower bounds for the parasitoid critical speed, I am also interested in comparing a few upper bounds. In Figure 9, I compare critical speed from Case 1 with the critical speed from two iterations of Yeh's iterated maximum row-sums. Since Case 1 describes ideal circumstances for the parasitoid, the critical speed from Case 1 is an upper bound on the critical speed for Case 2. For context, I also include the predicted critical speeds from Reimer et al.'s modification and geometric symmetrization. For sufficiently small δ p , Case 1 is the best upper bound for most sets of parameters, and it does especially well as R 0 increases. However, for even moderate values of parasitoid average dispersal-distance, δ p , the Case 1 approximation vastly overestimates the critical speed. speed.
Overall, geometric symmetrization from equation (54) is the most useful approximation for Φ n and Φ p to determine the parasitoid critical speed. It provides a relatively tight lower bound and qualitatively matches the critical speed from simulations for varied parameters. For kernels that do not have an explicit expression for geometric symmetrization, iterated geometric symmetrization is the best option, providing a tighter lower bound on the parasitoid critical geometric symmetrization for the Laplace kernel.
Finally, in Figure 8, I show some of the comparisons of geometric symmetrization and one iteration of iterated geometric symmetrization. Both approximations provide lower bounds on the parasitoid critical speed, and iterated geometric symmetrization gives a tighter bound. Iterated geometric symmetrization is based on the discretization from Section 3.1. For sufficiently small dx, the method cannot run on a laptop with 16 GB of memory. It is also much slower than determining critical speed using the analytic expression for geometric symmetrization is still the preferable approximation method. Figure 7 compares the parasitoid critical speeds determined by using maximum matrix row-sums, two iterations of Yeh's iterated maximum and minimum row-sums, and geometric symmetrization. Approximating dominant eigenvalues with the minimum row-sum predicted parasitoid death for c = 0, so this method was useless. The maximum row-sum approximation always drastically overestimated parasitoid critical speed. The critical speeds resulting from two iterations of Yeh's iterated maximum and minimum row-sums provided upper and lower bounds on the critical speed. The lower bound from iterated minimum row-sums is tighter than the upper bound from maximum iterated row-sums. Geometric symmetrization produces a tighter bound than the lower bound from two iterations of Yeh's minimum iterated row-sums, so geometric symmetrization qualitatively captures the effect of each biological parameter on the critical speed. Critical speed for parasitoid persistence increases monotonically for increasing R 0 . Higher host net reproductive rate, R 0 , means that the host population 6.3 Critical speed for parasitoid survival, coexistence: Using geometric symmetrization Using geometric symmetrization and equation (44), one can study the impact of growth parameters R 0 and b and average dispersal distances δ n and δ p on the parasitoid's ability to persist in a moving habitat. Figure 10 shows critical speeds for parasitoid persistence for a large range of parameter combinations.  is better able to overcome increasing losses corresponding to increasing climate velocity because more individuals are produced within the patch. Higher host density in the patch provides increased oviposition opportunities for the parasitoid, directly resulting in higher numbers of parasitoid offspring and increasing the parasitoid population's ability to overcome losses due to higher climate velocity. From Figures 10b, 10c, and 10d, it does appear that the benefits of higher R 0 are capped, which makes sense in comparison to Figure  9a where Case 1 provides a horizontal line as an upper bound for parasitoid critical speed, independent of the value of R 0 .
As long as parasitoid average dispersal-distance is relatively small, Case 1 provides a decent upper bound on the critical speed for parasitoid persistence for Case 2, IDE system (20). For most sets of biological parameters, however, a more accurate prediction of parasitoid persistence requires modeling the hostparasitoid interactions in a model for which space is implicit, using difference-critical speed increases as δ p increases.

Discussion
Dynamics of real-world populations are significantly impacted by biotic interactions with other species and available habitat spatially. Integrodifferenceequation models have the capacity to model species interactions while explicitly incorporating space, including climate-driven habitat shifts. At the same time, moving-habitat IDE models are rarely analytically tractable, and model analysis must rely on numerical approaches and approximations. Case 1 is a model that reduces the two-species host-parasitoid interactions to the dynamics of only one of the species while explicitly accounting for space. Analysis of Case 2 takes the opposite approach: modeling both the host and parasitoid species but implicitly accounting for space. Both approaches can be useful.
The crossing lines for larger values of host dispersal-distance in Figure 10a are consistent with a result from Cobbold et al (2005) that for higher host average dispersal-distance, the parasitoid must disperse shorter distances than the host to survive. For shorter host dispersal-distances towards the left of Figure 10a, the parasitoid gains an advantage from dispersing further, and the survival decreases.
For both host and parasitoid average dispersal distances, critical speed initially increases with increased average dispersal-distance and then decreases. This suggests that increased dispersal-distance is initially beneficial to the population, but that beyond a certain point, hosts and parasitoids can overdisperse. If the host population overdisperses, the critical speed for parasitoid to benefits to parasitoid persistence as b increases.
Critical speed for parasitoid persistence also increases monotonically for increasing b. Recall from Section 2 that while b is effectively a parasitoid growth-rate parameter, it is defined as the product of parasitoid searching efficiency, parasitoid clutch size, and host carrying capacity, b = acK. Biologically, more efficient parasitoid searchers are better able to find available hosts and persist on faster-moving habitat patches. Larger clutch size, c, directly increases the parasitoid's ability to reproduce and persist. Similarly to R 0 , larger host carrying capacity, K, increases hosts available for oviposition and parasitism. Unlike the results for R 0 , my explorations show no apparent limit The spatially-implicit approximation techniques used in this chapter nevertheless have many avenues for continued study. Both parasitism and density-dependent effects on the host were modeled with relatively simple functional forms, equations (8) and (7), corresponding to aggregated parasitism and compensatory density-dependence. For the nonspatial host-parasitoid models studied in Marcinko and Kot (2020), Ricker growth (Ricker, 1954) and/or an exponential form of parasitism led to wide array of possible dynamics, including population cycles and chaos. It would be interesting to consider how the type of host density-dependence or the level of parasitoid aggregation impacts the persistence of the coexisting host and parasitoid populations. Certain insect interactions with other species.
While spatially-implicit models have widespread use, a lingering question remains: What is lost by approximating a spatially-explicit system by a spatially-implicit model? In this paper, I have explored the conditions for which the parasitoid species persists by looking at a difference-equation system for the dynamics of average host and parasitoid populations. Simulations reveal that averaging over the patch loses some of the interesting behavior demonstrated by simulations. For example, under certain conditions, habitat patch movement leads to a higher maximum parasitoid-density than the maximum density on the stationary habitat patch, as seen in Figure 11. A spike or outbreak in the parasitoid population that moves through a landscape with the moving habitat could potentially have consequences for ecosystems and of habitat fragmentation on cyclic population dynamics.
Approximating the dynamics of a spatially-explicit model with a spatiallyimplicit model is certainly not a new idea. Spatially-implicit differenceequation models have been used to approximate IDE models with a network of domain patches (Fagan and Lutscher, 2006;Reimer et al, 2016), habitat fragmentation (Lutscher and Lewis, 2004), stage-structured populations (Fagan and Lutscher, 2006;Lutscher and Lewis, 2004;Reimer et al, 2016), alongshorecurrents (Reimer et al, 2016), a host-parasitoid system (Cobbold et al, 2005), a host-parasitoid system in a patchy landscape (Hughes et al, 2015), and two competing species on a moving habitat (Rinnan, 2018b). In the complementary realm of continuous-space models, 2012 reduced a PDE model for predator and prey dynamics to a system of spatially-implicit ODEs to explore the impact computational cost.
If a more accurate estimate for critical speed is needed for Case 2, a tighter lower bound can be achieved with iterated geometric symmetrization, though at a computational cost. For an upper bound, Case 1 is appropriate for some circumstances or Yeh's iterated maximum row-sums can be used, also at a a parasitoid species could withstand. Initial explorations further suggest that geometric symmetrization is flexible enough to permit modifications to the host and parasitoid habitat assumptions, similar to Rinnan's work (2018b) for two competing species. All of the results in this paper assumed a Laplace kernel for dispersal of both the host and the parasitoid. For some insect species, kernels with fatter tails than the Laplace kernel have experimentally been shown to better approximate dispersal, especially when individuals can disperse long distances (Nathan et al, 2012). The results in Sections 6 were based on average dispersaldistances for hosts and parasitoids assuming a Laplace dispersal kernel. For a fixed average dispersal-distance, it could be informative to consider the impact to persist as its habitat shifts due to climate change?
Alternately, for koinobiont parasitoids where parasitism occurs during early host larval stages, emergence time of the parasitoid larvae impacts population dynamics (Cobbold et al, 2009). In this paper, I assumed that density dependence precedes parasitism in the life cycle of the host. Does the relative timing of parasitism in the life cycle of the host affect a parasitoid population's ability species exhibit periodic cycles and outbreaks (Varley et al, 1974), and the interaction of spiking and dropping population densities over time could certainly complicate a population's ability to persist as its habitat shifts. Fig. 11: A simulation of the population distributions from Case 2 where the host and parasitoid occupy the same habitat that is moving to the right at speed c. These simulations used the same host and parasitoid parameters and set-up as Figure 2, but with c = 0.1 instead of c = 0.115. The population distributions for the moving habitat are shown every 10 generations. After the patch begins moving, the maximum parasitoid-density increases. The total parasitoid population of the steady-state distribution on the moving patch is smaller than the total parasitoid population prior to patch movement.
which gives the exclusion equilibrium, (n * , p * ) = Φ n R 0 − 1 (A1) 1 = Φ n u(n * , 0), Appendix A Appendix: Fixed points and stability for Case 2 system of difference equations The fixed points for the difference-equation approximation for Case 2, system (39), are found by setting n t+1 = n t = n * and p t+1 = p t = p * . Extinction of both species is an equilibrium solution for system (39). Additionally, if p * = 0, then the equilibrium is found by solving equation (39a) for n * , ability to persist despite climate-driven habitat shifts.
Finally, this chapter provides additional support for the use of geometric symmetrization for spatially-implicit modeling of IDEs with asymmetric kernels, specifically for two interacting species. Even though an explicit mathematical connection between geometric symmetrization and the dominant eigenvalue of the integral operator is currently an open question (Lutscher, 2019), geometric systematization provides a valuable tool for modeling species' of tail fatness on the critical speed for host and parasitoid coexistence, building on the single-species results from 2013. Host-Parasitoid Dynamics and Climate-Driven Range Shifts Springer Nature 2021 L A T E X template (n * , p * ) = Φ n is the coexistence equilibrium. For compensatory density-dependence and clumped or aggregated parasitism, the solution to this system is (A3a) Finally, the solution to the system 1 = Φ n u(n * , p * ), 1 = Φ p v(n * , p * ), where the subscripts on u and v denote partial derivatives. The partial derivatives of u and v are given in the first column of Table A1.
J(n, p) = (A5) , Φ p (pv p + v) Φ p pv n Φ n (nu n + u) Φ n nū p To determine the stability for each of the three fixed points, I also need the Jacobian of system (39), which is bΦ p + (R 0 − 1)Φ n .
(A4) Table A1: Partial derivatives used to simplify Jacobian expressions for the exclusion equilibrium and the coexistence equilibrium.
At this equilibrium, it also holds that Φ n g(n * ) = 1 from equation (A1). Use this to simplify the Jacobian to Springer Nature 2021 L A T E X template g(n * ) = 1 Further simplify by evaluating g(n * ) at point (A2) and simplify. Then .
Since R 0 and Φ n are positive, this criteria is most succinctly expressed as (A12) and |λ 2 | < 1. For the first eigenvalue, this criteria is When inequality (A13) is satisfied, λ 2 is positive. Further, satisfaction of inequality (A13) also implies that R 0 > 1 since 0 < Φ n < 1. Thus, stability criterion (A14) reduces to (A14) For the second eigenvalue, the stability criteria is where parameters pertaining to the parasitoid species are on the left and parameters pertaining to the host species are on the right. Φ n R 0 − 1 (A15) , n * p * R 0 − 1 A.3.1 Jury conditions for stability of the coexistence equilibrium Using the expressions from Table A1, the first Jury condition, inequality (A19) simplifies to Satisfaction of these Jury conditions guarantees asymptotic stability of the coexistence fixed point. A detailed discussion of the meaning of the Jury conditions for a similar host-parasitoid system can be found in Marcinko and Kot (2020).

A.3 Stability of the coexistence equilibrium
Now consider the coexistence equilibrium, point (A4). The Jacobian is Host-Parasitoid Dynamics and Climate-Driven Range Shifts 38 Springer Nature 2021 L A T E X template R 0 g(n * )h(p * ) + h(p * ) g(n * ) This is the same first Jury condition as Model 1 in my nonspatial work (Marcinko and Kot, 2020). If n * > 0, p * > 0, this reduces to h(p * ) R 0 n * (A22) > 0. n * (A23) > 0. τ = 2 + Φ n n * u n + Φ p p * v p = 2 − R 0 − 1 and Kot (2020), consider the sign of τ , This condition is satisfied whenever the coexistence fixed point is in the interior of the first quadrant. In the next section, I determine the criteria for the coexistence point to be in the interior of the first quadrant. Now look at the second Jury condition. Using a strategy from Marcinko 1 + 1 − R 0 to if n * > 0, p * > 0, and R 0 > 1. Finally, consider the third Jury condition, inequality (A21), which simplifies The expression for τ in equation (A24) is the same as τ for Model 1 in my previous work (Marcinko and Kot, 2020) and was shown to be positive for n * > 0, p * > 0, and R 0 > 1. Therefore, the second Jury condition is a sufficient condition for the first Jury condition for coexistence equilibrium (A4) and the first Jury condition is a sufficient criterion for the second Jury condition.
The n-coordinate of the coexistence equilibrium, point (A4), is which is true for R 0 > 1 for the equilibrium in the interior of the first quadrant.
The coexistence equilibrium value of the parasitoid for the averaged system