The impact of short- and long-range perception on population movements

Journal of Theoretical Biology 460 (2019) 227–242

Contents lists available at ScienceDirect

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

The impact of short- and long-range perception on population movements S.T. Johnston a,b, K.J. Painter c,d,∗ a

Systems Biology Laboratory, School of Mathematics and Statistics, and Department of Biomedical Engineering, University of Melbourne, Parkville, Victoria 3010, Australia b ARC Centre of Excellence in Convergent Bio-Nano Science and Technology, Melbourne School of Engineering, University of Melbourne, Parkville, Victoria 3010, Australia c Department of Mathematics and Maxwell Institute for Mathematical Sciences, Heriot-Watt University, Edinburgh, UK d Dipartimento di Scienze Matematiche, Politecnico di Torino, Torino, Italy

a r t i c l e

i n f o

Article history: Received 7 August 2018 Revised 10 October 2018 Accepted 12 October 2018 Available online 16 October 2018 Keywords: Animal navigation Cell migration Nonlocal sampling Hilltopping Perceptual range

a b s t r a c t Navigation of cells and organisms is typically achieved by detecting and processing orienteering cues. Occasionally, a cue may be assessed over a much larger range than the individual’s body size, as in visual scanning for landmarks. In this paper we formulate models that account for orientation in response to short- or long-range cue evaluation. Starting from an underlying random walk movement model, where a generic cue is evaluated locally or nonlocally to determine a preferred direction, we state corresponding macroscopic partial differential equations to describe population movements. Under certain approximations, these models reduce to well-known local and nonlocal biological transport equations, including those of Keller–Segel type. We consider a case-study application: “hilltopping” in Lepidoptera and other insects, a phenomenon in which populations accumulate at summits to improve encounter/mating rates. Nonlocal responses are shown to efficiently filter out the natural noisiness (or roughness) of typical landscapes and allow the population to preferentially accumulate at a subset of hilltopping locations, in line with field studies. Moreover, according to the timescale of movement, optimal responses may occur for different perceptual ranges. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Navigation and migration play numerous crucial roles, for example allowing circulating immune cells to seek and destroy infections, and populations to collect at breeding grounds. A key ingredient for successful navigation lies in the availability of external orienteering cues, ranging from chemicals to visual landmarks, interpreted via (for organisms) a variety of sensory organs (eyes, ears, noses, lateral lines etc.) or (for cells) highly specific transmembrane receptors. The means by which cues are detected, processed and integrated are, naturally, subject to considerable speculation (Gould and Gould, 2012). Numerous models have been formulated to describe oriented movement, most frequently for chemical gradient responses (chemotaxis). Given that movement data is typically recorded at an individual level (e.g. cell/organism tracking), a natural approach

∗ Corresponding author at: Department of Mathematics and Maxwell Institute for Mathematical Sciences, Heriot-Watt University, Edinburgh, UK. E-mail address: [email protected] (K.J. Painter).

https://doi.org/10.1016/j.jtbi.2018.10.031 0022-5193/© 2018 Elsevier Ltd. All rights reserved.

is to model an individual’s movement as a velocity-jump random walk (VJRW) (Othmer et al., 1988). In a VJRW, an individual undergoes an alternating sequence of smooth “runs” with constant velocity and “turns” where a new velocity is chosen. We note that the latter provides a point for including orientation. In a groundbreaking paper, Patlak (1953) derived a governing partial differential equation (PDE) for motions of this form, a model of advection-diffusion type in which an external bias drifts the population in a specific direction. This model foreshadowed phenomenological models such as the well-known Keller and Segel (1970) model for chemotactic guidance, a popular choice for incorporating taxis-type orientation into models for biological movement (Painter, 2018). Many models implicitly assume that guidance information is local to the individual’s position. For an adult fruitfly, capable of detecting and orienting according to a (spatial) chemical gradient computed across its antennae (Duistermars et al., 2009), the length scale of antenna separation is microscopic in comparison to its movement in the landscape at large. The scale of gradient detection is, effectively, local. Yet in other cases the sampling region for cue detection may be much larger, covering a significant

228

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

portion of the individual’s environs. An obvious example lies in visual sensing, where individuals can assess information within eyesight. Here, the perceptual range – the “distance from which a particular landscape element can be perceived as such (or detected) by a given animal” Lima and Zollner (1996) – is a key concept and has been estimated at distances ranging from metres to kilometres, according to species, habitat and conditions (Flaherty et al., 2008; Fletcher et al., 2013; Olden et al., 2004; Pe’er et al., 2004; Prevedello et al., 2011; Todd et al., 2007; Zollner and Lima, 1997). Auditory cues are also detectable at a distance: a few kilometres for elephants (McComb et al., 2003) and possibly hundreds of kilometres for baleen whales (Payne and Webb, 1971), the latter possible via the properties of sound propagation in water. Less obviously, fish and amphibians detect the electric fields generated by nearby organisms and objects, a process termed electrolocation (Hopkins, 2005). Long-range assessments can also extend down to the single cell, with certain cells generating membrane protrusions (filopodia, cytonemes) that stretch multiple cell diameters to acquire inherently nonlocal information about their tissue landscape (Kornberg and Roy, 2014; Teddy and Kulesa, 2004). In this paper we explore models, founded on a simplified VJRW, that incorporate local and nonlocal evaluation of the environment. Orientation enters by biasing turns into certain directions, according to an underlying navigation field that can be adapted to account for various modes of sampling the environment. As prototype models we consider three simple responses to a scalar cue: a local gradient-based orientation, a nonlocal gradient-based orientation and a nonlocal response in which orientation is according to the maximum cue value over the perceptual range. The corresponding macroscopic models for the VJRW are of drift-anisotropic diffusion type. Via approximating expansions, we show how the simple local gradient model reduces to a Keller–Segel equation while the nonlocal models reduce to forms related to commonly-employed integro-PDE equations, such as those used to describe nonlocal chemotaxis responses (Hillen et al., 2007), cell interactions (Painter et al., 2015), swarming/flocking (Mogilner and EdelsteinKeshet, 1999) and, pertinently, information gathering over some perceptual range (Fagan et al., 2017). While these simplified forms strictly apply only when the orientation is weak, higher order approximations can extend the range of application. Separate approximations can be applied for very strong orientation responses, generating a pure drift equation in the simplest case. To demonstrate how nonlocal assessment contributes to population structuring, we consider a specific case-study application: hilltopping in butterflies, moths and other insects. In hilltopping, population members move up slope and localise at peaks and summits, a “lekking” behaviour that increases mating encounters (Alcock, 1987; Scott, 1968; Shields, 1967). Given that elevation is (primarily) sensed visually, we explore how perceptual range and assessment mode impacts on the macroscopic population movement. Gradually introducing noise into an idealised elevation profile, we show how nonlocal responses can allow a population to overcome roughness. Using terrain data for the Bodega Marine Reserve, a recent site of field studies (Grof-Tisza et al., 2017), these findings are shown to extend to natural environments, with nonlocal sampling allowing the population to accumulate on a select subset of prominent peaks. We discuss the results both generally and within the specific context of hilltopping.

work we focus on two-dimensional navigation (l = 2). This type of navigation is relevant for movement across a landscape. Individuals move via a VJRW with instantaneous waiting time between reorientations, requiring two probability distributions: a runtime distribution over R+ that dictates the reorientation rate and a turning distribution over velocity space V that describes the new velocity. The former is taken to be a standard Poisson process, that is, exponentially distributed runtimes with constant mean runtime τ (or turning rate λ = 1/τ ). We note that other choices of runtime can be made, potentially giving rise to subdiffusive or superdiffusive behaviour (Estrada-Rodriguez et al., 2018; Friedrich et al., 2006; Taylor-King et al., 2015). For the turning distribution we assume: (i) individuals move with a fixed speed s (i.e. V = sSl−1 ,; (ii) the new heading does not depend on the previous heading. Assumption (i) means that the turning distribution can be defined in terms of a directional distribution on the unit sphere, q(n|t, x), specifying the probability of choosing direction n ∈ Sl−1 following a turn, where n is the directional heading on the unit sphere Sl−1 . For the above VJRW, given stochastically independent walkers, we can obtain an evolution equation for the macroscopic population density at position x and time t, u(x, t). The process is to first write down the analogous transport (or kinetic) equation and, via scaling, obtain a macroscopic equation in an appropriate limit. We refer to Hillen and Painter (2013) (and references therein) for details and simply state that applying a moment closure method1 generates a macroscopic model in the form of the Drift-Anisotropic Diffusion (DAD) equation

u(t, x )t + ∇ · (a(t , x )u(t , x )) = ∇∇ : (D(t , x )u(t , x )).

(1)

In the above, D(t, x ) is a l × l diffusion tensor matrix and a(t, x ) is the l-dimensional advection velocity. The colon (: ) denotes the contraction of two tensors, resulting in a summation across all second order derivatives

∇∇ : Du =

l 

∂ ∂ Du . ∂ xi ∂ x j i, j=1

The macroscopic quantities a and D are determined from the statistical properties of the individual-level model: mean speed, mean runtime and the turning distribution. Specifically,

a(t, x ) = s

D(t, x ) =



τ

Sl−1

nq(n|t, x )dn,

(2)

(sn − a )(sn − a )T q(n|t, x )dn .

(3)

 Sl−1

Thus, the drift is proportional to the expectation of q and the anisotropic diffusion is proportional to its variance-covariance matrix. Environmental guidance is encoded in q, via biasing turns into specific directions. For the two-dimensional case considered here, the von Mises distribution (Mardia and Jupp, 20 0 0) is a natural choice:

q(n|k, ν ) =

1 e k n ·ν , 2π I0 (k )

(4)

where Ij (k) denotes the modified Bessel function of first kind of order j. Eq. (4) plays an analogous role to the normal distribution of linear statistics, and generates a predominance of turns into

2. Model framework We first outline the framework, referring to Othmer et al. (1988), Perthame (2006), Hillen and Painter (2013) and Wang and Potts (2017) for further details. We let t denote time, x ∈ Rl the spatial coordinate and v ∈ Rl the velocity. While the following derivation is general, in this

1 Our moment closure relies on two frequently used assumptions: (i) the population’s variance is computable from the equilibrium distribution, and (ii) a fast flux relaxation. The latter effectively states that, at the space/time scales of interest, the particle responds instantaneously to the information obtained at its current position. In later simulations we compare directly the distributions from repeated simulations of the stochastic model and the macroscopic model, confirming their validity for the present problem.

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

229

Fig. 1. (a-c) Illustration of (M1)-(M3) in the context of hilltopping behaviour. Red circles and blue stars indicate starting and expected final locations of the individual; green arrows represent the initial direction of the bias; red and blue lines indicate the sampling region at the start and at the end of movement. (a) Local slope response; (b) Individual responds to the elevation angle, with θ R > θ L implying a bias to the right; (c) Individual responds to the absolute elevation, with EL > ER implying a bias to the left. (d) Illustration of nonlocal sampling. An individual at x samples the environmental cue E at points within some compact region (here taken to be a circle centred on x and characterised by the perceptual range, R), averaging the response into a preferred directional bias (green arrow). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

preferred direction ν ∈ S1 with certainty increasing with concentration parameter k. We note that ν(x, t) and k(x, t) are functions of space and time to reflect variation in the direction and strength of the orienteering cue, but we omit the dependences for clarity of presentation. As k → 0, q tends to a uniform distribution, while as k → ∞, q approaches a singular distribution in which the preferred direction is always selected. Expectations and variance-covariance matrices of (4) can be explicitly calculated (Hillen et al., 2017), yielding the following driftvelocity and diffusion tensor

a(k, φ ) = s

D(k, φ ) =

I1 (k ) ν, I0 (k )

τ s2



2

(5)



I2 (k ) 1− I − τ s2 I0 (k )





I2 (k ) I1 (k )2 − A. I0 (k ) I0 (k )2

(6)

Intuitively, the drift is in the dominant direction. The anisotropic diffusion tensor has been decomposed into isotropic (proportional to the identity matrix I) and anisotropic (proportional to the singular matrix A = ννT ) components. The relative strength of diffusion to advection depends on k, via the modified Bessel functions. We note that I1 (k ) , I0 (k ) I1 (k ) • , I0 (k ) •

tion;

•

I1 (k )2 I0 (k )2

I2 (k ) I0 (k ) I2 (k ) I0 (k )

>

→ 0 as k → 0 and we converge to isotropic diffusion; → 1 as k → ∞ and we converge to a pure-drift equa-

I2 (k ) I0 (k )

for all k, and hence the dominating axis of diffu-

sion lies orthogonal to that of drift. Modified Bessel functions arise in numerous applications and a correspondingly large literature has emerged (Olver, 2010). In particular, approximating expansions are available for small and large arguments, as presented in Appendix A. 3. Orientation: Local and nonlocal sampling Orienteering information is encoded into a navigation vector field, denoted w(x, t ) ∈ Rl . The directions and vector lengths of w (x,t ) form the inputs for the von Mises distribution: ν (x, t ) = |w w (x,t )|

and k(x, t ) = |w(x, t )|. The navigation field in turn depends on a navigating factor, represented as a single scalar intensity function E(x, t), such as chemical concentration, elevation, temperature or magnetic intensity. We define w to have one of two generic forms: •

Local sampling – the navigation field is computed according to information obtained strictly at an individual’s current position;

•

Nonlocal sampling – the navigation field is computed over some spatially-extended sampling region, centred on the current position.

For both sampling mechanisms we consider a general form before tailoring into the following three model prototypes: (M1) a local and linear gradient model; (M2) a nonlocal and nonlinear gradient model; (M3) a nonlocal and nonlinear maximum model. The above have widespread applicability, but are particularly intuitive in the context of hilltopping. In (M1), responses are to the local gradient, as highlighted in Fig. 1(a). An individual, initially at the red circle, moves up the slope until coming to rest at the nearest local peak (blue star). In (M2), presented in Fig. 1(b), evaluations are made up to a given perceptual range (red/blue lines). The bias is weighted in the direction of the maximum elevation angle from the current position and, while all three peaks are observed from its initial position, the individual typically moves towards the closer of the two higher peaks: despite its lower overall height, it has the larger elevation angle from the initial location. In (M3), the individual estimates the true height of peaks and biases according to the highest peak within perceptual range, as highlighted by the schematic in Fig. 1(c). 3.1. Local sampling Local models implicitly assume pointwise sampling (at the scale of the overall environment). To generate a directional response we assume the local gradient of E can be determined. This results in a (local) tropotaxis response in which w points in the direction of ∇ E and the magnitude depends on E and |∇ E| (and other factors, if necessary). Thus,

w(x, t ) = k(E , |∇ E | )

∇E . |∇ E |

where k(E, |∇ E|) describes the strength and form of response: positive (negative) values of k imply attraction (repulsion). We note that k(E, |∇ E|) must satisfy k = 0 for |∇ E | = 0. Given data, k can be estimated through fitting, for example in leukocyte chemotaxis where chemokine responses have been calculated for different concentrations and gradients (Jeon et al., 2002). To generate the prototype model (M1) we invoke the simple choice of linear gradient dependence: k = κ1 |∇ E | with associated response coefficient κ 1 . The navigation field is therefore directly proportional to ∇ E, yielding

w = κ1 ∇ E .

(M1 )

230

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

We remark that similar models have been proposed in a cell biology context, where E arises from local crowding between cells (Binny et al., 2016). Movement in these models can either be repulsive and hence proportional to −∇ E, or attractive and hence proportional to ∇ E (Binny et al., 2016). We note that such models consider navigation based from interactions between individuals, whereas we focus on processes where individuals do not interact with each other. Instead, navigation information is obtained from an underlying environmental cue. 3.2. Nonlocal sampling We consider the schematic in Fig. 1(d) and suppose an individual at x computes a movement response by directly sampling the environment at positions x + ri for i = 1 . . . L, where L is total number of perceived locations. Sampling could occur through (for an organism) visual focus on a distant point or (for a cell) extending a filopodium. We define a probability distribution, (r|x, t), that denotes the probability that position x + r is sampled from position x at the time of reorientation t. Logically, has compact support to reflect a maximum perceptual range, but analytically convenient forms with decaying tails (such as Gaussians) may also be reasonable. We suppose the strength of attraction towards (or repulsion from) x + r is encoded into a response function g(E (x + r, t ), E (x, t ), |r| ), potentially depending on the cue at both the current and sampled points, as well as the sampling distance. The attraction vector pr (x, t) in direction r/|r| due to information at x + r is hence

pr (x, t ) = (r )g(E (x + r, t ), E (x, t ), |r| )

r

|r|

L 1 r (r|x, t )g(E (x + r, t ), E (x, t ), |r| ) . L |r|

If a large number of points are sampled we can approximate w via the integral form

w ( x, t ) =

Rl

(r|x, t )g(E (x + r, t ), E (x, t ), |r| )

r

|r|

dr .

(7)

This generic formulation is easily adapted and we next consider the simple prototype models. Note that here, for simplicity, we exclusively concentrate on uniform sampling over a circular region of radius R, where we define R to be the perceptual range

(r|x, t ) =

1

VR

0

E ( x + r, t )n − E ( x, t )n . |r|

The parameter n gives additional control, increasing the weighting of attraction bias as n is increased. This gives the following navigation field (M2)

w(x, t ) = α2



E ( x + r, t )n r dr . |r| |r| BR ( 0 )

(M2 )

We note that we have absorbed various constants into α2 =

κ2 ω 1 , where VR ωn

    E ( x + r, t )i r  ωi =  dr. r r | | | | BR ( 0 ) The scaling normalises the navigation field, so that the magnitude remains essentially the same as n is altered. Note that Taylor expansion again yields (M1) in a small R approximation. 3.2.2. Maximum value evaluation Our second prototype nonlocal model takes g to depend on the absolute size of the cue at E (x + r, t )

if r ∈ BR (0 ) , otherwise

where κ 3 is the sensitivity coefficient. Substituting into (7) along with (8) now yields

w ( x, t ) =

κ3



VR

E ( x + r, t )

BR ( 0 )

r dr . |r|

The distinction with the gradient-based evaluation is exemplified through a Taylor expansion, where we find

i=1



g(E (x + r, t ), E (x, t ), |r| ) ∝

g(E (x + r, t ), E (x, t ), |r| ) = κ3 E (x + r, t ),

.

We suppose w is formed from the net attraction, obtained by sampling over different locations and calculating the average. Thus, we take

w ( x, t ) =

Hence, as R → 0 the nonlocal gradient model converges to the local gradient model with κ2 = 2κ1 . It is straightforward to incorporate nonlinear processing of the cue, and in particular we consider

(8)

where BR (0) is the ball of radius R, centred on the origin, and VR is its volume (i.e. VR = π R2 for l = 2).

w(x, t ) = κ3 R∇ E + O(R3 ) . Hence, as R → 0 the navigation vector field shrinks to zero length and orienteering information is lost. As before, we can extend to a nonlinear form to obtain (M3)

w(x, t ) = α3 where α3 =

 BR ( 0 )

κ3 ω 1

VR ωn

E ( x + r, t )n

r

|r|

dr,

(M3 )

and

    i r  ωi =  E ( x + r, t ) dr, r || BR ( 0 ) provides the normalisation as n is altered.

3.2.1. Nonlocal gradient sampling Suppose attraction is according to the nonlocal gradient of E

g(E (x + r, t ), E (x, t ), |r| ) = κ2

E ( x + r, t ) − E ( x, t ) , |r|

with sensitivity coefficient κ 2 . Substituting both the above and (8) into (7), we obtain

w ( x, t ) =

κ2 VR



BR ( 0 )

E ( x + r, t ) r dr . |r| |r|

Notably, a Taylor series expansion about x applied to the RHS yields

w ( x, t ) =

κ2 2

∇ E + O ( R2 ) .

4. Approximations A choice must be made about whether to perform direct stochastic simulations of the VJRW model or solve the full DAD model (1) with their Bessel function dependent coefficients (5) and (6). Both approaches have associated costs. Simulations of the stochastic VJRW are expensive, scaling with population size; evaluations of the modified Bessel functions and implementing anisotropic diffusion terms are also costly. In an analytic context, anisotropic diffusion terms may be more difficult to deal with than isotropic counterparts. Consequently, it is worthwhile investigating simplifying reductions.

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

231

Fig. 2. Four different idealised elevation profiles.

4.1. Local sampling approximations We consider model (M1). Expansions of the modified Bessel functions (see Appendix A) allow a series of approximations to (1) to be obtained, relevant for small or large values of k (weak or strong navigating cues, respectively). For weak cues the lowest order approximation (O(k1 )) yields the well known Keller–Segel equation for (chemo)taxis

ut = d∇ 2 u − ∇ · (χ u∇ E ),

(9)

τ s2 /2)

with constant diffusion (d = and chemotactic (χ = sκ1 /2) coefficients2 . The above is strictly valid only for weak cues, but a series of higher order “corrections” extends the range of k values for which the approximation accurately describes the full DAD model. We tabulate up to O(k4 ) in Appendix B. For example, the O(k2 ) approximation in two dimensions introduces a correction for each of the diffusion and chemotaxis terms, such that



ut = D∇ ∇ :



1−



κ12 |∇ E |2 κ12 |∇ E |2 I − 2 2 2 8 + 2κ12 |∇ E | 4 + 2κ12 |∇ E |  

κ12 |∇ E |2 A u 2 8 + 2κ12 |∇ E |     2 2 + κ12 |∇ E | −∇ · χ 1 − u∇ E . 2 4 + κ12 |∇ E | −

(10)

Notably, this correction introduces the anisotropic diffusion component. We remark that diffusion coefficients remain positive for all κ 1 , so the system is well-behaved. As k becomes larger the approximations become less useful: more terms are required to achieve reasonable accuracy, and it would be natural to instead revert to the full DAD model (1). However, very large k can alternatively be dealt with through large argument approximations (Appendix B). For example, for sufficiently large k (corresponding to strong navigation fields) we can approximate the full DAD model (1) via the pure drift-equation



ut + ∇ ·

s u∇ E |∇ E |



= 0,

which simply generates direct movement with speed s in the direction of the steepest gradient. Corrections to this again yield a sequence of chemotaxis/anisotropic diffusion style models, tabulated in Appendix B. For example, the O(k−1 ) approximation in 2D generates

   1 τ s2 s 2κ1 |∇ E | − 1 ut = ∇∇ : ( I − A )u − ∇ · u∇ E . κ1 |∇ E | |∇ E | 2κ1 |∇ E | 2 We remark the general choice k would instead yield a general form Keller–Segel equation in the first order approximation:



ut = d∇ 2 u − ∇ · k(E , ∇ E , |∇ E | )u∇

E

|∇ E |



.

The anisotropic diffusion in the above corresponds to a degenerate/singular form in which diffusion is one-dimensional and along the axis orthogonal to the preferred direction. 4.2. Nonlocal sampling approximations The process is identical for nonlocal sampling scenarios, and we therefore only state the first order weak cue form for the general navigation field (7), and for both (M2) and (M3). In the general case, the first order approximation under a weak cue yields a nonlocal integro-PDE equation

   r ut =d∇ u−∇ · σ u (r|x, t )g(E (x + r, t ), E (x, t ), |r| ) dr , |r| Rl 2

(11)

τ s2 /2

for d = and σ = s/2. The above is closely related to various nonlocal PDE models. For example, if E is directly based on the population’s own density distribution the above is similar to those employed to describe swarming (Mogilner and EdelsteinKeshet, 1999) and cell interactions (Painter et al., 2015). Moreover, the above is closely related to the model proposed in Fagan et al. (2017) for perceived gradient following. If we consider the two prototype models, the above becomes

ut = d∇ 2 u − ∇ · for (M2), and

ut = d∇ 2 u − ∇ ·

  σ α2 u   σ α3 u



E n ( x + r, t ) r dr , |r| |r| BR ( 0 )



BR ( 0 )

E n ( x + r, t )

r dr , |r|

for (M3). 4.3. Numerical validation We restrict our validation to local gradient sampling under an environmental cue that does not vary with time. The population is initially distributed uniformly across the rectangular region 0 ≤ x ≤ 1, 0 ≤ y ≤ 2. We generate a quasi-one-dimensional scenario by assuming the cue varies only with y, specifically

E (y ) = e−β (y−1) . 2

(12)

The corresponding cue profile is presented in Fig. 2(a). The dominant direction for individuals with locations y < 1 (y > 1) will be π /2 (−π /2) and subsequently individuals are expected accumulate along the ridgeline y = 1. Comparisons are presented for (i) the average behaviour of multiple realisations of the velocity-jump process, (ii) numerical solution of the full DAD model (1) under (5) and (6), and (iii) the O(k) to O(k4 ) approximations. These approximations rely on an assumption that the strength of the navigating cue, k, is small. Numerical methods are described in the Supplementary Materials, Section S3. Briefly, an implicit timestep method and an operator splitting method are used to numerically solve the full DAD model and the corresponding O(k) to O(k4 ) approximations. The trajectories of individuals in the velocity-jump

232

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

Fig. 3. (a),(d),(g) Comparison between the drift-velocity function (red), and the corresponding O(k) to O(k4 ) approximations (orange, purple, green, blue, respectively) for (a) κ1 = 0.5, (d) κ1 = 1, (g) κ1 = 2. (b),(e),(h) Comparison between the diffusivity function (red), and the corresponding O(k) to O(k4 ) approximations (orange, purple, green, blue, respectively) for (b) κ1 = 0.5, (e) κ1 = 1, (h) κ1 = 2. (c),(f),(i) Comparison of the average velocity-jump behaviour (black, dashed), and the numerical solutions to the full DAD model (red) and the O(k) to O(k4 ) approximations (orange, purple, green, blue, respectively) for (c) κ1 = 0.5, (f) κ1 = 1, (i) κ1 = 2. Random walk parameters are set at s = 0.001 and τ = 1; elevation is set by (12) with β = 10; simulations are based on averaging over 106 realisations of the VJRW and a spatial discretisation of y = 0.001 for the numerical approximation of the continuous models. Simulations are computed until tend = 100. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

process are simulated directly as described in Painter (2014), and additionally detailed in the Supplementary Material, Section S3. The Bessel function approximations are independent of s and we therefore expect that the difference between approximations and the full DAD model (1) is constant as s changes. Thus, differences between solutions under increasing s should stem only from the assumptions made when deriving the macroscopic model from the velocity-jump process. The difference between the approximations and either the full DAD model solution or the average behaviour of the velocity-jump process decreases with the order of the approximation. Further, both the full DAD model solution and approximations prove to be less accurate against the average velocity-jump behaviour as s increases (Supplementary Material, Section S1), consistent with previous investigations (Painter, 2014). Generally, whether (1) provides a reasonable description of the underlying random walk depends on whether the spatial and temporal scales are appropriately macroscopic, for example see Codling et al. (2008) and Hillen and Painter (2013) for further discussion. Specifically, average run lengths need to be short, compared to the overall domain dimensions, and individuals need to make numerous turns over the study time. We next explore the match for fixed s = 0.001 and varying response coefficient κ 1 . Varying κ 1 directly influences the concentration parameter k, so increasing κ 1 should impact on the matching ability of the approximations. Note that for a

quasi-one-dimensional case the diffusion tensor reduces to a function depending on y and k, denoted d(y), while the drift reduces to a “chemotaxis”-type function, denoted a(y). As expected, truncating at higher order terms yields a better fit for these functions, see Fig. 3(a,d,g) for a(y) and Fig. 3(b,e,h) for d(y). For sufficiently small κ 1 (e.g. κ1 = 0.5) even the lowest order approximations prove highly acceptable, while larger values (e.g. κ1 = 2) demand a higher order for a good fit. Turning to full numerical solutions, we compare O(k) to O(k4 ) approximations with those of the full DAD model and the average velocity-jump behaviour. Note that for the parameter values under consideration, the full DAD model closely describes the average velocity-jump behaviour. For κ1 = 0.5, Fig. 3(c), all approximations are generally successful, although O(k) and O(k2 ) approximations show some inaccuracy in regions where the diffusivity function approximation is at its worst. Increasing κ 1 generates a poorer fit: under this cue distribution and κ1 = 2, the maximum value of k ≈ 5.43 is beyond the expected range of k values for reasonably valid approximations. Nevertheless, the O(k4 ) still proves reasonably accurate. Overall, numerical comparisons validate the approximations, but only when k remains within relevant regions. Thus, if it is known that a particular cue provides only a weak bias (say, a cell migrating in a fairly stable and shallow gradient) it is reasonable to use the simplified PDE, such as the local or nonlocal

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

233

Keller–Segel equation. If, rather, cue strength cue varies widely, one should resort to employing the full DAD Eq. (1), assuming the spatial/temporal scales are appropriately macroscopic.

results and further discussion are provided in the Supplementary Material, Section S2.1. Thus, all proposed mechanisms are plausible and we shift to a more sophisticated analysis of the mechanisms.

5. Short and long-range sampling: Hilltopping case study

5.1.1. Perceptual range and nonlocal response order Here we investigate the influence of the perceptual range, R, and nonlocal response order, n. We consider a landscape containing two peaks of different heights, as highlighted in Fig. 2(b). For model (M3) we observe a significant change between R = 1.5 and R = 2.0. For smaller radii the population is split into separated subpopulations, one on each hilltop with the majority of the individuals on the higher of the two peaks. In contrast, for R = 2.0 the individuals merge into a single contiguous population at the higher peak. Comparing the dominant direction for the two radii, Fig. 4(c) and (f), a transition is observed such that under smaller radii each peak has its own “zone of attraction”. For larger radii, however, the zone of attraction for the larger peak expands to cover the entire domain. Note that for this example, n does not qualitatively change the results. Under model (M2) with a high nonlocal order (n = 5), two distinct sub-populations are maintained even for large R. That is, even when the higher peak is fully within perceptual range of members on the lower peak, the two sub-populations remain separate. Here, higher weighting is applied according to the largest elevation angles. Consequently, the lower peak can always be perceived as more attractive for individuals located sufficiently close. As the nonlinear order is decreased, however, the sub-population that initially formed on the lower peak begins to migrate towards the higher peak. This is reflected in the comparison of the dominant direction, Fig. 4(l), where we observe a collapse from two zones of attraction to one. As the nonlocal gradient model effectively places a weight on higher environmental cue values close to an individual, it is difficult for an individual to move away from the peak if that weighting is further compounded by a high nonlinear order. To examine whether migration is unidirectional (from the lower to higher peak), we consider the migration proportion for two initially separated populations: one uniformly distributed around the lower peak and the second about the higher peak (see Supplementary Material, Section S2.1). Notably, when peak to peak migration occurs it is only for those initially located on the lower peak, which translocate to the higher one. We therefore investigate how the combination of perceptual range, ratio of peak heights and nonlocal order influences the proportion of the lower peak population that have migrated by a fixed time, set at t = 50 0 0. The influence of peak height ratio is determined via two elevation profiles, Fig. 2(b) and (c), with peaks located at the same position but with distinct heights to generate ratios of 1.5 and 2.0 respectively. Fig. 5(a) summarises the results for (M3). For all scenarios, there is negligible migration for R = 1 and full migration by R = 2. In between, there is a steep transition with its location changing according to nonlocal order and peak height ratio: increases in either will lower the perceptual range required to shift the population. Effectively, an increase in either of these properties increase the attractiveness of the higher peak and the population migrates accordingly. Fig. 5(b) summarises the results for (M2). Notably, there is now a significant difference between n = 1 and n = 5 for both profiles. As discussed previously, the extra weighting toward higher environmental values close to an individual, introduced by the 1/|r| factor in the nonlocal gradient term, influences the dominant direction (Fig. 4(o)). Substantially altering the disparity between peak heights, though, can overcome this additional weighting.

The accumulation of butterflies, moths and other insects at hilltops has been widely documented (Alcock, 1987; Ehrlich and Wheye, 1986; Grof-Tisza et al., 2017; Pe’er et al., 2004; Scott, 1968; Shields, 1967). Many insect populations fluctuate over different locations and seasons, and strategies that offer robustness against extinction events are necessary. Hilltopping is perceived as one such mechanism, due to an improved rate of mating encounters. The relatively local scale over which movements occur, the precision to which terrain can be documented and the (relative) ease of tracking/counting population members combine to form an ideal system for determining how environments structure a population. Given visual assessment of topography, it is natural to query how perceptual range impacts on hilltopping. Natural terrains are typically “rough”, so responses based strictly on the local gradient could easily trap a large percentage of the population, scattered through small accumulations on many local peaks. Instead, studies of a hilltopping butterfly (Melitaea trivia) indicate that their orientation may be according to both the local slope and higher peaks within a longer range ( ∼ 50 metres) (Pe’er et al., 2004). Moreover, a recent capture-mark-recapture study of the tiger moth (Arctia virginalis) revealed that populations accumulated only on a subset of summits, with the largest numbers found at the highest elevations (Grof-Tisza et al., 2017). Motivated by these findings, we use this section to explore how short- and long-range topographical sensing alters hilltopping outcomes. We consider our three prototype models for moving up slopes, effectively describing populations that respond to: (M1) the strictly local gradient; (M2) the nonlocal gradient, that is, according to the elevation angle of terrain within perceptual range; (M3) the true height of land within perceptual range. We first generate idealised terrains in order to test principles before considering an application using elevation data for the Bodega Marine Reserve, as used in a recent hilltopping survey conducted by GrofTisza et al. (2017). We note that several previous simulation studies have explored hilltopping. For example, in Pe’er et al. (20 06, 20 05) an agentbased modelling approach was used to simulate butterfly movement paths, based on field study data (Pe’er et al., 2004). Topography was shown to channel movement paths along “virtual corridors”. Further agent-based studies in Pe’er et al. (2013) highlighted the importance of an element of randomness in the movement paths to avoid the above described trapping. In Painter (2014) the multiscale framework used here was employed to demonstrated how accumulations formed by hilltopping could optimise mating for low-density populations, but may be disadvantageous for more abundant populations. The current study here extends that work, specifically incorporating long-range perception to explore its impact on population distribution. 5.1. Idealised study We first explore how the prototype orienteering mechanisms impact hilltopping in a series of idealised studies. We assume a quasi one-dimensional terrain and, unless stated otherwise, initially suppose the population to be uniformly randomly distributed across the terrain. We first examine the suitability of models (M1)–(M3) to describe hilltopping, considering movement on a terrain containing a single peak centred at y = 1, as presented in Fig. 2(a). For succinctness we simply note that all model prototypes allow the population to relocate to the peak by t = 100;

5.1.2. Influence of noisy environments So far we have considered smoothly varying cues. In reality, environments are expected to be considerably rougher, due to small-scale variations in topography. Since a strictly local response

234

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

Fig. 4. Comparison of different R and n for the full DAD model (solid) and the average velocity-jump process (dashed), respectively, using (a)-(f) (M3) and (g)-(l) (M2). Cyan/blue corresponds to n = 1, orange/black corresponds to n = 5. Initially the population is uniformly distributed and the elevation profile corresponds to Profile Two in Fig. 2. Parameters used are as follows: (a)-(f) κ3 = 3, (g)-(l) κ2 = 0.3; random walk parameters are set at τ = 1, s = 0.01; numerical simulations based on averaging over 2 × 105 realisations of the VJRW and using a spatial discretisation of y = 0.001 for the numerical approximation of the continuous models. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. Proportion of population initially located on the lower peak that migrates to the higher peak by t = 50 0 0 for various R, n and peak height ratios. (a) Results for (M3), using κ3 = 3; (b) Results for (M2), using κ2 = 0.2. Other parameters are set at y = 0.003, τ = 1, s = 0.01.

depends on the gradient at a specific point, this variation may significantly impact on an individual’s ability to move across some environment. To explore this, we introduce a form of gradient noise into our environmental cue: we refer to Appendix C for details, but remark that the added noise contains local spatial correlation to ensure a differentiable underlying environment. We examine the average displacement of a population moving in response to a linear topographical profile, Fig. 2(d), subject to the addition of various levels of noise. Initially, the population is uniformly distributed between y = 0 and y = 1.2 (and zero elsewhere). Average displacement is calculated by numerically solving the full DAD Eq. (1) under twenty identically-prepared noise terms. That is, we generate twenty different noise terms accord-

ing to the same process, with the only difference stemming from the numbers sampled from the (pseudo)-random number generator. We then solve (1) for each of the noise terms and calculate the overall displacement of the population, defined as the difference between the median of the population at t = 0 and t = 500. We repeat this process for different levels of noise and scale according to the average displacement under zero noise. The results for models (M1) and (M3) are summarised in Fig. 6 (a), where for the nonlocal models we consider R = 0.2, R = 0.5 and R = 1.0. As might be expected, increasing the strength of the noise compared to the strength of the environmental cue results in a decrease in the average displacement. Local sensing (M1) is the most susceptible to the added noise, with increasing noise signifi-

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

235

Fig. 6. (a) Scaled average displacement by t = 500 for a population of individuals initially uniformly distributed between y = 0 and y = 1.2 for both local (solid) and nonlocal (dashed) responses and three R values. For each parameter combination the average displacement is obtained from 20 simulations on identically-generated terrain. Parameters used are κ1 = 0.3, κ2 = 0.6, κ3 = 2.0, y = 0.003, τ = 1, s = 0.01, n = 1. (b) Average proportion of the population initially located on the lower peak in Profile Two that migrates to the higher peak by t = 50 0 0 for M3 and R = 1.75 (cyan), M3 and R = 2.0 (orange), M2 and R = 2.2 (green), and M2 and R = 2.8 (purple). For each parameter combination the average migratory proportion is obtained from twenty simulations on identically-generated terrain. Parameters used are κ2 = 0.6, κ3 = 1.5, y = 0.003, τ = 1, s = 0.01, n = 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

cantly reducing the displacement. Nonlocal models (M2) and (M3) prove considerably more resilient, particularly for large perceptual ranges. This result is intuitive, as the nonlocal response involves averaging the environmental cue over the perceptual region, reducing the contribution of any noise. Note that as R → 0 for (M2), we observe convergence to (M1), consistent with our earlier findings. Overall, given that environments in nature are non-smooth and exhibit variation, nonlocal responses may be necessary for successful migration. We now return to the example with two peaks, with two initially distinct populations, and examine how the presence of noise influences the successful migration from the lower to higher peak. We note that a small amount of noise can elicit a large change in the proportion of the population that undergoes migration (Supplementary Material, Section S2.1). To determine how different levels of variation influence hilltopping, we append twenty identically-prepared noise terms to Profile Two. Using two distinct R values we examine the proportion of individuals initially located on the lower peak that migrate to the higher peak by t = 50 0 0, presenting the results in Fig. 6 (b) for (M2) and (M3). As expected, we find that increasing the amount of noise generally reduces the ability for individuals to translocate to the higher peak for both models: for combinations of R and κ that allow the entire population to relocate in the absence of noise, introducing noise inhibits a portion of the population from doing so. Consistent with previous findings, larger perceptual ranges offset noise-induced trapping. As naturally occurring environments do not typically exhibit the smoothness of idealised models, this suggests that the mechanism by which sensory information is processed must be able to account for the noise. 5.2. Bodega Marine Reserve study 5.2.1. Study site and elevation data We extend our case study to a natural environment, utilising topographical data for the Bodega Marine Reserve, a site of recent hilltopping studies (Grof-Tisza et al., 2017). High spatial resolution LIDAR data was obtained from the United State Geological Sur-

vey’s “The National Map” project over the co-ordinate range of the study site, and longitude, latitude, elevation and classification (water, ground etc.) are extracted at each data point. Interpolation onto a regular grid is performed, and the resulting topographical structure is mapped in Fig. 7(a). This topographical structure allows individuals to obtain navigation information in a manner that mimics the behaviour of a population undergoing directed motion according to the topography in the Bodega Marine Reserve. A population of butterflies/moths is uniformly randomly distributed across the landmass at t = 0 and their movements are followed over 20 h. Of course, more realistic distributions should account for oviposition/larval sites but a uniform distribution limits any subsequent bias in the final distributions. Simulation methods are described in the Supplementary Methods. Note that we impose boundary conditions that effectively confine the population within the bounds of the region plotted in Fig. 7(a), as described in the Supplementary Methods, Section S3. For later analyses we subdivide the land region using two methods. The first is a peak-centred scheme, Fig. 7(b), where 200 m by 200 m (4 hectare) non-overlapping regions are defined according to topography: (i) we find the complete set of local maxima; (ii) zone 1 is centred on the highest maximum and both this maximum and any others within distance are excluded; (iii) zone 2 is centred on the next highest non-excluded maximum and so forth. The second subdivision, Fig. 7(c), is a regular tessellation that intersects with the land coverage. The use of two methods is to reduce bias due to zoning. We consider the three prototype models for topographic sensing, (M1)–(M3). Note that in order to concentrate on how altering the navigation field and its associated parameters impacts on population structuring, we fix the speed, turning rate and sensitivity coefficients at values similar to those in an earlier study (Painter, 2014). Specifically, we set s = 4 m/min and a mean run time of τ = 1 min; note that these values account for pauses/resting periods during flight. Sensitivity coefficients are set at κ1 = 5, κ2 = 10 and κ3 = 0.2; the relationship between κ 1 and κ 2 ensures convergence from (M2) to (M1) as R → 0; the relationship between κ 2 and κ 3 is to generate

236

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

Fig. 7. (a) Topography of Bodega Marine Reserve study site. (b)-(c) Zones used for population counts, following either (b) peak-centred or (c) tessellated gridding.

Fig. 8. Representation of orientation information, for w given by (M2) with (a) R = 0, (b) R = 25, (c) R = 100 and (d) R = 500. Note that the R = 0 case is effectively the local gradient model (M1). For these plots we use n = 5 as the nonlinear coefficient.

quasi-equivalent navigation fields for (M2) and (M3) when R = 50 m. To provide context to these values, a butterfly located on a linear ramp of gradient 10% would choose an uphill direction on just over 75% of occasions, rising to slightly more than 90% for a 20% gradient.

5.3. Navigation field variation We first consider how orientation information changes as we shift from local to nonlocal models. For ease of illustration we plot only a portion of the full study site (the top left region in Fig. 7(a)). Plots are presented in the form of a “direction field”, where red arrows indicate the local direction of w at selected points and blue lines indicate the trajectory traced out if an individual exactly follows these arrows. Light blue circles indicate trajectory end points, and can be interpreted as attractors that potentially lead to local accumulations of the population.

The plots in Fig. 8 show results for (M1), i.e. R = 0, and (M2) with R = 25, R = 100 and R = 500; note that these plots assume n = 5. Under local sampling we observe significant fluctuations in arrow directions, with some arrows pointing in the direction of prominent peaks but others directed elsewhere. Numerous local maxima exist, so that trajectories travel short distances before becoming trapped at an attractor. Transitioning to the nonlocal model, however, dramatically reduces the number of attractors. The majority of trajectories now converge on a select few summits (examples given by black arrowheads) and trajectories can be significantly longer than the perceptual range. Effectively, movement is first according to the highest elevation angle within range, but subsequent higher peaks may be observed. For example, in Fig. 8 (c) we see trajectories converging on the top-left attractor over a distance of approximately 1 km, an order of magnitude greater than the perceptual range. Increasing R steadily reduces the number of attractors, such that by R = 500 all trajectories converge on one of just two attractors within the plotted region,

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

237

Fig. 9. Simulations of the VJRW and full DAD model (1) using elevation data for the Bodega Marine Reserve. Navigation field given by (left to right) the local model (M1), i.e. R = 0, and (M2) with n = 5 and R = 25, R = 100 and R = 500. (a)-(d) Simulations of the VJRW. Blue circles represent end locations and red dots indicate trajectories, for 250 individuals that are initially randomly uniformly distributed across the landmass. (e)-(p). Corresponding population distribution predicted by (1) at (e)-(h) 2, (i)-(l) 8 and (m)-(p) 20 h. Density plotted according to the colourbar in (e), where M denotes the density if the population is uniformly distributed across the landmass. Parameter values are as described in the main text. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

corresponding to prominent peaks. These peaks are not necessarily the two highest peaks across the studied zone, rather they correspond to the highest within their region. Further analyses have been performed for (M2) with n = 1 and (M3) with both n = 1 and n = 5 (see Supplementary Materials, Section S2.2). Qualitatively we observe similar phenomena, with a reduction in the number of attractors as R is increased. However, a few subtleties emerge: (i) higher nonlocal orders promote the persistence of attractors at lower prominences, consistent with our earlier idealised study; (ii) for n = 1, attractors can occasionally form some distance away from the nearest local peak. This latter observation stems from the fact that, under a lower nonlocal order, the weighting to summit points is diluted and evaluation is more according to a broad assessment of the terrain. However, broadly the results are consistent and we note that perceptual range, rather than the precise sensing model or degree of nonlinearity, appears to be the primary determinant in the number of attractors. For the remainder of our investigation we therefore restrict to presenting the results for (M2) under n = 5, noting that results for the other cases are qualitatively consistent.

5.4. Comparison of the velocity jump and macroscopic models We proceed to simulate the VJRW and the full (unapproximated) macroscopic model (1) across a range of R, using (M1), i.e. R = 0, and (M2) with R = 25, R = 100 and R = 500. Results are presented in Fig. 9. Agent-based model simulations of the VJRW are presented in Fig. 9(a)–(d), where blue circles denote final locations (after 20 h) and red dots indicate trajectories. The results parallel the previous direction field plots. For the local model, some larger collections occur at the most prominent peaks but the overall population remains scattered with numerous small accumulations trapped at local maxima. Expanding the perceptual range reduces trapping, allowing the population to coalesce into fewer and larger aggregates that correlate with the prominent peaks. The capacity for long-range information to influence movement is particularly illustrated by the individuals that move across the bays, attracted to peaks on the other side (arrowheads in Fig. 9(d)). Under each agent-based simulation we show the corresponding population distributions predicted by the macroscopic model

238

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

Fig. 10. Comparison between the full DAD model (1), (a) and (c), and the first order approximation obtained under the assumption of a weak cue, (b) and (d): (a,b) t = 4 h; (c,d) t = 20 h. Simulation details as in Fig. 9, where for this comparison we use (M2) with n = 5, κ2 = 10 and R = 100.

Fig. 11. Results of zoning analysis for (M1) and (M2) under various R. (a) At the end of the simulation (t = 20 h) we calculate the number of zones for which the population has increased to × 2, 5, 10, 20 the mean density. (b) We plot the normalised encounter rate as a function of time for (M1) and (M2) under various R. All results use the peak-centred zoning method while model parameters are as described in the text.

after 2, 8 and 20 h. Note that the distributions correlate well with those of the agent-based model, suggesting that the macroscopic model accurately captures individual-level behaviour at the spatial and temporal scales of study. We subsequently focus on the macroscopic model, exploiting its computational advantages. The simulations confirm the results of the agent-based model, with large sampling radii allowing the population to coalesce into fewer and larger aggregates: for example, for R = 500 we observe the population has coalesced into just 3 principle populations. We also test the validity of using an approximation, comparing the full DAD model with its first order approximation (assuming a weak cue) in Fig. 10. Subtle difference emerge in regions corresponding to the strongest field strengths but, overall, the first order approximation performs remarkably well: differences between the approximation and full DAD model are likely to be small in comparison to errors arising from parameter estimation. Moreover, differences are primarily noticeable at early times, with the results in the quasi-steady state distribution, presented in Fig. 10(c)-(d), visually indistinguishable. We perform an analysis of coalescence by counting population densities present in different zones, under the two schemes described previously. Results are consistent regardless of zoning scheme, and we therefore only present the data for the peakcentred scheme. In Fig. 11(a) we present bar charts, which display the number of zones (at 20 h) where the zone population density is at × 2, × 5, × 10 and × 20 the mean density level M, defined as the density if the population is uniformly distributed across the landmass. Local/short range sampling generates numerous sites with moderately raised densities, but very few with significantly increased densities. For abundant populations these moderate increases may offer mating benefits, but for rare species they may have minimal effect on the overall encounter rate. Longer range sampling, however, generates relatively few sites with raised densities, but those sites exhibit substantial increases in population

density. The conferred advantage is likely to be particularly acute for the rarest populations, substantially raising densities and increasing the likelihood of encounters. We further calculate the normalised encounter rate, which we define as

ER(t ) =

Z  ui (t )2 . u2 i=1 max

In the above, Z is the total number of zones, ui (t) is the density of individuals within zone i at time t and umax is the maximum zone density (i.e. the density if the entire population is concentrated into a single zone). Note that the above assumes that male/female distributions are equivalent and that encounter rates are proportional to the product of male and female densities. The resulting value ranges between 0 and 1. Plotted as a function of time, nuances are revealed in terms of how the perceptual range impacts on the encounter rate. Over short time ranges lower to medium values of R appears to be more advantageous; at longer times, larger R values are optimal. Intuitively, larger perceptual range can encourage direct movements towards distant summits and less frequent on route encounters; lower to medium R may facilitate earlier encounters. This is reflected by comparing the trajectories indicated by Fig. 9(c) and (d). For R = 100 we observe trajectories (indicated by red dots) densely channelled along “virtual corridors” (Pe’er et al., 2005), which subsequently become dispersed for R = 500. Given the short lifespans of many insects, and the potential need to mate rapidly, sampling across a very large portion of the environment may therefore not always be the optimal strategy. 6. Discussion and conclusions The ability to detect information located distant to a cell or organism can be crucial for its navigation and movement. Larger organisms can achieve this through various means: vision and

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

hearing provide obvious, pervasive examples; electrolocation by fish through their lateral line organs presents a more specialised case (Hopkins, 2005). Certain migrating cells can extend long protrusions (filopodia, cytonemes) that span multiple cell diameters (Teddy and Kulesa, 2004). The majority of theoretical models, however, implicitly assume that the underlying cue is local (or effectively local) with respect to the individual’s position: for example, a chemical gradient computed across the individual’s basic body dimension. Here we have described a multiscale framework, based on a biased velocity-jump random walk, where the orienting environmental cue is either local or distant to the individual’s current location. Our underlying random walk is standard and, when scales of interest are suitably macroscopic, its behaviour can be described by a general form nonlinear drift-anisotropic diffusion equation (Hillen and Painter, 2013), the parameters of which depend on the navigation field/environmental cue. Via standard expansions of the modified Bessel functions, we show that in certain parameter regimes we can approximate the general form equation as a Keller–Segel type equation (for local gradient sensing) and a nonlocal integro-partial differential equation (under nonlocal sensing). Nonlocal integro-PDEs have become popular in movement models, in both ecological (Eftimie, 2012; Fagan et al., 2017; Mogilner and Edelstein-Keshet, 1999) and cellular systems (Hillen et al., 2007; Othmer and Hillen, 2002; Painter et al., 2015). The work here provides a means of connecting these models with microscopic/individual level behaviour; for other recent approaches see Middleton et al. (2014) and Buttenschoen et al. (2018). Of course, the form of continuous model is intrinsically linked to the VJRW assumptions and their alteration could demand recalculation. For example, the default assumption of exponentially distributed runtimes does not universally apply: under certain scenarios, occasional long transits are observed, a so-called Levy process described by long-tailed power-law distributions. Generalising VJRWs to include other runtime distributions or resting periods between runs has been considered in a number of studies (e.g. see Estrada-Rodriguez et al., 2018; Hillen, 2003; Taylor-King et al., 2015 amongst others). Even more challenging is the assumption of stochastically-independent walkers: while reasonable for scattered/dispersed population, this becomes questionable if the individuals aggregate/cluster tightly. Explicit descriptions of (ecological) perceptual range have been included in a variety of theoretical/computational models, including agent-based (Pe’er and Kramer-Schadt, 2008), stochastic mechanistic resource selection models (Barnett and Moorcroft, 2008) and integro-PDEs of drift/diffusion type (Fagan et al., 2017). Our modelling here has similarities, in particular, to the latter two studies. Similar to Barnett and Moorcroft (2008), our modelling is founded on stochastic random walk movements of an individual evaluating landscape across a perceptual range; in Moorcroft and Barnett (2008) the authors similarly derive diffusion-advection approximations for their stochastic model. The integro-PDE model studied in Fagan et al. (2017) includes a drift component that features nonlocal gradient following to some resource, similar to the form (11) suggested as a first order approximation to our general model. The utility of the model was demonstrated through a specific application to hilltopping behaviour. Under both idealised and genuine terrain profiles we demonstrated that nonlocal sensing allows a population to overcome terrain noise, so that individuals transit from lower to higher peaks; local sensing, on the other hand, can lead to trapping of population members at local maxima close to their starting location. Since topographical data will typically possess numerous local maxima, this suggests that an effective hilltopping mechanism would be based on highly nonlocal sensing. Applied to the topographical data of the Bodega Marine

239

Reserve, nonlocal sensing substantially increases population densities at a subset of “prominent peaks” across the studied area. As such, the model broadly reproduces field findings on the distribution of tiger moths (Arctia virginalis) at the same site (GrofTisza et al., 2017). A direct translation between our results and the field findings of Grof-Tisza et al. (2017) would require additional refinements, for example accounting for the biases arising from larval site/vegetation distribution. Here we have not accounted for this, allowing us to focus exclusively on the role played by nonlocal assessment on population distribution. Perceptual range, rather than the exact mode of sensing or nonlinear order, appears to be the primary determinant in the number of congregation sites. Large perception ranges are likely to be particularly helpful for the scarcest populations, bringing the dispersed population into concentrations at a small number of sites. Yet, the advantages become less clear-cut at shorter timescales, where the more dispersed travel routes taken to reach distant sites may reduce the likelihood of earlier mating: effectively, the narrow “virtual corridors” (Pe’er et al., 2005) encouraged by lower perceptual ranges become dispersed for larger R. Of course, the present modelling neglects to account for adaptive strategies, where individuals may first preferentially move to a nearby (lower) maximum and subsequently move to one more distant if no encounters are formed. Such hypotheses can easily be tested through more sophisticated nonlocal sensing rules. Here we consider two prototype models for nonlocal sensing: a nonlinear gradient and a nonlinear maximum model. In the context of hilltopping, these have clear and direct interpretations regarding whether the bias is according to the highest elevation angle or elevation within perceptual range. Either model appears to be an effective mechanism for overcoming terrain roughness and generating fewer, higher populated accumulations, although we remark that our earlier idealised experiments suggest that the nonlinear maximum model may be slightly more effective for overcoming noise. From a human perspective, elevation angle is easily gauged but estimating true height would require an additional capacity for depth perception. Overall, though, these prototypes can serve as simple models for describing various forms of sensing, such as visual, auditory, or chemical information. However, focused applications would clearly demand closer inspection of their foundation. For example, the choice of uniform is somewhat simplistic, where we may expect a more gradual drop-off in the ability to detect at a distance. Clearly the more general form (7) allows for a lot of tailoring. The hilltopping application features a fixed environmental cue, however, often the cue will vary with time. This could occur independently of the population, such an odour transported by turbulent currents, or according to the population distribution, such as a cell population that internalises a chemical signal or whales communicating through song. A key study would be to extend the analysis here to examine how nonlocal sensing impacts on oriented movements under dynamic variation of the cue; steps in this direction have been made in Fagan et al. (2017), where the potential beneficial effects of nonlocal sensing were illustrated. Related to this, here we have restricted ourselves to consider only fixed noise. Time-varying noises in the cue, such as molecular fluctuations in a chemical signal, are of significant interest: fluctuations may further restrict the effectiveness of local gradient-based sensing mechanisms and highlight the importance of nonlocal sensing. Studies of this nature may provide insight into some key questions, including the impact of human activity on ecological populations. For example, the effective perceptual range of cetacean auditory communication has been heavily eroded by noise pollution, such as commercial fishing, naval sonars and oil/gas exploration (Payne and Webb, 1971; Weilgart, 2007). Detailed modelling of such problems via our framework may provide important

240

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

predictions on how such unwanted noise can impact on population dynamics.

STJ thanks the Australian Mathematics Society for the award of a Lift-Off Fellowship to visit KJP. KJP acknowledges the Dipartimento di Scienze Matematiche and the Politecnico di Torino for a Visiting Professorship award.

Fig. A.12 presents a comparison between a numerical approximation of the relevant ratios modified Bessel functions and the O(kn ) approximation of each of the ratios. For k < 2 the ratios are adequately approximated by at least one of the truncated power series and, as expected, including additional terms results in a better approximation. As k increases the distance between the approximations are less accurate, suggesting that it would be inappropriate to implement these approximations in (1) above some threshold k value.

Appendix A. Modified Bessel function expansions

A2. Large argument approximations

A1. Small argument expansions

For larger k we can instead exploit the following large argument approximations for the modified Bessel functions of first kind (Olver, 2010), where for order j we have

Acknowledgement

Under nonnegative integer orders, modified Bessel functions of the first kind have the following power series expansions (Olver, 2010)

 j

I j (k ) =

k 2

∞  m=0

1 m!( j + m + 1 )

 2m k 2

.

(A.1)

Given (5) and (6), we are particularly interested in j = 0, 1, 2 and the ratios I1 (k)/I0 (k), I12 (k )/I02 (k ) and I2 (k)/I0 (k). Using (A.1), we first approximate the modified Bessel functions where j = 0, 1, 2 as

k2 k4 k6 + + + ··· , 4 64 2284

(A.2)

k k3 k5 k7 + + + + ··· , 2 16 384 18256

(A.3)

k2 k4 k6 k8 + + + + ··· . 8 96 3056 182560

(A.4)

I0 (k ) = 1 +

I1 (k ) =

I2 (k ) =

I12 (k ) I02

(k )

=

8k2 + 2k4 + . . . , 32 + 16k2 + 3k4 + · · ·

I2 (k ) 24k + 2k + . . . = . I0 (k ) 192 + 48k2 + 3k4 + · · · 2

(A.5)

(A.8)

Again, we are primarily interested in ratios I1 (k)/I0 (k), I12 (k )/I02 (k ) and I2 (k)/I0 (k). As the summation component of (A.8) is less than one for large k, we can approximate the relevant ratios as

1 1 1 25 I1 (k ) = 1− − 2 − 3 − + ··· , I0 (k ) 2k 8k 8k 128k4 I12 (k ) I02 (k )

To obtain the relevant ratios, I1 (k)/I0 (k), I12 (k )/I02 (k ) and I2 (k)/I0 (k), we truncate the polynomials on the numerator and denominator such that terms up to O(kn ) are included. We obtain

I1 (k ) 32k + 4k3 + · · · , = I0 (k ) 64 + 16k2 + k4 + · · ·

 ∞ n  ek (−1 )n 

2 2 I j (k ) = √ 1+ 4 j − ( 2 m − 1 ) . n! ( 8k )n 2π k n=0 m=1

=1−

1 1 1 − 3 − 4 + ··· , k 8k 4k

2 I2 (k ) 1 1 1 = 1 − + 2 + 3 + 4 + ··· . I0 (k ) k k 4k 4k

(A.9)

(A.10)

(A.11)

A comparison between the approximations truncated at O(kn ) for n = 0, −1, · · · , −4 and the numerical approximation of each of the ratios of modified Bessel functions is presented in Fig. A.13. We observe that all approximations for I1 (k)/I0 (k) are indistinguishable from the numerical approximation for k > 101 , whereas for both I12 (k )/I02 (k ) and I2 (k)/I0 (k) the approximations are near exact if k > 102 . Interestingly, the O(k−1 ) − O(k−4 ) converge to each other before converging to the numerical solution, which suggests that there is limited benefit to higher order approximations for large k values.

(A.6) Appendix B. Equation approximations

4

(A.7)

In general, according to the size of k we can truncate the above calculations at the kn th term to obtain a reasonable approximation:

Here we only state approximations for model (M1); similar approximations can be formulated for (M2) and (M3) following an analogous process. Note that we suppress subscripts on the κ sensitivity coefficient for clarity of presentation.

Fig. A12. Comparison between the ratio of modified Bessel functions (red) (a) I1 (k)/I0 (k), (b) I12 (k )/I02 (k ), (c) I2 (k)/I0 (k), and the O(k) (orange), O(k2 ) (purple), O(k3 ) (green) and O(k4 ) (blue) approximations of the ratios. Note that the O(k2 ) and O(k3 ) approximations are the same for I12 (k )/I02 (k ) and I2 (k)/I0 (k). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

241

Fig. A13. Comparison between the ratio of modified Bessel functions (red) (a) I1 (k)/I0 (k), (b) I12 (k )/I02 (k ), (c) I2 (k)/I0 (k), and the O(1), O(k−1 ) (orange), O(k−2 ) (purple), O(k−3 ) (green) and O(k−4 ) (blue) approximations of the ratios. Note that the O(k−1 ) and O(k−2 ) approximations are the same for I12 (k )/I02 (k ). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

B1. (M1) under small arguments

Appendix C. Terrain noise

Assuming small k (shallow gradient ∇ E/small response coefficient κ ), we approximate (5) and (6) using the first few terms of the expansions (A.5)–(A.7). The series of O(k ), O(k2 ), O(k3 ), . . . expansions are as follows: O ( k1 )

O ( k2 )

O ( k3 )

O (k ) 4

a1

= s2κ ∇ E,

D1

2 = τ 2s I;

a2

= s2κ 1 −

D2

2 = τ 2s 1 −

a3

=

D3

=

a4

=

D4

=





−



κ 2 |∇ E |

I − τ s2

8+2κ 2 |∇ E |2 2κ 2 |∇ E |2



κ 2 |∇ E |2

4+2κ 2 |∇ E |2

−

κ 2 |∇ E |2



8+2κ 2 |∇ E |2

A;

sκ 2 1 − 16+4κ 2 |∇ E |2 ∇ E,   2 2 κ |∇ E | κ 2 |∇ E |2 τ s2 1 − κ 2 |∇ E |2 I − τ s2 − A; 2 2 2 2 2 2 2 8+2κ |∇ E | 8+2κ |∇ E |  4+2κ |∇ E | 8κ 2 |∇ E |2 +κ 4 |∇ E |4 sκ 2 1 − 64+16κ 2 |∇ E |2 +κ 4 |∇ E |4 ∇ E,   8κ 2 |∇ E |2 +2κ 4 |∇ E |4 τ s2 1 − 24κ 2 |∇ E |2 +2κ 4 |∇ E |4 I − τ s2 2 192+48κ 2 |∇ E |2 +3κ 4 |∇ E |4 32+16κ 2 |∇ E |2 +3κ 4 |∇ E |4 2 4 2 4 24κ

|∇ E | +2κ |∇ E |

A;

192+48κ 2 |∇ E |2 +3κ 4 |∇ E |4

. . .

. . .

∇ E, 2

2+κ 2 |∇ E |2 4+κ 2 |∇ E |2

. . .

For large arguments k (shallow gradient ∇ E/small response coefficient κ ), we can utilise the first few terms of the expansions (A.9)–(A.11). Substituting into 5(6) generates a sequence of O(k0 ), O(k−1 ), O(k−2 ), . . . expansions as follows: a0

E = s |∇ ∇ E| ,

D0

= 0;

a1

s 2κ|∇ E −1 E = ( κ|∇ E|| ) |∇ ∇ E| ,

D1

1 = τ s2 κ|∇ I − A ); E| (

a2

=

D2

= τ s2

O(k−1 )

O(k−2 )

O(k−3 )







s 8κ 2 |∇ E |2 −4κ|∇ E |−1





8κ 2 |∇ E |2

κ|∇ E |−1 κ 2 |∇ E |2

=

D3

= τ s2



∇E |∇ E | ,

(I − A ) +

8κ 3 |∇ E |3 8κ 2 |∇ E |2 −8κ|∇ E |−3 8κ 3 |∇ E |3



τ s2 2κ 2 |∇ E |2





I;

∇E |∇ E | ,

(I − A ) +

τ s2 2κ 2 |∇ E |2

s 128κ 4 |∇ E |4 −64κ 3 |∇ E |3 −16κ 2 |∇ E |2 −16κ|∇ E |−25

a4

=

D4

= τ s2



128κ 4 |∇ E |4 8κ 3 |∇ E |3 −8κ 2 |∇ E |2 −3κ|∇ E |−4

2 + 2τ s 2 2 κ |∇ E |

. . .





s 8κ 3 |∇ E |3 −4κ 2 |∇ E |2 −κ|∇ E |−1

a3

O(k−4 )

Supplementary material Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.jtbi.2018.10.031. References

B2. (M1) under large arguments

O ( k0 )

To generate smoothly varying noisy terrain we consider a form of fractal noise (Perlin, 1985). We first sample numbers from the normal distribution with mean zero and standard deviation one at m points. We interpolate these numbers onto a grid with 2(m − 1 ) + 1 points using spline interpolation, and take the first m points from the grid. This ensures that any two noise points are correlated. We repeat the process, generating n random values and interpolate these values onto a grid with 4(m − 1 ) + 1 points, taking the first m points from this new grid, and adding this noise to the previously generated noise. This process is repeated log2 (y/(2f y)) times. Here we take f = 3. Finally we apply a three point hat filter to smooth the noise further and scale it such that it lies between − 0.5 and 0.5.

. . .

. . .



8κ 4 |∇ E |4

4κ 2 |∇ E |2 +2κ|∇ E |−3 4κ 2 |∇ E |2



I;



(I − A )



2κ|∇ E |+1 2κ|∇ E |

∇E |∇ E | ,



I;

Alcock, J., 1987. Leks and hilltopping in insects. J. Nat. Hist. 21, 319–328. Barnett, A., Moorcroft, P., 2008. Analytic steady-state space use patterns and rapid computations in mechanistic home range analysis. J. Math. Biol. 57, 139–159. Binny, R.N., Haridas, P., James, A., Law, R., Simpson, M., Plank, M., 2016. Spatial structure arising from neighbour-dependent bias in collective cell movement. PeerJ 4, e1689. Buttenschoen, A., Hillen, T., Gerisch, A., Painter, K.J., 2018. A space-jump derivation for non-local models of cell–cell adhesion and non-local chemotaxis. J. Math. Biol. 76, 429–456. Codling, E.A., Plank, M.J., Benhamou, S., 2008. Random walk models in biology. J. Roy. Soc. Interface 5, 813–834. Duistermars, B.J., Chow, D.M., Frye, M.A., 2009. Flies require bilateral sensory input to track odor gradients in flight. Curr. Biol. 19, 1301–1307. Eftimie, R., 2012. Hyperbolic and kinetic models for self-organized biological aggregations and movement: a brief review. J. Math. Biol. 65, 35–75. Ehrlich, P., Wheye, D., 1986. “Nonadaptive” hilltopping behavior in male checkerspot butterflies (euphydryas editha). Amer. Nat. 127, 477–483. Estrada-Rodriguez, G., Gimperlein, H., Painter, K.J., 2018. Fractional Patlak–Keller–Segel equations for chemotactic superdiffusion. SIAM J. Appl. Math. 78, 1155–1173. Fagan, W.F., Gurarie, E., Bewick, S., Howard, A., Cantrell, R.S., Cosner, C., 2017. Perceptual ranges, information gathering, and foraging success in dynamic landscapes. Am. Nat. 189, 474–489. Flaherty, E.A., Smith, W.P., Pyare, S., Ben-David, M., 2008. Experimental trials of the northern flying squirrel (glaucomys sabrinus) traversing managed rainforest landscapes: perceptual range and fine-scale movements. Can. J. Zool. 86, 1050–1058. Fletcher, R.J., Maxwell, C.W., Andrews, J.E., Helmey-Hartman, W.L., 2013. Signal detection theory clarifies the concept of perceptual range and its relevance to landscape connectivity. Landsc. Ecol. 28, 57–67. Friedrich, R., Jenko, F., Baule, A., Eule, S., 2006. Exact solution of a generalized Kramers–Fokker–Planck equation retaining retardation effects. Phys. Rev. E 74, 041103. Gould, J.L., Gould, C.G., 2012. Nature’s Compass: The Mystery of Animal Navigation. Princeton University Press.

242

S.T. Johnston, K.J. Painter / Journal of Theoretical Biology 460 (2019) 227–242

Grof-Tisza, P., Steel, Z., Cole, E.M., Holyoak, M., Karban, R., 2017. Testing predictions of movement behaviour in a hilltopping moth. Anim. Behav. 133, 161–168. Hillen, T., 2003. Transport equations with resting phases. Eur. J. Appl. Math. 14 (5), 613–636. Hillen, T., Painter, K., Schmeiser, C., 2007. Global existence for chemotaxis with finite sampling radius. Disc. Cont. Dyn. Sys. B 7, 125–144. Hillen, T., Painter, K.J., 2013. Dispersal, Individual Movement and Spatial Ecology: A Mathematical Perspective. Springer, pp. 177–222. Transport and anisotropic diffusion models for movement in oriented habitats Hillen, T., Painter, K.J., Swan, A.C., Murtha, A.D., 2017. Moments of von Mises and Fisher distributions and applications. Math. Biosci. Eng. 14 (3), 673–694. Hopkins, C.D., 2005. Passive electrolocation and the sensory guidance of oriented behavior. In: Electroreception. Springer, pp. 264–289. Jeon, N.L., Baskaran, H., Dertinger, S., Whitesides, G.M., Van De Water, L., Toner, M., 2002. Neutrophil chemotaxis in linear and complex gradients of interleukin-8 formed in a microfabricated device. Nat. Biotech. 20, 826–830. Keller, E., Segel, L., 1970. Initiation of slime mold aggregation viewed as an instability. J. Theor. Biol. 26, 399–415. Kornberg, T.B., Roy, S., 2014. Cytonemes as specialized signaling filopodia. Development 141, 729–736. Lima, S.L., Zollner, P.A., 1996. Towards a behavioral ecology of ecological landscapes. Trends Ecol. Evol. 11, 131–135. Mardia, K.V., Jupp, P.E., 20 0 0. Directional Statistics. Wiley, New York. McComb, K., Reby, D., Baker, L., Moss, C., Sayialel, S., 2003. Long-distance communication of acoustic cues to social identity in african elephants. Anim. Behav. 65, 317–329. Middleton, A.M., Fleck, C., Grima, R., 2014. A continuum approximation to an off-lattice individual-cell based model of cell migration and adhesion. J. Theor. Biol. 359, 220–232. Mogilner, A., Edelstein-Keshet, L., 1999. A non-local model for a swarm. J. Math. Biol. 38, 534–570. Moorcroft, P.R., Barnett, A., 2008. Mechanistic home range models and resource selection analysis: a reconciliation and unification. Ecology 89, 1112–1119. Olden, J.D., Schooley, R.L., Monroe, J.B., Poff, N.L., 2004. Context-dependent perceptual ranges and their relevance to animal movements in landscapes. J. Anim. Ecol. 73, 1190–1194. Olver, F.W.J., 2010. NIST Handbook of Mathematical Functions. Cambridge University Press. Othmer, H.G., Dunbar, S.R., Alt, W., 1988. Models of dispersal in biological systems. J. Math. Biol. 26, 263–298. Othmer, H.G., Hillen, T., 2002. The diffusion limit of transport equations ii: chemotaxis equations. SIAM J. Appl. Math. 62, 1222–1250. Painter, K.J., 2014. Multiscale models for movement in oriented environments and their application to hilltopping in butterflies. Theor. Ecol. 7, 53–75.

Painter, K.J., 2018. Mathematical models for chemotaxis and their applications in self-organisation phenomena. J. Theor. Biol. Article in Press. Painter, K.J., Bloomfield, J.M., Sherratt, J.A., Gerisch, A., 2015. A nonlocal model for contact attraction and repulsion in heterogeneous cell populations. Bull. Math. Biol. 77, 1132–1165. Patlak, C.S., 1953. Random walk with persistence and external bias. Bull. Math. Biol. 15, 311–338. Payne, R., Webb, D., 1971. Orientation by means of long range acoustic signaling in baleen whales. Ann. N.Y. Acad. Sci. 188, 110–141. Pe’er, G., Heinz, S., Frank, K., 2006. Connectivity in heterogeneous landscapes: analyzing the effect of topography. Land. Ecol. 21, 47–61. Pe’er, G., Kramer-Schadt, S., 2008. Incorporating the perceptual range of animals into connectivity models. Ecol. Model. 213, 73–85. Pe’er, G., Saltz, D., Frank, K., 2005. Virtual corridors for conservation management. Conserv. Biol. 19, 1997–2003. Pe’er, G., Saltz, D., Münkemüller, T., Matsinos, Y.G., Thulke, H.-H., 2013. Simple rules for complex landscapes: the case of hilltopping movements and topography. Oikos 122, 1483–1495. Pe’er, G., Saltz, D., Thulke, H.H., Motro, U., 2004. Response to topography in a hilltopping butterfly and implications for modelling nonrandom dispersal. Anim. Behav. 68, 825–839. Perlin, K., 1985. An image synthesizer. ACM SIGGRAPH Comput. Gr. 19, 287–296. Perthame, B., 2006. Transport Equations in Biology. Springer Science & Business Media. Prevedello, J.A., Forero-Medina, G., Vieira, M.V., 2011. Does land use affect perceptual range? evidence from two marsupials of the atlantic forest. J. Zool. 284, 53–59. Scott, J., 1968. Hilltopping as a mating mechanism to aid the survival of low density species. J. Res. Lepid. 7, 191–204. Shields, O., 1967. Hilltopping. J. Res. Lepid 6, 69–178. Taylor-King, J., van Loon, E., Rosser, G., Chapman, S., 2015. From birds to bacteria: generalised velocity jump processes with resting states. Bull. Math. Biol. 77, 1213–1236. Teddy, J.M., Kulesa, P.M., 2004. In vivo evidence for short-and long-range cell communication in cranial neural crest cells. Development 131, 6141–6151. Todd, L., Poulin, R., Brigham, R., Bayne, E., Wellicome, T., 2007. Pre-migratory movements by juvenile burrowing owls in a patchy landscape. Avian Cons. Ecol. 2. Wang, Y.-S., Potts, J., 2017. Partial differential equation techniques for analysing animal movement: a comparison of different methods. J. Theor. Biol. 416, 52–67. Weilgart, L.S., 2007. The impacts of anthropogenic ocean noise on cetaceans and implications for management. Can. J. Zool. 85, 1091–1116. Zollner, P.A., Lima, S.L., 1997. Landscape-level perceptual abilities in white-footed mice: perceptual range and the detection of forested habitat. Oikos 80, 51–60.