Heterogeneity in Behaviour and Movement can Influence the Stability of Predator–Prey Periodic Travelling Waves

Cyclic predator–prey systems are often observed in nature. In a spatial setting, these can manifest as periodic traveling waves (PTW). Environmental change and direct human activity are known to, among other effects, increase the heterogeneity of the physical environment, which prey and predator inhabit. Aiming to understand the effects of heterogeneity on predator–prey PTWs, we consider a one-dimensional infinite landscape Rosenzweig–MacArthur reaction–diffusion model, with alternating patch types, and study the PTWs in this system. Applying the method of homogenisation, we show how heterogeneity can affect the stability of PTW solutions. We illustrate how the effects of heterogeneity can be understood and interpreted using Turchin’s concept of residence index (encapsuling diffusion rate and patch preference). In particular, our results show that prey heterogeneity acts to modulate the effects of predator heterogeneity, by this we mean that as prey increasingly spend more time in one patch type over another the stability of the PTWs becomes more sensitive to heterogeneity in predator movement and behaviour.


Introduction
Population cycles are one of the most studied aspects of populations dynamics (Barraquand et al. 2017), and such cycles can occur uniformly along the space that the species inhabit, or, alternatively, more complex spatiotemporal phenomena can be observed. One example of the latter happens in the red grouse population in north-east Scotland where population fluctuations in each location exhibited a decrease in synchronicity the further apart these locations were (Moss et al. (2000)). Another particular example studied by Lambin et al. (1998) and Berthier et al. (2014) for field voles in Kielder Forest in northern England, and in France, is a type of travelling wave observed in cyclic populations. In this case, populations in different spatial locations oscillate with the same temporal period, but exhibit a phase shift in space, giving the appearance of a wave with temporal oscillations in its wake. Such patterns are called periodic travelling waves (PTWs). PTWs happen when the populations move in a way that present periodicity both in space and time, with peaks and troughs of individuals moving across space at constant speed. A useful analogy is the wave-like phenomenon that spectators perform in a stadium during a crowded sport event, where they raise and lower their arms at similar moments to their neighbours in such a way that the net effect is a wave of raising and falling arms propagating around the stadium. Tenow et al. (2013) describes spatiotemporal data consistent with PTW for winter moth (Operophtera brumata) distributions in continental Europe. The study was able to estimate wavelength, speed and direction of the PTW, which exhibited decade-long periodic outbreaks. Sherratt and Smith (2008) review other past studies reporting populations possibly exhibiting PTW behaviour.
Spatiotemporal data is often challenging to obtain (Bennett and Sherratt 2017), since it involves a laborious and costly monitoring at several locations through usually extensive periods of time (Berthier et al. 2014). Hence, a motivation behind understanding the properties of PTWs comes from the possibility of improving our ability to track and predict the distribution of populations across space, e.g. in conservation field studies. Pest management can also benefit from an improved ability to track and control future outbreaks by increasing the efficiency of costly data collection (Petrovskii et al. 2014).
There are several frameworks for studying population dynamics in space, among which we highlight reaction-diffusion equations (Cantrell and Cosner 2004). They have been playing a fundamental role in modelling population movement since the work of Skellam (1951). Populations exhibiting PTWs can then be studied as solutions of reaction-diffusion models. Such solutions can be classified as either stable, with a regular spatiotemporal structure as described above, and unstable, which can manifest though irregular oscillations with no distinguishable shape (spatiotemporal chaos (Sherratt et al. 1995)). Populations displaying irregular spatiotemporal oscillations can be difficult to monitor and detect (Lambin et al. 1998). Therefore, by studying such reaction-diffusion models it is possible to better understand the underlying causes of spatiotemporal oscillations of natural populations. The seminal work of Kopell and Howard (1973) first determined the conditions for existence and stability of PTW solutions of coupled reaction-diffusion equations.
Amid the myriad of factors that may affect PTW stability, we investigate the effect of heterogeneity in the environment. In line with this thought, Johnson et al. (2004) argued how a mosaic-like heterogeneous environment alone can induce travelling waves in larch budmoth. There are many ways in which environmental heterogeneity affects populations. Different habitat types can be more or less suitable for population growth, and boundaries between habitat types can affect behaviour and movement (Schultz and Crone 2001;Bélisle and Desrochers 2002). For example, the grey shrike-thrush (Colluricincla harmonica) and the white-throated treecreeper (Cormobates leucophaeus) were found to act differently in forests when compared to how they act in open areas, even if they are adjacent, trying to avoid the latter moving through them more quickly (Robertson and Radford 2009). The predation risk of forest birds' nests by a nest predator (red squirrel, Tamiasciurus hudsonicus) can also vary spatially depending on habitat quality (from the point of view of the predator) (Martin and Joron 2003).
Nonetheless, there have been a small number of studies examining the effect of environmental heterogeneity specifically on PTWs. We are then motivated by the question left open by Sherratt and Smith (2008): "What is the effect of spatial heterogeneity on periodic travelling waves?" Sieber et al. (2010) showed that small temporal noise in reaction-diffusion equations is able to supress PTWs. In contrast, we are interested in considering spatial heterogeneity. Connected with such interest is the work of Sherratt et al. (2003), who studied heterogeneity generated by obstacles in the landscape which could not be traversed by the species. They found that the shape and size of obstacles affect the wavelength of PTWs in a two-dimensional reaction-diffusion predator-prey system.
However, spatial heterogeneity can be observed in different ways, such as in habitat patches, i.e. regions with relative homogeneity but significantly different from their immediate surroundings, being more or less suitable for each of the species inhabiting them. For example, Gauduchon et al. (2013) studied the repercussions of heterogeneity by comparing the dynamics of two predator-prey models in a habitat composed of three one-dimensional patches. The three patches differed by the species' parameters values from patch to patch. Their results showed that fragmentation via habitat loss can decrease cycle amplitude and average density of prey and predator populations. Shigesada et al. (1986) demonstrated that patch environmental heterogeneity in a single-species reaction-diffusion system could generate PTW solutions. Their analysis resulted in conditions for the population to either be able to invade the environment through a PTW or become extinct. Alternatively, Kay and Sherratt (2000) showed that environmental heterogeneity (through random spatial variation on the kinetic parameters) in a predator-prey system can allow persistence of PTWs that would otherwise die out in a finite domain. However, they did not consider the spatial heterogeneity in movement or the effects of habitat boundaries.
Many of the studies that consider heterogeneity via a collection of patches contemplate what we refer to as coarser grain environmental heterogeneity. By that we mean that patch size was not necessarily assumed to be small compared to the scale of movement rates. However, species can often encounter many habitat types during their lifetime and rapidly move through a landscape. Such a scenario can be thought of as finer grain heterogeneity, in which the effects of environmental heterogeneity can be studied through a landscape composed of a large number of small patches of different habitat types. One way to theoretically deal with such multiscale heterogeneity is via asymptotic homogenisation (Holmes 2012).
Our work is based on the previous studies of Yurk and Cobbold (2018) and Cobbold et al. (2022), who developed an asymptotic homogenisation framework for reaction-diffusion equations on a landscape composed of an infinite sequence of one-dimensional patches of two types. The patches differed in population dynam-ics, movement and the organism's response to patch boundaries. They studied how environmental heterogeneities on the patch level (e.g. in the scale of meters) have an impact on landscape level (e.g. in the scale of kilometres) population densities. The focus species are assumed to move through the environment in a way that each individual passes through a substantial number of different patches during its lifespan, behaving differently in each patch type. This framework essentially corresponds to finding a systematic way of "averaging" the patch level heterogeneity, resulting in an approximate but potentially much simpler system of equations that describes the aggregate effect of the many different patches. In particular, Cobbold et al. (2022), Yurk and Cobbold (2018) developed the homogenisation theory that could handle discontinuities in population density at patch interfaces. Such discontinuities arise when individuals show patch preference when they reach a patch boundary (Maciel and Lutscher 2013;Schultz 1998).
Hence, the present work aims to understand how environmental heterogeneity affects PTW solutions in predator-prey systems, using asymptotic homogenisation. The PTW solutions of the homogenised model were analysed through the same approach as Sherratt et al. (2003): the application of a normal form transformation as an approximation to small amplitude PTWs. In Sect. 2, we present the homogenised predator-prey model originally presented by Cobbold et al. (2022) and apply normal formal analysis to study the PTW solutions. The approach provides conditions for stability of the small-amplitude PTW solutions. In Sect. 3, we analyse the effects of heterogeneity on both the period and amplitude of spatially uniform solutions and on the stability of small amplitude PTW solutions. In particular, we consider how PTW stability is affected by the interplay between heterogeneity in behaviour and movement.

The Patchy Landscape Model
The model describes a predator-prey system inhabiting an infinite one-dimensional landscape composed of a sequence (indexed by i) of two alternating patch types, 1 and 2, with lengths l 1 and l 2 . The prey (u(x, t)) and predator (v(x, t)) population densities evolve in space (x) and time (t) according to the following reaction-diffusion equations where x i denote the boundaries of the patches. Since Eq. (1) describes the dynamics in each patch, we refer to them as the patchlevel equations. The diffusion rates of u and v inside patch i are constant within the patch and are denoted by Fig. 1 Example of the landscape structure composed of an infinite sequence of alternating patch types (1 and 2). The arrows represent the possibilities each prey individual has at the edge of the patch: it moves to patch type 1 with probability α u and to patch type 2 with probability 1 − α u D u i = D u 2 , and analogously for v. The same holds for all the parameters indexed by i. When a prey individual encounters the interface between the two patches, we assume it chooses patch type 1 with probability α u and patch type 2 with probability 1 − α u (analogously, probabilities α v and 1 − α v for a predator individual). In Fig. 1, we illustrate the patch structure. The interface conditions for the prey population density (u) are the same as used by Yurk and Cobbold (2018), which are where and x + i , x − i correspond to the right and left boundaries of patches i and i + 1, respectively. Conditions (2) are derived from the assumption that the flux of individuals must be constant across any two patches. Different assumptions over individual behaviour at patch interface in the underlying random-walk model (Ovaskainen and Cornell 2003) lead to different interface conditions. Our model assumes scenario 3 from (Maciel and Lutscher 2013), in which individuals choose each patch with different probabilities as described, but keep the step size constant once inside the patches. An analogous set of conditions are assigned to the predator population density v.
We assume the reaction terms on each patch in equation (1) follow the Rosenzweig-MacArthur model (Rosenzweig and MacArthur 1963): The Rosenzweig-MacArthur model is a predator-prey model with cyclic solutions for a well-known range of parameters (see Kot 2001). On each patch, the prey (u) grows logistically at a per capita rate λ i , with the strength of the intraspecific competition denoted by μ i . The predator (v) attacks the prey at rate a i and the consumption of prey saturates (for fixed v) as u increases, as a consequence of each predator individual expending a nonzero amount of time (the handling time h i ) to consume each prey individual. The prey-to-predator conversion coefficient in this prey consumption is γ , and the predator is removed from the system at a per capita rate m i (mortality rate).
The assumption that the landscape is periodic may seem relatively restrictive. However, from a sufficiently large scale, many real landscapes can be seen as collection of repeating homogenous patches (Fitzgibbon et al. 2001). Therefore, the periodic landscape described above can be a good approximation to many real ecological environments (Garlick et al. 2011). The framework utilised in this paper (asymptotic homogenisation, Sect. 2.2) is valid in more general heterogeneous environments. The periodicity assumption is, however, mathematically convenient, providing simplifications that allows for analytical results.

Asymptotic Homogenisation
The patchy landscape model described by equations (1), (4) and (2) is computationally expensive to numerically simulate for our case of interest, where we must consider large landscapes in order to observe PTW solutions. Asymptotic homogenisation (Holmes 2012) consists of averaging over the landscape heterogeneity, resulting in an approximated but much simpler version of the PDEs. We follow closely the work from Cobbold et al. (2022), where the homogenisation technique was developed for two-species reaction-diffusion systems on a patchy landscape, as presented in 2.1, with discontinuities at patch boundaries. The core underlying assumption for the asymptotic homogenisation to hold is that the patches are small enough that each individual is expected to go through a large number of patches during its lifetime. Through considerations involving the dynamic level of the populations (see Turchin 1998), it is possible to derive partial differential equations for the leading order of the power series expansion of the two populations' densities (leading orders named U for the prey, V for the predator). Cobbold et al. (2022) obtain the set of homogenised equations for U and V , withf and f 1 , f 2 given by Eq. (4). An analogous expression holds for g as an average of g 1 and g 2 . The "averaging" brackets correspond to the formula · = · 1 l 1 + · 2 l 2 l 1 + l 2 .
Residence index in patch type 1 Homogenised attack rate in patch type i Homogenised handling time in patch type î Averaged prey growth ratê Averaged predator mortality rate The parameter ρ u,v i , given by is called the residence index. At each given location in space, it can be interpreted as a quantity, which is proportional to the steady-state population density at that location if this population was only subject to diffusion. It is a "relative measure of the average time that an organism spends between entering and leaving a unit area" Turchin (1998). In our case, these unit areas are the individual patches of the landscape, and therefore, the residence indices (ρ u,v 1 , ρ u,v 2 ) are proportional to the average time each organism spends in each patch type. Equation (7) combines movement through diffusion (D u,v i ) and behaviour through patch preference (ξ u,v 1 = (1 − α u,v ) and ξ u,v 2 = α u,v ). If the residence index is large, it can be either because the individuals move slowly in that location or they have a high preference for that patch type.
The homogenised diffusion rates of Eq. (5) are given by: By expanding Eq. (6), the analogous expression forĝ and redefining the parameters according to Table 1, we obtain the following expressions for the reaction terms of Eq. (5):f The homogenised equations for U and V are similar to the ones at the patch level (1) for u and v, but due to spatial heterogeneity, have different reaction termsf and g. Indeed, in the limiting case where all patch level parameters have the same values between patch types 1 and 2 and there is no patch preference (i.e. the limiting case of a homogeneous landscape), the reaction terms simplify and become identical to the . In this limiting case, our model equations match the ones studied by Sherratt et al. (2003).
It is worth highlighting that the factors ρ u i ρ u l i l 1 +l 2 appearing as averaging weights in the expression for the homogenised parameters (Table 1) correspond to the relative amount of time that prey spends in type i patches if the system was only governed by diffusion . A completely analogous argument holds for the predator and the expression In order to be able to directly compare with the work of Sherratt et al. (2003) and to simplify the analysis, we perform a non-dimensionalisation of (5) using the following rescaled parameters: The resulting rescaled equations are: where h(X , T ) and p(X , T ) correspond to the rescaled prey and predator densities, respectively. Note that if η = 1, Eq. (11) simplifies and the reaction terms have the same functional form as the Rosenzweig-MacArthur model.

Normal Form Analysis
Equation (11) describes, with good approximation, the dynamics of the original patchy landscape model presented in Sect. 2.1, but are considerably easier to analyse and to numerically simulate. Equation (11) has up to 3 nonnegative spatially uniform equibria: extinction, prey-only, and coexistence. Cobbold et al. (2022) analysed the stability of these equilibria and showed that stable limit cycles exist. We are interested studying the region of parameter space where stable limit cycles occur. When our model is considered in the limiting case of a homogeneous landscape, the parameter used as a bifurcation parameter corresponds to the same used by Sherratt et al. (2003). Therefore, throughout this paper we use C as the bifurcation parameter to allow comparison between our work and Sherratt et al. (2003). There is a critical value of C (C = C crit (E 1 , η, β)) at which the kinetics undergo a Hopf bifurcation, with stable limit cycles for values of C above C crit . In the reaction-diffusion system, when C > C crit , Eq. (11) can exhibit PTW solutions and standard analysis of these solutions is possible. We follow the script provided by the work of Sherratt et al. (2003), performing a reduction to normal form Guckenheimer and Holmes (2013). In the case of δ = 1, the kinetics in Eq. (11) can be approximated by Hopf normal form giving equations of lambda-omega type (see, e.g. Murray 2001): whereh andp are nonlinear combinations of h and p (determined by the reduction to normal form ) and r 2 =h 2 +p 2 . With the aid of a computer algebra package, ω 0 and ω 1 can be written in terms of the homogenised parameters E 1 , B 1 , β, η, C − C crit . The full Mathematica notebook for the normal form transformation and subsequent stability analysis are provided at the open-access repository (Andrade and Cobbold 2022). The calculations involved in finding ω 1 (E 1 , B 1 , η, β) were performed with Mathematica (Wolfram Research 2019). All plots produced in this paper were made with Python's Matplotlib library (Hunter 2007). The reduction to normal form consists of a nonlinear transformation of Eq. (11) into a simpler set of equations that approximate the original equations when C is close to C crit , which we refer to as the small-amplitude regime. Thus, the stability of original PTW solutions can be inferred from the stability of the approximated PTW solutions that satisfy (12). Albeit limited to an approximation of small-amplitude solutions, the lambda-omega system provides insights into the dependence of the stability of PTW solutions on the ecological parameters. The work of Kopell and Howard (1973) shows that the lambda-omega system (Eq. 12) has a family of PTW solutions of the form: parameterised by R, the amplitude of the PTWs.
Typically PTWs are generated by either boundary conditions or population invasions. In order to study PTW stability analytically, we therefore consider a semi-infinite domain and we follow the approach from Sherratt et al. (2003), with zero-Dirichlet boundary condition at the origin. The asymptotic homogenisation and the reduction to normal form are still valid if the model described in Sect. 2.1 is semi-infinite. Sherratt et al. (2003) obtained a stability condition for the PTW solutions of the family (13) generated by Dirichlet boundary conditions if C is sufficiently close to C crit . The derivation is overviewed in Appendix A. The condition is that a PTW solution (13) generated by a zero-Dirichlet boundary conditions at the origin is stable if and only if Therefore, by transforming system (11) into its approximate version (12), we obtain an expression for ω 1 in terms of the homogenised parameters E 1 , B 1 , η and β. A different criterion holds for stability of PTWs generated by invasion (Sherratt 1998(Sherratt , 2001. Criterion (14) is then used to classify regions of the parameter space separating stable PTW and unstable PTW (spatiotemporal chaos). Figures 5, 7, 8, 9 and 10 in Sect. 3.2 are made using criterion (14) to classify each point of the paramater space as corresponding to a stable or an unstable PTW solution.
The reduction to normal form can only be applied to a particular form of Eq. (11), where δ = 1. This corresponds to requiring the diffusion coefficients in equation (5) to be identical (see Eq. 8): This constraint means that we have one less degree of freedom when choosing values for the parameters related to movement and space, We explore the biological implications of this constraint in the discussion.

Results
Before analysing the behaviour of the homogenised model, we compare the PTW solutions of the homogenised equations (5) to solutions of the patchy landscape model (1). Figure 2 illustrates that the homogenised equations provide a very close approximation. In Fig. 2 (left), we illustrate a stable PTW moving in the positive x-direction.
Changing parameter values leads to unstable PTW, giving rise to irregular oscillations and spatiotemporal chaos (Fig. 2 right). In the unstable PTW regime, for larger times the homogenised solution begins to deviate from the corresponding patchy landscape solution. This is a consequence of the spatiotemporal chaos associated with the unstable PTW, where small differences in the initial conditions due to the homogenisation approximation result in the long-term deviation of the two solutions (Sherratt et al. 1995;Sherratt 1996). The stable PTW solutions remain close for all times.

Cyclic Properties of the Spatially Uniform Solutions
Before studying the PTW solutions of (11), we examine how patch heterogeneity affects the stable limit cycle solutions of the kinetic ordinary differential equations, which correspond to the spatially uniform equivalent of Eq. (11). At a particular value of C = C crit , the non-trivial steady state undergoes a Hopf bifurcation and for C > C crit the non-trivial steady state (h * , p * ) is locally unstable and the kinetics have a stable limit cycle solution. C crit is determined by standard linear stability theory Initial conditions for the homogenised model are U (x, 0) = 4(1+sin(2π x/500)) for the prey and V (x, 0) = 2(1+sin(2π x/500)) for the predator. The corresponding prey initial conditions for the patchy landscape model were ρ u 1 U (x, 0)/ ρ u in patch type 1 and ρ u 2 U (x, 0)/ ρ u in patch type 2. Analogous conditions are used for the predator (V ). Boundary conditions: zero-Dirichlet at x = 0 and zero-flux at x = 500. Other parameters: The numerical simulations are carried out using the method described in Yurk and Cobbold (2018). The Strang splitting is applied, with diffusion terms implemented with Crank-Nicholson and fourth-order Runge-Kutta to update the kinetic step. For the patchy landscape model, derivatives that appear in the patch interface conditions were implemented using second-order forward or backward difference. Discretisation: (Glendinning 1994). The Hopf bifurcation occurs when the eigenvalues, λ * , of the Jacobian matrix (J ) of the kinetic ordinary differential equations evaluated at (h * , p * ) are purely imaginary, or equivalently when Solving this equation for C gives an expression for C crit in terms of η, β, E 1 , B 1 . The expression for C crit is algebraically messy (see section 1 of the Mathematica notebook (Andrade and Cobbold 2022)), but in the special case of η = 1 it simplifies to Close to the Hopf bifurcation (Tr[J (h * , p * )] ≈ 0), the temporal frequency of the limit cycle solutions is approximated by the imaginary part of λ * : √ det(J (h * , p * )). Therefore, the temporal period (τ ) of such solutions is approximated by: At η = 1, and C = C crit the temporal period simplifies to In the special case of η = 1, equation (18) demonstrates that for small amplitude cycles increasing B 1 increases period, while increasing E 1 has a nonlinear effect on period (Fig. 3). Our focus is on understanding the effects of heterogeneity, so we interpret how changes in E 1 and B 1 relate to patch differences. From expressions 10 we see that E 1 depends onm, the patch averaged predator mortality rate. Thus, increasing E 1 can correspond to decreasingm, or equivalently decreasing predator mortality rate on one or both patches. Sincem is a weighted average of the patch mortality rates, the weights (the proportion of time spent on each patch) control which patch most strongly influences the value ofm.
In addition to letting η = 1, if we assume predator mortality rate on both patches is the same (m 1 = m 2 ) and also assume that prey handling time is constant across the landscape (h 1 = h 2 ) then Eq. (18), and hence, cycle period is now independent of predator residence index (see Appendix B). By expressing η in terms of dimensional parameters as the expression for η gives us insights into the trade-offs between patch level parameters of the system. Consider an example in which predator attack rate varies between patch types. One might expect that if a predator spends more time in locations where attack rate is high this would increase total predator density and decrease prey density, impacting on population cycles, but this is not the case. Instead, the η = 1 constraint ensures that any difference in patch attack rates is balanced by differences in prey residence index. Prey spend less time on patches with high attack rates, so predators would gain no advantage by spending more time in high attack rate patches and therefore cycle period will be independent of predator residence index.
Hence, η = 1 presents a special case in which predator residence index does not affect population cycles, but prey movement does. In Fig. 3, we illustrate what happens as we relax this constraint on η. As η increases both period and amplitude increases, and for η = 1 predator residence index does affect population cycles (Fig. 4). Higher values of η correspond to larger patch differences in predator attack rate or handling time, which are not balanced by changes in prey movement. To fully understand these effects of heterogeneity, it is useful to consider the dimensional model (5) and the dimensional parameters (Table 1), allowing us to explicitly study the effects of prey and predator residence index, and patch size. The dimensional model also allows us to consider variation in patch prey growth rate, which we could not study in the non-dimensional model as average prey growth rate was used to scale time (see expressions 10).
In Fig. 4, we consider a scenario where prey growth rate on patch 1 (good patch) is greater than on patch 2 (bad patch), λ 1 = 10 > λ 2 = 5. We then vary the parameters associated with movement: l 1 /l 2 (the relative size of patch 1 compared to patch 2) and prey and predator residence index. Predator amplitudes are shown in Appendix C. From the plots in Fig. 4, three main results should be highlighted.
Firstly, the overall trend shown in Fig. 4 is that the amplitudes and periods of the predator-prey cycles increase with the proportion of 'good' patch sizes (increasing l 1 /l 2 ). This effect is particularly clear for small l 1 /l 2 and can be explained by the fact that, if the good patches occupy a smaller proportion of the habitat, prey spend more time in the bad patches, where prey growth rate is lower, leading to smaller population cycle amplitudes. For larger l 1 /l 2 , our results suggest that both the period and amplitude saturate at a fixed value, which varies depending on the value of prey and predator residence indices. The only exception to the overall trend is shown at the rightmost amplitude plot (ρ u 1 = 6.0, ρ u 1 = 2.0). For small l 1 /l 2 , increasing l 1 /l 2 can lead to a decrease in amplitude without an increase in period if ρ v 1 is close to ρ v 2 . This exception can be explained in the following way. Since ρ u 1 is bigger than ρ u 2 , prey spend more time in the good patches (patch 1). However, unless ρ v 1 is also considerably bigger than ρ v 2 , predator do not spend as much time in patch type 1 as the prey, allowing prey to increase its population to larger values. The phenomenon of a decrease in amplitude with l 1 /l 2 vanishes, as the proportion of good patches increases.
Secondly, we observe the effects of heterogeneity in predator residence index alone. In the centre and right plots of Fig. 4, the amplitude and periods of the cycles increase as we increase ρ v 1 and decrease ρ v 2 (solid to dotted to dashed to dashed-dotted lines). We obtain a similar qualitative result by just increasing ρ v 1 or just decreasing ρ v 2 (not shown). As the predator spends proportionally more time in the good patch (patch 1), where the prey also spends more time, there is an increase in the expected number of prey-predator encounters. Thus, predation is intensified and predator population grows to larger amplitudes. This causes prey oscillation amplitude to also be increased, with a corresponding increase in the period of the system.
Finally, we highlight the interplay between heterogeneity in prey and predator residence index. We find that the magnitude of the changes in period and amplitude is increased when there is heterogeneity in prey residence index (compare left and centre plots in Fig. 4). Moreover, the greater the heterogeneity in prey residence index, the greater the effect of heterogeneity in predator residence index (compare centre and right plots). In particular, if ρ u 1 = ρ u 2 , predator residence index has no effect on period  Fig. 4 Amplitude (upper) and period (lower) of prey oscillations in the spatially uniform solution of (11) as a function of relative patch size for different values of predator residence index. Patch type 1 is the 'good' patch type (λ 1 = 10) and type 2, the 'bad' patch type (λ 2 = 5). Patch sizes l 1 , l 2 satisfy l 1 + l 2 = 1. Left: ρ u 1 = 4, ρ u 2 = 4 (curves corresponding to each predator residence index pair are indistinguishable ). Centre: ρ u 1 = 5, ρ u 2 = 3. Right: ρ u 1 = 6, ρ u 2 = 2. Other parameters h 1 = h 2 = 0.5, γ = 0.9, m 1 = m 2 = 1.0, μ 1 = μ 2 = 1.75, a 1 = a 2 = 3.0. and amplitude of prey cycles. This indicates that the time populations spend in each patch type are an important factor in determining the properties of the population cycles.

Stability Boundaries of the Non-dimensional Homogenised Model
In this section, we consider the stability of the PTW solutions of Eq. (11). We show how heterogeneity could affect the stability of real biological systems exhibiting PTWs. In Fig. 5, we illustrate the stability regions for different values of η with respect to E 1 , B 1 , in order to directly compare to the results from Sherratt et al. (2003). The plot illustrates that stable PTWs are present for intermediate values of predator birth/death rates (E 1 ) and when the ratio of prey and predator maximum birth rate (B 1 ) is sufficiently high. Increasing the value of η (moving from the grey to the green and yellow curves) corresponds to increasing landscape heterogeneity by making the two patch types more distinct, either through increasing differences in prey residence index or through increasing differences in attack rate. If we consider the three biological examples discussed by Sherratt (2001), we find that increasing η shifts the hare-lynx system from a stable to an unstable PTW, and the zooplankton-phytoplankton system shifts from the unstable to the stable PTW region. On the other hand, the weasel-vole is predicted to remain in the unstable regime. This suggests that landscape heterogeneity could influence the dynamics observed in natural systems.
Such predictions are valid for small amplitude solutions. However, ecologists are generally interested in larger-amplitude PTWs. In Fig. 6, we illustrate the shift in stability with numerical simulations of the homogenised equations (11) for three amplitude values, measured in terms of the difference between the bifurcation parameter C Stability boundaries in terms of the dimensionless parameters E 1 , B 1 for different η, for β = 1.0 when C is close to C crit . The areas above the curves correspond to parameter regions where the PTW solutions are stable. The regions below each curve correspond to unstable PTW solutions. We plot with respect to B 1 (prey/predator maximum birth rate) and E 1 (predator birth/death rate) so that in the limiting case (η = 1.0 corresponding to spatial homogeneity) we obtain the corresponding stability boundary discussed by Sherratt et al. (2003). The crosses denote the parameter sets given by Sherratt (2001) for three example predator-prey systems: hare-lynx, zooplankton-phytoplankton and weasel-vole (Color Figure Online) and its critical value C crit ( = C − C crit = 0.01, 0.5, 1.0: upper, centre and lower plots, respectively). We consider the parameters of the zooplankton-phytoplankton system (E 1 = 1.8, B 1 = 4.0 (Sherratt 2001)) for two different values of η : 1 (left plots) and 7 (right plots). We find that as predicted from Fig. 5, the solutions for E 1 = 1.8, B 1 = 4.0 are unstable for η = 1 and stable for η = 7. Moreover, such result holds both for smaller ( = 0.01) and for larger amplitudes ( = 0.5, 1.0). However, the stability shift predicted for small amplitude does not necessarily hold for larger amplitudes for all points of the (E 1 , B 1 )-parameter space. In fact, the stability of larger-amplitude PTW solutions is still an open problem in the literature. Moreover, the direct classification of PTW stability through numerical integrations is additionally challenging since simulating Eq. (11) for E 1 , B 1 values close to the stability boundaries can require longer domains and simulation times to differentiate between stable and unstable solutions. This is illustrated in Appendix D, where it is shown plots for the numerical simulations for E 1 , B 1 of the hare-lynx and weasel-vole systems from Fig. 5.

Heterogeneity in Prey Residence Index
The remainder of the results section examines PTW solutions in terms of the parameters of the dimensional model (5). This allows us to isolate the effects of various sources of landscape heterogeneity (e.g. patch variation in movement, attack rates, etc.) that cannot be teased apart by looking at the non-dimensional model. We study the (λ −m)-parameter space, as these parameters appear in the expressions for E 1 and B 1 , respectively. Throughout the results section, we consider the case where prey growth rates and predator mortalities are identical between the two patches (λ 1 = λ 2 and m 1 = m 2 ), so thatλ andm do not depend on ρ u i and ρ v i . Therefore, we assume the only parameters which vary between patches are prey residence index (ρ u i ), predator residence index (ρ v i ) and attack rate (a i ). The normal form analysis constrains us to select prey and predator residence indices in such a way that the ratio of the homogenised predator and prey diffusion coefficients is 1 (Eq. 15). For example, if D v is fixed, increasingly higher values of ρ u 1 must be accompanied of, e.g. decreasingly smaller values of ρ u 2 in such a way thatD u has the same value for every pair (ρ u 1 , ρ u 2 ) considered. In this subsection, we consider the effects of prey residence index alone. In Fig. 7, it is shown that the PTW stability boundary shifts towards smaller predator mortalities and becomes narrower as heterogeneity in prey residence increases (increasing The shift can be interpreted in the following way: increasing ρ u 1 ρ u 2 means that the prey tend to spend more time in patch type 1. The predator, whose residence index is unchanged, encounters proportionally less prey in patch type 2. This causes lower predator growth rate across the landscape. Hence, for stable periodic travelling waves to be obtained, the predator needs to survive long enough in order to spend a sufficient amount of time in the locations where the prey has a larger residence index (patch type 1). Longer predator survival corresponds to smaller values ofm. Therefore, the stability boundaries shift towards smaller values of predator mortality rate as λ (prey growth) UNSTABLE Homogeneous ρ u 1 = 6.0, ρ u 2 = 2.0 ρ u 1 = 7.0, ρ u 2 = 1.0 ρ u 1 = 7.6, ρ u 2 = 0.4 Fig. 7 Stability boundaries for different values of prey residence index (coloured curves) in terms of the per capita prey growth rate (λ) and the per capita predator mortality rate (m), when the system is close to the Hopf bifurcation (when C is close to C crit ). The homogeneous boundary corresponds to do patch variation in prey residence index (ρ u 1 = 4.0, ρ u 2 = 4.0). The areas above the curves correspond to regions of parameters where the PTW solutions are stable, whereas the areas below the curves correspond to regions where the PTW solutions are unstable. λ 1 = λ 2 and m 1 = m 2 , henceλ andm do not depend on ρ v,u i . The other parameters are ρ v 1 = ρ v 2 = 4.0, a 1 = a 2 = 0.5, h 1 = h 2 = 0.5, γ = 0.9, l 1 = l 2 = 0.5. (Color Figure Online) Figure 8 shows stability boundaries for varying values of attack rate in patch type 2 (coloured curves). The main result is that heterogeneity in prey residence index can compensate for the effects of heterogeneity in attack rate. This is precisely the tradeoff discussed in Sect. 3.1. The different curves in Fig. 8 correspond to different values of attack rates in patch type 2 (a 2 ), whereas a 1 is fixed. In the left plot, prey residence index is the same on both patch types. In the right plot, prey residence index is larger in patch type 1 and smaller in patch type 2. If we increase prey residence index in patch type 1 enough, it is possible to counter act the effect of heterogeneity in attack rate, as we show in the right plot, where the green curve corresponds exactly to the homogeneous landscape stability boundary in Fig. 7 since a 2 ,a 1 ,ρ u 1 ,ρ u 2 are such that η = 1.

Heterogeneity of Attack Rate
The opposing effects of heterogeneity in attack rate and prey residence index can be interpreted as follows: if the predator attacks disproportionally more in one of the patch types, it has an equivalent effect to the prey occupying the other patch type proportionally more, where it is being preyed upon proportionally less often.

Heterogeneity of Predator and Prey Residence Indices
In Fig. 9, we show stability boundaries for varying predator residence indices (coloured curves) in patch type 2. The overall result is that strong heterogeneity in prey residence index increases the effect that heterogeneity in predator residence index has upon the location of the stability boundary. As before, prey and predator residence indices are selected in such a way that constraint (15)  UNSTABLE ρ u 1 = 7.00, ρ u 2 = 1.00 Fig. 8 The effect of varying the attack rate on patch type 2 (a 2 , coloured curves) in the stability boundaries for different combinations of prey residence index (ρ u 2 , ρ u 1 ). Plots with respect to per capita prey growth rate (λ) and the per capita predator mortality rate (m), when the system is close to the Hopf bifurcation (when C is close to C crit ). The areas above the curves correspond to regions of parameters where the PTW solutions are stable, whereas the areas below the curves correspond to regions where the solutions are unstable. λ 1 = λ 2 and m 1 = m 2 ; hence,λ andm do not depend on ρ v,u i . The other parameters are ρ v 1 = ρ v 2 = 4.0, a 1 = 0.5, h 1 = h 2 = 0.5, γ = 0.9, l 1 = l 2 = 0.5. (Color Figure Online) residence index has no effect on stability if prey residence index does not vary with patch type (other parameters being fixed). This is consistent with the result from Sect. 3.1. As we move from the left to right plot in Fig. 9, increasing heterogeneity in prey residence index, we notice an increased sensitivity to heterogeneity in predator residence index. The curves for different ρ v 1 , ρ v 2 pairs move further away from the ρ v 1 = ρ v 2 curve as ρ u 2 /ρ u 1 gets large. We saw a similar interaction between the effects of prey and predator residence indices in Sect. 3.1.
In the particular case where prey tend to spend equal time in both patches (Fig. 9, leftmost plot), the stability of the PTW solutions is unaffected by heterogeneity in predator residence index. This last result is explained from the fact that ω 1 does not depend on ρ v i (see Appendix B) when η = 1 and predator mortalities and handling times are identical between the two patches (m 1 = m 2 , h 1 = h 2 ).
In Fig. 10, we fix prey residence index so that it does not vary with patch type and attack rate is chosen to be larger in patch type 2 (a 2 > a 1 ). Here we highlight how heterogeneity in predator residence index interplays with heterogeneity in attack rate and predator mortality. In particular, heterogeneity in residence index can act to compensate or intensify the shift caused by heterogeneity in attack rate. The efficiency of the predator in the landscape depends on where it spends its time, how frequently it attacks prey in each location and how long it is alive. Hence, we see a trade-off of these quantities when we examine PTW stability boundaries. Decreasing ρ u 1 = 4.00, ρ u 2 = 4.00 ρ u 1 = 1.00, ρ u 2 = 7.00 Fig. 9 The effect of varying predator residence index ((ρ v 2 , ρ v 1 ), coloured curves) on the stability boundaries for different combinations of prey residence index (ρ u 2 , ρ u 1 ). Plots with respect to per capita prey growth rate (λ) and the per capita predator mortality rate (m), when the system is close to the Hopf bifurcation (when C is close to C crit ). The areas above the curves correspond to regions of parameters where the PTW solutions are stable, whereas the areas below the curves correspond to regions where the solutions are unstable. λ 1 = λ 2 and m 1 = m 2 , henceλ andm do not depend on ρ v,u i . The other parameters are a 1 = a 2 = 0.5, h 1 = h 2 = 0.5, γ = 0.9, l 1 = l 2 = 0.5. (Color Figure Online

Fig. 10
Stability boundaries for different levels of heterogeneity on predator residence index in terms of the per capita prey growth rate (λ) and the per capita predator mortality rate (m), when the system is close to the Hopf bifurcation (when C is close to C crit ). The areas above the curves correspond to regions of parameters where the PTW solutions are stable, whereas the areas below the curves correspond to regions where the solutions are unstable. λ 1 = λ 2 and m 1 = m 2 , henceλ andm do not depend on ρ v,u i . The other parameters are ρ u 1 = ρ u 2 = 4.0, a 1 = 0.5, a 2 = 3.5, h 1 = h 2 = 0.5, γ = 0.9, l 1 = l 2 = 0.5. (Color Figure Online) same predator efficacy could be obtained by a long-lived predator with a low prey attack rate, or a short-lived predator with a high attack rate. This last claim is supported by the shift from yellow to blue and red curves as we increase ρ v 2 /ρ v 1 . As the predators tend to occupy patch type 2 more frequently, increased predator mortality rates can maintain PTW stability (right shifting curves). In summary, stability can be controlled by patch differences in predator residence index.

Discussion
PTWs are often observed in ecological systems (Sherratt 2008) and were the focus of works that discussed their potential drivers (Kaitala and Ranta 1998;Petrovskii and Malchow 1999). As PTWs are an inherently spatial phenomena, it is natural to ask what are the implications of landscape alterations for such spatiotemporal dynamics. One way by which such alterations can happen is through the introduction of heterogeneity, which is expected to have an effect on ecological systems (Fahrig 2003). Environmental heterogeneity is known to affect PTWs (Shigesada et al. 1986;Sherratt et al. 2003;Kay and Sherratt 2000), but there is a lack of studies assessing the effects of heterogeneity that manifests through species movement and behaviour. Theoretical research can help to elucidate such effects and help identify when we expect to observe stable PTWs.
Our approach is the application of asymptotic homogenisation to study PTW solutions in a predator-prey model on a patchy landscape. Through a reduction to normal form, the stability of the PTW solutions is then analysed in terms of heterogeneity of residence index and attack rate. Our main result is that heterogeneity in prey and predator residence indices alone can shift biological systems from stable to unstable PTW regime and vice versa. This has implications for land management. For example, deforestation can change the proportion of forested to open areas. Many species show preference for forested habitats (Bélisle and Desrochers 2002), resulting in individuals spending more time in forested patches of land. This shift in residence index can be accompanied by a shift in the stability of PTWs. A given predator-prey exhibiting PTWs could see these destabilised resulting in irregular spatiotemporal oscillations. For the well-studied lynx-hare system (Turchin 2003;Krebs et al. 2018), our results suggest that introduction of environmental heterogeneity, such that hare spend more time in some locations than in others while lynx movement is unchanged, could result in the system being less likely to exhibit stable PTWs.
More generally, we contribute to the body of theoretical research that investigates the effects of environmental heterogeneity on predator-prey interactions (Ryall and Fahrig 2006;Stone and He 2007;Vitense et al. 2016). We showed that if prey spend disproportionally more time in one patch, the PTW stability boundary curve shift towards smaller predator mortalities. Such a shift can be reversed if the predator either spends more time or attacks more prey in the patches of high prey residence index. Our findings then illustrate the importance of the residence index as a useful tool to understand landscape-level phenomena. Since the residence index combines the effects of patch preference and diffusion, we stress the importance of patch-level properties on landscape-level patterns. The significance of patch properties has also been echoed by others, including Maciel and Lutscher (2013), who found that increased preference for a suitable (good) patch over an inhospitable (bad) surrounding patch decreased the minimal size of the good patch required for survival of the focal species.
The findings we obtain can be related to the work of Kay and Sherratt (2000) who considered PTWs generated by invasion on a finite domain. They showed that sufficient spatial noise in the kinetic parameters can allow the persistence of PTWs. However, if the noise is present at higher levels, irregular oscillations persist. In our work, we show that the heterogeneity in movement and individual behaviour can have complex effects on PTW stability, either fostering or suppressing stability of PTWs, depending on the values of the other kinetic parameters.
We show that the effects of predator residence index heterogeneity on PTW stability are increased as we increase the strength of heterogeneity in prey residence index. Our results suggest that a key factor affecting cyclic phenomena in general is the difference between the time prey and predator spend in each location. We argue that prey residence index should receive special attention in studies aiming to track and monitor predator-prey systems potentially exhibiting PTWs. We predict that if prey spend time relatively homogeneously across the landscape, heterogeneity in predator residence index should have a relatively small effect on the stability of the PTWs.
When considering the temporal cycles of spatially uniform solutions, we found a general increase in cycle period and amplitude with the relative good patch size. This result agrees with the overall trend obtained by previous reaction-diffusion predatorprey models (Maciel and Kraenkel 2014;Strohm and Tyson 2009). In a study involving coarse grain heterogeneity Gauduchon et al. (2013) showed that cycle amplitude in the middle of the patch decreases with predator patch preference for the good patch in the Rosenzweig-MacArthur model, but increases with predator patch preference in the May model. Despite our patch level equations having Rosenzweig-MacArthur reaction terms, our model shows a similar result to the May model studied by Gauduchon et al. (2013). The reasons for this mismatch remain unclear, but likely relates to the differing effects of coarser and finer grain heterogeneity on the Rosenzweig-MacArthur model.
Our analysis of the stability of the PTWs is constrained to the case where the homogenised prey and predator diffusion rates were equal. In biological systems, predator diffusion is typically larger than prey diffusion, so the δ = 1 constraint may appear overly constrictive. Indeed, the ratio between prey and predator diffusion is known to have an effect on the stability boundary of the PTWs in homogeneous systems (Smith and Sherratt 2007;Bennett and Sherratt 2017). However, δ = 1 can be satisfied while still allowing predator diffusion to be larger than prey diffusion at the patch level. For example, we can consider a scenario where the predator moves faster than the prey in both patches (namely, D u 1 < D v 1 and D u 2 < D v 2 ). This can be achieved while keeping δ = 1 if prey and predator prefer different patch types (see Appendix E). Differences in prey and predator habitat preference have been observed in the Vicuñas-Puma system (Smith et al. 2019) when prey chose habitat to avoid being attacked. Hence, while the δ = 1 constraint is restrictive it still allows the study biologically plausible scenarios. Future research with the use of alternative techniques would help elucidate the effects of relaxing the δ = 1 constraint. Among such techniques, we highlight numerical continuation, which can be used to numerically determine the essential spectra of reaction-diffusion operators (Rademacher et al. 2007). The stability of PTWs solutions can then be assessed through the study of such spectra (Sherratt 2013).
Our results are restricted to the scenario where asymptotic homogenisation is valid. Namely, when patch sizes are small compared to species diffusion. Other approaches would be necessary to investigate the effects of coarser grain heterogeneity and with larger patch sizes (e.g. as done by Gauduchon et al. (2013), Maciel and Kraenkel (2014), and Cobbold and Lutscher (2014)). However, even under this limitation, we are able to show that heterogeneity in movement and behaviour is able to affect PTWs in complex ways.
The greater part of our study is performed in the small-amplitude approximation. In general, ecologists are interested in larger-amplitude population cycles. Merchant and Nagata (2010), Merchant and Nagata (2015) developed approaches to study particular PTWs solutions of predator-prey systems further away from the Hopf bifurcation.
However, the analytical study of general PTW solutions for larger-amplitude oscillations remains an open problem in the literature and is therefore a possible step for future work. A natural question is whether our conclusions involving the effects of heterogeneity in movement and behaviour on PTW stability continue to hold for larger-amplitude PTWs. Nonetheless, the investigation of small amplitude is still used in the literature Sherratt 2017, 2019) as it provides useful analytical insights and limiting case predictions. Moreover, it produces a valuable framework to numerically investigate larger amplitude scenarios and a guideline to explore the parameter space through other non-analytical methods (e.g. numerical continuation).
Furthermore, another topic that deserves systematic study is the exploration of the effects of heterogeneity on PTW speed. Predicting the speed at which populations travel across the landscape is naturally a useful tool for monitoring and conservation studies but was beyond the scope of this work. As argued by Cobbold et al. (2022) and through our results, fine scale heterogeneity can lead to significantly different predictions from those found in a homogeneous environment. Therefore, one may expect that heterogeneity will also affect PTW speed.
In summary, our results shed light on factors affecting the stability of PTWs. In particular, we have shown that movement and behaviour in a heterogeneous environment can change the stability of PTWs in regions of parameter space where real-world examples of predator-prey dynamics can be found. We highlight that the strength of spatial variation in predator and prey residence indices can have a strong influence on PTW stability. Hence, future empirical research aiming to predict the effects of environmental heterogeneity (e.g. caused by human activity or climate change) on population spatiotemporal oscillations should take into consideration its repercussions on variation of individual movement and behaviour across the landscape.
where R corresponds to the amplitude of the PTW. Sherratt (2003) then showed that the boundary conditionh =p = 0 at x = 0 on Eq. (12) generates a particular solution of the family (13), in which R takes the value of Solutions under boundary conditionh =p = 0 at x = 0 correspond to solutions of Eq. (11) under Dirichlet boundary conditions h = h * and p = p * at X = 0, where (h * , p * ) is the steady-state solution. Therefore, combining (20) and (21), stability condition for the PTW solutions of Eq. (11) under Dirichlet boundary conditions h = h * and p = p * at X = 0 is which can be manipulated to obtain |ω 1 | < 1.110468. Finally, Sherratt et al. (2003) give solid evidence that such condition is valid for general Dirichlet boundary conditions at X = 0 provided the bifurcation parameter is sufficiently close to its critical value. Therefore, condition (14) determines the stability of solutions of Eq. (11) under zero-Dirichlet boundary condition at X = 0.
and the dependency of expressions (22) and (23) on ρ v 1,2 vanishes. In the special case of η = 1, the conditions for PTW stability also simplify. The reduction to normal form of Eq. (11) (with δ = 1) to Eq. (12) follows the script fully described by Sherratt et al. (2003) based on Guckenheimer and Holmes (2013). The first step in this script is to convert the linear part of the equations to normal form to give with H and P being linearly transformed versions of h and p. F(H , P) and G(H , P) determine the stability of the PTW solutions (13) and are obtained with the aid of Mathematica [see section 2 of Andrade and Cobbold (2022)]. The variable = C − C crit is assumed to be small (small-amplitude approximation). The expressions for 0 , 0 , 1 are computed from the eigenvalues of the Jacobian matrix associated with Eq. (11) evaluated at the coexistence steady state.

Appendix D: Extra Numerical Integrations of the Homogenised System
In this section, we exhibit numerical integrations for other points of the plot from Fig. 5, complementing the plots in Fig. 6. Figure 12 shows numerical simulations of the homogenised equations (11) with parameters values estimated for the hare-lynx system, predicted to be stable for η = 1 and unstable for η = 7. These plots illustrate how longer domains are necessary to determine the stability of the PTWs, as considering short domains can make unstable solutions to be misclassified as stable (compare the solution of Fig. 12 for η = 7, = 0.5 for X ∈ [0, 800] to the same solution for X ∈ [0, 7000]), both for small and large amplitudes. However, the simulations  Fig. 12 Numerical simulations of the homogenised equations (11) for the demographic parameters of the hare-lynx (E 1 = 1.2, B 1 = 4.4) system from Fig. 5. We consider the η = 1.0 (left plots) and η = 7.0 (right plots) cases, for different values of the bifurcation parameter C with respect to its critical value C crit : = C − C crit = 0.01 (upper plots), 0.5 (centre plots), 1.0 (lower plots). The solutions are plotted at T = 500,000. The same initial conditions, numerical scheme and discretisation from Fig. 6 were used. Boundary conditions: Dirichlet (h, p) = (h * , p * ) at X = 0 (steady state-dashed lines) and zero-flux at X = L = 7000. Other parameters: β = 1, δ = 1 (Color Figure Online (11) for E 1 = 4.0, B 1 = 10.0. We consider the η = 1.0 (left plots) and η = 7.0 (right plots) cases, for different values of the bifurcation parameter C with respect to its critical value C crit : = C − C crit = 0.01 (upper plots), 0.5 (centre plots), 1.0 (lower plots). The solutions are plotted at T = 500,000. The same initial conditions, numerical scheme and discretisation from Fig. 6 were used. Boundary conditions: Dirichlet (h, p) = (h * , p * ) at X = 0 (steady state-dashed lines) and zero-flux at X = L = 2500. Other parameters: β = 1, δ = 1 (Color Figure  Online) of longer domains require very long computation times in order to transient effects to vanish. Similar challenges arise when simulating the equations for the parameter values of the weasel-vole system (Fig. 13), predicted to be unstable for η = 1 and η = 7. Moreover, the (E 1 , B 1 ) pairs shown in Figs. 12 and 13 are close to the stability boundaries of η = 1 and η = 7 in Fig. 5. This causes the solutions to exhibit long transients, which require even longer simulations. In contrast, in Fig. 14 we simulate Eq. (11) for a point in the (E 1 , B 1 )-space further away from the stability boundaries. By comparing left and right plots in Fig. 14, there is clear shift in stability even for a shorter domain, matching the prediction from Fig. 5 both for small and large amplitudes: unstable for η = 1 and stable for η = 7.