## Abstract

Many natural populations undergo multi-year cycles, and field studies have shown that these can be organized into periodic travelling waves (PTWs). Mathematical studies have shown that large-scale landscape obstacles represent a natural mechanism for wave generation. Here, we investigate how the amplitude and wavelength of the selected waves depend on the obstacle size. We firstly consider a large circular obstacle in an infinite domain for a reaction–diffusion system of ‘*λ*–*ω*’ type. We use perturbation theory to derive a leading order approximation to the wave generated by the obstacle. This shows the dependence of the wave properties on both parameter values and obstacle size. We find that the limiting values of the amplitude and wavelength are approached algebraically with distance from the obstacle edge, rather than exponentially in the case of a flat boundary. We use our results to predict the properties of waves generated by a large circular obstacle for an oscillatory predator–prey system, via a reduction of the predator–prey model to normal form close to Hopf bifurcation. Our predictions compare well with numerical simulations. We also discuss the implications of these results for wave stability and briefly investigate the effects of obstacles with elliptical geometries.

## 1. Introduction

Many natural populations are cyclic with multi-year oscillations in population density. In recent years, attention has focused on the possibility of a spatial component to these oscillations, and a number of spatio-temporal field studies have shown that they are organized into periodic travelling waves (PTWs). Examples include field voles in Northern UK (Lambin *et al*. 1998; MacKinnon *et al*. 2001) and in Fennoscandia (Ranta & Kaitala 1997), red grouse in Northeast Scotland (Moss *et al*. 2000) and larch budmoth outbreaks in the European Alps (Bjørnstad *et al*. 2002; Johnson *et al*. 2004). PTWs are a generic solution type in self-oscillatory systems, and in the case of oscillatory reaction–diffusion equations, there is a large body of literature spanning the last three decades (e.g. Kopell & Howard 1973; Ermentrout 1981; Ermentrout *et al*. 1997; Romero *et al*. 2000; Blowey & Garvie 2005; Garvie 2007). In particular, any oscillatory reaction–diffusion system has a one-parameter family of PTWs, with wave speed, wavelength or amplitude being possible parameters. In the context of the PTWs seen in ecological field data, there are two central questions to be addressed by mathematical models. Firstly, what is the underlying mechanism that drives the system to PTWs rather than, say, spatially homogeneous oscillations? And secondly, which member of the PTW family is selected by this mechanism?

Sherratt *et al*. (2002, 2003) have shown that a natural potential mechanism for PTW generation is the presence of large-scale landscape features for which the appropriate boundary conditions are of Dirichlet type. For PTWs generated in this way in real ecological systems, a critical issue is the effect of obstacle size on wave selection. Sherratt *et al*. (2003) performed a numerical study of this for one particular cyclic predator–prey model, showing that wavelength decreased with obstacle size, approaching a limiting value as the size tended to infinity. Our objective in this paper is a systematic study of this size dependence. We consider the following reaction–diffusion system of ‘*λ*–*ω*’ type(1.1a)(1.1b)with *ω*_{0} and *ω*_{0}−*ω*_{1} having the same sign. The amplitude . These equations are the normal form of an oscillatory reaction–diffusion system with scalar diffusion near a supercritical Hopf bifurcation, and as such our results can be expected to generalize to a wide range of oscillatory reaction–diffusion systems. In ecological applications, the variables *u* and *v* reflect the deviation of population densities from a coexistence steady state; hence *u* and *v* can be positive or negative. PTW solutions of (1.1*a*) and (1.1*b*) have been studied extensively since their first description by Kopell & Howard (1973). The wave family has the simple form(1.2a)(1.2b)where the amplitude *r*^{*} lies between 0 and 1.

The simplest context in which to consider the effect of obstacle size on PTW generation is a circular obstacle in an infinite domain. Provided that the initial conditions have circular symmetry, the whole solution will do too, and will satisfy(1.3a)(1.3b)where *x* (≥0) is the distance from the centre of the obstacle. The two terms with spatial derivatives in (1.3*a*) and (1.3*b*) come from converting the spatial diffusion terms in (1.1*a*) and (1.1*b*) to polar coordinates.

Defining the obstacle radius as *X*, we consider the boundary condition *u*=*v*=0 at *x*=*X*. In ecological applications, this would correspond to populations being fixed at a coexistence steady state. Such a condition can be appropriate in situations where there is an abrupt change in habitat; a full discussion of this is given in §3. For numerical solutions, we also need a boundary condition at *x*=*x*_{∞}, a suitably large domain radius, and we use *u*_{x}=*v*_{x}=0. Numerical solutions of (1.3*a*) and (1.3*b*) subject to these boundary conditions show ‘target pattern’ waves throughout the domain, except for a thin region close to the *x*=*x*_{∞} boundary (figure 1); the waves can move either towards or away from the obstacle, depending on the parameters *ω*_{0} and *ω*_{1}. Close to the obstacle, the wave properties are influenced by their curvature, but at large values of *x* the waves are essentially linear, and thus the amplitude at large *x* represents a simple means of comparing solutions for different values of the obstacle radius *X*. Our numerical solutions show that the wave amplitude decreases as a function of *X* (figure 2). As *X*→0 the amplitude approaches 1, corresponding to spatially homogeneous oscillations (the wavelength→∞); this is the behaviour that would develop when there is no obstacle. Conversely as *X*→∞, the amplitude approaches a limiting value that we denote by *a*. Note that the generation of PTWs in solutions such as that illustrated in figure 1 depends crucially on the Dirichlet boundary condition at the obstacle edge. Altering this to a Robin condition, for example, leads to a different PTW being selected (Sherratt submitted), while for a zero-flux (Neumann) condition, the long-term behaviour is spatially uniform oscillations.

The results illustrated in figure 2 mirror those reported by Sherratt *et al*. (2003) for a cyclic predator–prey model, suggesting that this variation in wave properties with obstacle size may be a general feature of oscillatory reaction–diffusion systems. In practice, the waves generated by small obstacles are not of interest in applications, since their wavelengths would be large compared with typical ecological domain sizes. For example, parameter estimates for the interaction between field voles (*Microtus agrestis*) and weasels (*Mustela nivalis*) suggest that linear obstacles generate a periodic wave of wavelength approximately 16 km (Sherratt *et al*. 2003), while an obstacle of radius 1 km generates a wave of wavelength approximately 25 km. In comparison, the largest field vole–weasel habitat in the UK (Kielder Forest) is approximately 30 km across. This illustrates the general conclusion that a PTW will be undetectable in practice unless the obstacle generating it is large.

In this paper, we study the case of large obstacle radius in detail for the *λ*−*ω* system (1.1*a*) and (1.1*b*). Rewriting (1.3*a*) and (1.3*b*) in terms of polar coordinates *r*−*θ* in the *u*−*v* plane gives(1.4a)(1.4b)In terms of these new variables, the PTW solutions (1.2*a*) and (1.2*b*) are simply(1.5)The simplicity of this new formulation means that the *λ*–*ω* system can be studied much more easily than most oscillatory reaction–diffusion equations. In particular, Sherratt (2003) derived a form for the PTW amplitude generated by the zero-Dirichlet boundary condition in the limiting case of *X*=∞ (i.e. a linear obstacle). In §2, we will extend this to derive the rate at which the wave amplitude approaches this limiting value when *X* is large but finite.

## 2. Wave selection for large obstacle radius in the *λ*–*ω* equations

Detailed investigation of the numerical solutions of (1.1*a*) and (1.1*b*) subject to *u*=*v*=0 at *x*=*X* shows that at large times, *r* and *θ*_{x} approach steady-state solution profiles. Therefore, we look for solutions of (1.4*a*) and (1.4*b*) in which *r* and *θ*_{x} are functions of *x* only. The latter implies that , where *κ* is a constant of integration. Substituting this solution form and *r*=*R*(*x*) into (1.4*a*) and (1.4*b*) give a third-order system of ODEs. Since we are focusing on large values of the obstacle radius *X*, we define *ϵ*=1/*X*. We also make two other notational changes that simplify subsequent algebra: and . The equations for the steady-state solution profiles then become(2.1a)(2.1b)where , with the boundary condition at the obstacle edge(2.2a)Since we are looking for solutions that approach a PTW far from the obstacle, we also require(2.2b)for some value of the wave amplitude *A*. Substituting this behaviour at infinity into (2.1*a*) and (2.1*b*) implies that the wave amplitude *A* is related to the constant via , and for convenience we use *A* rather than henceforth.

In the case of a linear obstacle (*ϵ*=0), Sherratt (2003) has shown that the large-time solutions for *R* and *ϕ* have the simple forms(2.3)whereWe will use this solution as the basis for a perturbation theory study of (2.1*a*) and (2.1*b*). For algebraic convenience, we will use *a* rather than *ω*_{1} as a parameter in subsequent calculations. Equation (2.3) implies that 0<*a*≤1. In the remainder of the paper, we will assume that *a*<1, excluding the special case *a*=1, which corresponds to *ω*_{1}=0. This case can be dealt with separately in a straightforward manner. When *ω*_{1}=0 , (2.1*b*) implies thatwhere *μ* is a constant of integration. The boundary condition (2.2*a*) then implies that *μ*=0, i.e. *R*^{2}*ϕ*=0. From (2.2*b*), it then follows that either *A*=0 or 1. The former would mean that *R*≡0, the spatially uniform equilibrium, which is unstable and thus not a potential large-time solution. Hence *A*=1, independent of *ϵ*. Therefore when *ω*_{1}=0, the PTW has the degenerate form of a spatially uniform oscillation for any obstacle radius.

Our aim in the remainder of this section is to calculate the wave amplitude to first order in *ϵ*≠0 when *a*<1, in the form . We require different expansions of (2.1*a*) and (2.1*b*) depending on whether or 1=*O*(*ϵy*), and we refer to these as the near-field and far-field solutions, respectively. The expansion of the far-field solution will allow us to determine appropriate end conditions for the near-field solution. We will then use the near-field solution to determine the leading order correction, *A*_{1}, to the wave amplitude.

### (a) Perturbation analysis of far-field solution

We begin by considering (2.1*a*) and (2.1*b*) when 1=*O*(*ϵy*); this scaling enables us to track the effect of boundary curvature at *y*=0 forward to *y*=∞. We substitute the rescaling into (2.1*a*) and (2.1*b*), and also , which follows from (2.3). This gives(2.4a)(2.4b)where we use bars to distinguish the far-field solution variables from the near-field solution. We now expand , and *A* as power series in *ϵ* as follows:(2.5a)(2.5b)(2.5c)Substituting the expansions (2.5*a*), (2.5*b*) and (2.5*c*) into (2.4*a*) and (2.4*b*) givesfor the *O*(1) terms, andfor the *O*(*ϵ*) terms. (Recall that we are assuming *a*<1.) These solutions satisfy the required conditions (2.2*b*) as *y*→∞ with , as demanded by the leading order solution (2.3).

One immediate consequence of the form of the leading order correction , to the solution in the far field is that the amplitude *R* approaches its limiting value *A* very slowly as *y*→∞. For fixed *ϵ* (note that we have not yet determined *A*_{1}),(2.6)This algebraic rate of approach is in marked contrast to the case of a linear obstacle (*ϵ*=0), for which (2.3) implies an exponential approach to the periodic wave amplitude: as *y*→∞. This difference is clearly visible in the numerical solutions of (1.3*a*) and (1.3*b*) on sufficiently large domains, as illustrated in figure 3.

### (b) Perturbation analysis of near-field solution

We now consider (2.1*a*) and (2.1*b*) when . Again we expand *R* and *ϕ* as power series in *ϵ*, and in this case the leading order solution is known from Sherratt (2003), which are as follows:(2.7a)(2.7b)Substituting the expansions (2.7*a*), (2.7*b*) and (2.5*c*) into (2.1*a*) and (2.1*b*) gives(2.8a)(2.8b)for the *O*(*ϵ*) terms. These equations have similarities to those studied by Sherratt (submitted), who used perturbation theory to investigate a variant on the boundary condition at the edge of a linear obstacle. Following Sherratt (submitted), we convert (2.8*a*) and (2.8*b*) into a single third-order equation for *R*_{1} by differentiating (2.8*a*) with respect to *y* and eliminating *ϕ*_{1}, giving(2.9)where(2.10)The third-order equation (2.9) can be reduced to second order using the derivative of the leading order solution. Substituting gives(2.11)In appendix A, we show that the corresponding homogeneous equation(2.12)has linearly independent solutions(2.13)Here *α*, *β* and *γ* are functions of *a* given in equations (A 2*a*), (A 2*b*) and (A 2*c*) in appendix A, and *F* is the hypergeometric function. (For a review of this special function, see ch. 15 of Abramowitz & Stegun (1970).) A general solution of (2.9) is therefore given by(2.14)Here, *K*_{w}(*a*) is defined by ; it is given explicitly as a function of *a* in equation (A 4) in appendix A.

### (c) Matching

To proceed further, we require a boundary condition on the near-field solutions *R*_{1} and *ϕ*_{1} as *y*→∞. This is found by matching the near-field and far-field solutions. We rescale *y* and such that both can be expressed in terms of an intermediate variable, *z*,(2.15)Expanding for small *ϵ* in the near-field and far-field solutions with *z*=*O*(1) gives(2.16a)(2.16b)Therefore, matching requires that and as *y*→∞. In appendix A, we show that the first of these conditions is satisfied if and only if *C*_{1}=*C*_{2}=0 and the second then follows immediately from (2.8*a*).

### (d) Behaviour near *y*=0

To consider the boundary condition (2.2*a*) at *y*=0, it is necessary to investigate the behaviour of the solution (2.14) near *y*=0. We definewhereWith this notation, direct differentiation of (2.14) gives the Taylor series as *y*→0. Higher-order terms can then be found either by further direct differentiation of (2.14) or more easily by expanding the coefficients of (2.9) as power series in *y*, givingas *y*→0. The condition *R*_{1}(0)=0 requires *C*_{3}=0, and equation (2.8*a*) then implies that(2.17)as *y*→0. Although there is no formal boundary condition on *ϕ*_{1} at *y*=0, we do require that the coefficient of 1/*y*^{2} be zero, otherwise the expansion (2.7*b*) is invalid for . Therefore which gives the required expression for *A*_{1}(2.18a)where(2.18b)(2.18c)This formula can be evaluated numerically using the software package Maple (Monagan *et al*. 2007); a Maple procedure for this is given in appendix B. In figure 4, we plot the dependence of *A*_{1} against both *a* and the original model parameter *ω*_{1}; recall that *a* is the wave amplitude generated by a linear obstacle, which depends on *ω*_{1} according to the formula derived by Sherratt (2003), and given in (2.3). A simple check of the formula (2.18*a*), (2.18*b*) and (2.18*c*) is given by the limiting value as *a*→1^{−}. Although we have been unable to evaluate (2.18*a*), (2.18*b*) and (2.18*c*) analytically in this limit, the limiting value of *A*_{1} can be deduced from (2.14), for which all the integrals can be done exactly when *a*=1. This implies the limiting solution formas *a*→1^{−}, so that *R*_{1}(∞)→0. Therefore as *a*→1^{−}, in agreement with the numerical results illustrated in figure 4.

## 3. Application to a predator–prey model

In §2, we calculated the leading order correction to the PTW form generated by a Dirichlet condition on a boundary that is circular, with a large radius. This calculation is for the particular case of a reaction–diffusion system of *λ*–*ω* form, and exploits the simplicity of the *λ*–*ω* kinetics when expressed in terms of amplitude and phase gradient. Consequently, an analogous calculation is not possible for more general reaction–diffusion systems. However, our results can be applied to any reaction–diffusion system with scalar diffusion close to a supercritical Hopf bifurcation in the kinetics, since the *λ*–*ω* equations are the normal form for such a system.

As a specific example of this, we consider a standard model for the interaction of a predator population and its prey. Reaction–diffusion models have been widely used to study spatio-temporal patterns in predator–prey systems (Gurney *et al*. 1998; Medvinsky *et al*. 2002; Sherratt *et al*. 2003; Morozov *et al*. 2006; Pearce *et al*. 2006; Garvie 2007). The model we consider has the dimensionless form(3.1a)(3.1b)This model is presented in many mathematical biology papers and textbooks (see, e.g. Murray 2002; Britton 2003; Turchin 2003). *P* and *Q* denote predator and prey densities, respectively. The dimensionless parameters *b* and *c* have simple ecological interpretations: *b* is the ratio of predator birth and death rates and *c* is the ratio of prey and predator birth rates. The parameter *d* reflects the rate at which prey consumption per predator saturates as prey density increases. There are three uniform steady states when *b*>(1+*d*)/*d*: *Q*=*P*=0 (trivial); *Q*=1, *P*=0 (prey only); and , (coexistence). The coexistence steady state is unstable when , and the kinetics then have a stable limit cycle.

Low values of the parameter *d* would typically correspond to a generalist predator, which has a number of alternative prey species in addition to *Q*, while high *d* would correspond to a specialist predator. Spatial variation in predation type is well documented in a number of systems. For example, in Fennoscandia, small rodents have both specialist mammalian predators (stoats and weasels) and a range of generalist predators (foxes, common buzzards, cats, etc.). The latter increase greatly in number from North to South, and this has been proposed as the cause of the well-known latitudinal gradient in population cycle period (Hanski *et al*. 1991). Below 60° N, there are no clear multi-year cycles in small rodent populations, and the hypothesis is that the parameter *d* (with *P* being stoats/weasels) decreases below the critical value (*b*+1)/(*b*−1) at this latitude. A more sudden change in *d* could occur when there is an abrupt change in habitat, such as the edge of a wood, or a boundary between moorland and farmland. Such a change could result in specialist predators dominating in the study domain (the wood or moorland), with primarily generalist predation in the surrounding environment. The model would then apply in both regions of space, but with a relatively large value of *d* in the study domain and a much smaller value in the surroundings. The latter would typically imply stability of the coexistence steady state, and it would then be natural to consider the behaviour only in the study domain, subject to the boundary condition *Q*=*Q*_{s}, *P*=*P*_{s}. This is directly analogous to the boundary condition *u*=*v*=0 used in §2 for the *λ*–*ω* equations. Again, we consider this boundary condition at the edge of a circular obstacle in an infinite domain; we denote the radial space coordinate by *ξ*, with the obstacle edge at *ξ*=*ξ*_{0}. Numerical simulations show PTWs moving away from the boundary at *ξ*=*ξ*_{0}, in a manner very similar to that illustrated in figure 1. As in the *λ*–*ω* case, the periodic wave amplitude and wavelength at large values of *ξ* depend on the boundary radius *ξ*_{0}, and this is illustrated in figure 5.

When *d* is close to the Hopf bifurcation value (*b*+1)/(*b*−1), we can use our results from the previous section to predict the variation in wave properties with *ξ*_{0}, when this is large. The first step in this calculation is to determine the appropriate values of *ω*_{0} and *ω*_{1}, via a reduction of (3.1*a*) and (3.1*b*) to normal form, close to the Hopf bifurcation point. This calculation is described in detail by Sherratt *et al*. (2003) and gives(3.2a)(3.2b)Substituting (3.2*b*) into (2.3) gives a formula for the wave amplitude *a* generated by a linear boundary in terms of the ecological parameters *b* and *c*. This can then be used to calculate the appropriate value of *A*_{1}. The reduction to normal form involves a rescaling of space (Sherratt *et al*. 2003), which implies that the appropriate radius of curvature to substitute into our previous calculation is(3.3)With these values of *a*, *A*_{1} and *R*, the first-order approximation to the wave amplitude is then given by *a*+*A*_{1}/*R*.

In practice, it is wavelength rather than wave amplitude that is usually reported in ecological studies of populations exhibiting PTWs (Bjørnstad *et al*. 1999; Liebhold *et al*. 2004). Again, we must take into account the rescaling of space that is required for the reduction to Hopf normal form, and this implies that our prediction of the wavelength at large values of *ξ*_{0} is given bywhere *R* is related to *ξ*_{0} via (3.3). Owing to the algebraic approach of the solution amplitude (and hence wavelength) to this limiting value, it is impractical to solve (3.1*a*) and (3.1*b*) numerically on a domain large enough for effective convergence. However, the far-field solution for the *λ*–*ω* equations, derived in §1, can be used to derive the appropriate correction. If *ξ*^{*} is the position at which the wavelength is being calculated, then (2.6) implies that the appropriate correction to the wave amplitude isThis in turn implies the approximation(3.4)for the wavelength at *ξ*=*ξ*^{*}. Here, *a* is given by substituting (3.2*b*) into (2.3), and *A*_{1} is then given by (2.18*a*), (2.18*b*) and (2.18*c*).

In figure 5, we compare this formula with numerical results from simulations of the model equations (3.1*a*) and (3.1*b*). The formula provides a good approximation to the numerical results provided that both *ξ*_{0} is large and *d* is close to the Hopf bifurcation point (*b*+1)/(*b*−1). In appendix C, we give estimates for the errors due to various aspects of the numerical scheme; these imply an error of at most 0.8% in our numerical calculation of the wavelength. Therefore, the formula (3.4) and the numerical results agree to within a numerical error for values of the obstacle radius *ξ*_{0} greater than approximately 900 for both values of *d*. For smaller obstacle radii, the main source of error in (3.4) comes from neglecting the *O*(*ϵ*^{2}) terms in the solutions of the *λ*–*ω* equations. Note that *ξ*_{0}=900 corresponds to *ϵ*≈0.02.

## 4. Discussion

Our objective in this paper has been to generalize previous results on the generation of PTWs by Dirichlet boundary conditions on linear boundaries. In the case of the *λ*–*ω* system of equations, we have derived the leading order effect on wave selection of the boundary being circular with a large radius. We have then shown how this result can be applied to a specific ecological model close to Hopf bifurcation in the kinetics.

One interesting consequence of the variation in wave amplitude with obstacle radius is its implications for wave stability. In the family of PTW solutions of an oscillatory reaction–diffusion equation, some members are stable as reaction–diffusion solutions, while others are unstable. Here and in the following, we are referring to ‘essential stability’ that determines stability on an infinite domain (Sandstede & Scheel 2000; Sandstede 2002; Rademacher *et al*. 2007). In general, there is no analytical expression for the division of the wave family into stable and unstable cases, but for *λ*–*ω* systems, Kopell & Howard (1973) showed that the PTW (1.2*a*) and (1.2*b*) is stable if and only if *r*^{*} satisfies(3.5)If the wave amplitude *a* generated by a linear obstacle is greater than *r*_{stab}, then an obstacle of any radius will generate a stable wave. However, if *a*<*r*_{stab}, a linear boundary will generate an unstable wave, while obstacles of sufficiently low radius will generate a stable wave. (Recall that as obstacle radius →0^{+}, the wave amplitude →1^{−}, which corresponds to the limit cycle of the kinetics, which is stable as a reaction–diffusion solution for scalar diffusion (Kopell & Howard 1973).) Numerical simulations indicate that when the obstacle radius is such that a stable wave is predicted, the whole solution (not just the behaviour far from the obstacle) is stable. However, when an unstable wave is predicted, the long-term behaviour is irregular spatio-temporal oscillations, with a band of PTWs visible close to the boundary in some cases. The natural interpretation of this is that the boundary selects a wave, which then destabilizes, giving spatio-temporal irregularities. In figure 6, we illustrate the transition from regular to irregular long-term dynamics as obstacle radius is increased.

Note that for this figure, we have deliberately used initial conditions without circular (or other) symmetry. When a stable wave is selected, these initial asymmetries rapidly disappear, but in the unstable case they persist in the irregular oscillations.

We have shown that Dirichlet boundary conditions, on a circular obstacle, generate PTWs in equations (1.1*a*) and (1.1*b*), with wave amplitude decreasing as a function of obstacle radius. A natural extension to this study is to consider obstacles with other geometries. This is a significantly more difficult problem, because it is fundamentally two dimensional in space. A limited exploration of this problem was performed by Sherratt *et al*. (2003), who simulated equations (3.1*a*) and (3.1*b*) in two dimensions, with a rectangular obstacle in the centre of a large domain. Their results suggested that wave selection depended primarily on the largest obstacle dimension. We performed a similar, but more systematic, study to that of Sherratt *et al*. (2003), for elliptical obstacles, with major and minor axes of half-length *L*_{major} and *L*_{minor}, respectively (figure 7).

As in the case of a circular obstacle, we found that the long-term behaviour far from the obstacle was an essentially linear PTW; thus the solution amplitude and phase gradient approach limiting values as the distance from the obstacle tends to infinity. Fixing *L*_{minor}=1 and increasing *L*_{major} from 1 to 100 resulted in a decrease in amplitude that mirrors the effects of varying the radius of a circular obstacle. However, when *L*_{major} is greater than approximately 10, the resulting PTW amplitude is less than for a circular obstacle of radius *L*_{major}. These results suggest that the sizes of both the major and minor obstacle dimensions play a role in determining the PTW properties. To illustrate this, we fixed *L*_{major}=20 and varied *L*_{minor} from 1 to 40 (figure 7). Amplitude initially increases with *L*_{minor} up to a maximum near where *L*_{major}=*L*_{minor} and then declines again. This indicates a significantly more complex dependence of amplitude on obstacle geometry than that suggested by the limited study of Sherratt *et al*. (2003), and a detailed analytical and numerical investigation of this is a natural area for future work.

Traditionally, ecological field studies have focused on temporal rather than spatial dynamics. However, the last decade has seen the publication of results from a number of long-term spatio-temporal field studies. Several of these studies demonstrate PTWs (Ranta & Kaitala 1997; Lambin *et al*. 1998; Moss *et al*. 2000; MacKinnon *et al*. 2001; Bjørnstad *et al*. 2002; Johnson *et al*. 2004; Bierman *et al*. 2006), suggesting that this may be a widespread spatio-temporal structure in ecological systems. A mathematical understanding of PTW generation is a crucial accompaniment to this fieldwork, in view of the enormous time and expense required for each field study. Our work represents one step in improving this understanding, showing that the curvature of landscape obstacles may affect the travelling waves generated by them in ecological systems. Moreover, our results in §2*a* predict that the wave properties will vary much more gradually with distance from the obstacle edge than in the case of a flat boundary. The PTWs reported in empirical studies to date have wavelengths that are relatively large compared with the size of the habitat, and it seems unlikely that any ecological domain is sufficiently long to allow the development of wave trains of many wavelengths (say more than 10). We therefore predict that if curved obstacles are generating PTWs in ecological systems then the wave properties will vary to a measurable extent with distance from the obstacle edge. It may be important to consider this when analysing field data. For example, it may explain why two studies of the PTWs in abundance of field voles in Kielder Forest (Northern UK; Lambin *et al*. 1998; MacKinnon *et al*. 2001), differing in the spatial and temporal scales of the data analysed, reported different wave properties (wave speeds of 19 km yr^{−1} reported by Lambin *et al*. (1998) compared with 14 km yr^{−1} reported by MacKinnon *et al*. (2001); these correspond to wavelengths of approx. 76 and 56 km, respectively).

Significantly, all of our results in this paper apply only close to Hopf bifurcation in the kinetics. An investigation of behaviour far from Hopf bifurcation is a natural topic for future study. With a number of spatio-temporal field studies due to report over the next few years, such a study would be particularly timely.

## Acknowledgments

This work was supported in part by an Advanced Research Fellowship from EPSRC (J.A.S., N.J.A.), a Doctoral Training Award from EPSRC (N.J.A.) and the NERC Environmental Mathematics and Statistics Programme (M.J.S.). J.A.S. thanks Adri Olde Daalhuis (University of Edinburgh) for help with the hypergeometric function.

## Footnotes

- Received August 31, 2007.
- Accepted November 2, 2007.

- © 2007 The Royal Society