*Article* **Directionally Correlated Movement Can Drive Qualitative Changes in Emergent Population Distribution Patterns**

#### **Jonathan R. Potts**

School of Mathematics and Statistics, University of Sheffield, Hicks Building, Hounsfield Road, Sheffield S3 7RH, UK; j.potts@sheffield.ac.uk

Received: 26 June 2019; Accepted: 16 July 2019; Published: 18 July 2019

**Abstract:** A fundamental goal of ecology is to understand the spatial distribution of species. For moving animals, their location is crucially dependent on the movement mechanisms they employ to navigate the landscape. Animals across many taxa are known to exhibit directional correlation in their movement. This work explores the effect of such directional correlation on spatial pattern formation in a model of between-population taxis (i.e., movement of each population in response to the presence of the others). A telegrapher-taxis formalism is used, which generalises a previously studied diffusion-taxis system by incorporating a parameter *T*, measuring the characteristic time for directional persistence. The results give general criteria for determining when changes in *T* will drive qualitative changes in the predictions of linear pattern formation analysis for *N* ≥ 2 populations. As a specific example, the *N* = 2 case is explored in detail, showing that directional correlation can cause one population to 'chase' the other across the landscape while maintaining a non-constant spatial distribution. Overall, this study demonstrates the importance of accounting for directional correlation in movement for understanding both quantitative and qualitative aspects of species distributions.

**Keywords:** animal movement; correlated random walk; movement ecology; population dynamics; taxis; telegrapher's equation

#### **1. Introduction**

Understanding spatial distributions of animal species is a key concern for ecology, being of fundamental importance for a wide range of applications including conservation efforts [1], quantifying biodiversity [2], and invasive species research [3,4]. For mobile animals, decisions about where to move drive their individual locations. Consequently there has been a huge amount of effort in recent years to understand the causes and consequences of animal movement [5–7]. However, these movement decisions have an effect not only on individuals' locations but also on the spatial distribution of the whole population [8]. To understand this individual-to-population upscaling in a non-speculative way requires mathematical models that are built from the underlying movement processes of individuals [9]. Such models exhibit emergent phenomena on the population level that can be then quantified and related concretely to the underlying mechanisms of individual movement (e.g., [10,11]).

One recent study in this general area examines how the movement responses between individuals from *N* ≥ 2 different animal populations can drive spatio-temporal distribution patterns on temporal scales whereby births and deaths are negligible [12]. This takes spatial ecology in a slightly different direction to its tradition trajectory, whereby the combination of nonlinear birth-and-death terms (a.k.a. kinetics) combine with diffusive or cross-diffusive movement to drive spatio-temporal patterns [13–17]. Instead, the study of [12] shows that inter-population taxis (a form of cross-diffusion) can drive a wide range of complex patterns on its own, without the need for nonlinear kinetics.

However, the study by [12] assumes that individuals move diffusively in the absence of interactions. While this may be a reasonable approximation in many cases, in reality animals will always display some directional correlation in movement [18,19], if only due to the significant energetic costs of turning [20]. For some organisms, this correlation may only persist over quite small spatio-temporal scales, in which case diffusion is a reasonable model [21,22]. However, if the spatial scale of correlation is at a similar order of magnitude to the scale over which animals move in the time between successive interactions, diffusive assumptions may be less valid [23].

Here, I address this issue of directional correlation by extending the model of [12] from a diffusion-taxis system to a telegrapher-taxis system in 1D. The telegrapher's equation accurately models random movement with directional correlation in a 1D setting [24,25]. Moreover, it is a direct generalisation of the diffusion equation, simply requiring the introduction of a characteristic time-scale parameter, *T*, which is set to *T* = 0 in the diffusion limit.

The aim of this study is to examine the effect of introducing *T* > 0 on the linear pattern formation properties of the model in [12]. These patterning properties separate parameter space into three well-known categories [26]. First, patterns may not form at all from small non-constant perturbations of the steady state, i.e., the system is *stable* to such perturbations. Second, small non-constant perturbations of the steady state may grow over small times in a non-oscillatory fashion. This is classically known as a *Turing instability* after [27]. Third, these perturbations may both grow in magnitude and oscillate, which is sometimes known as a *Turing-Hopf instability*.

I give general criteria for when an increase *T* can cause a shift from one of these three patterning regions to another. When this shift occurs, it holds for all *T* > *T*<sup>∗</sup> where *T*<sup>∗</sup> is a threshold persistence time, which is quantified. As an example, I examine in detail the case *N* = 2 in the absence of self-aggregation. Here, the diffusion-taxis model (*T* = 0) can either be stable to non-constant perturbations or exhibit a Turing instability, dependent on the parameter values of the model. Furthermore, the parameter regimes where the system falls into these two patterning regions is known precisely [12]. However, the *N* = 2 and *T* = 0 case is never susceptible to a Turing-Hopf instability unless a self-aggregation process is in play [28]. I show that when *T* is increased sufficiently high, Turing-Hopf instabilities are possible without self-aggregation. Furthermore, they occur when the populations are engaged in a "pursuit-and-avoid" situation, with one population exhibiting taxis towards the other, and the latter exhibiting taxis away from the first. In contrast, for *T* = 0, this pursuit-and-avoid situation is always linearly stable to non-constant perturbations.

#### **2. The Modelling Framework and General Results**

This work focusses on a system of *N* populations, each of fixed size (i.e., no births or deaths). Individuals from each population move on a 1D landscape as biased, correlated random walkers (i.e., those where there is both persistence in movement and bias in a particular direction). This bias depends on the density of the various populations in the system, and is represented by a taxis term up or down the density gradient of the various populations. The correlated aspect of movement is modelled using a telegrapher's equation formalism [24,25].

Denoting by *ui*(*x*, *t*) the spatial probability density function of population *i* at time *t*, the study system is given by the following telegrapher-taxis equation for each population *i* (*i* ∈ {1, . . . , *N*})

$$T\frac{\partial^2 u\_i}{\partial t^2} + \frac{\partial u\_i}{\partial t} = d\_i \frac{\partial^2 u\_i}{\partial \mathbf{x}^2} - \frac{\partial}{\partial \mathbf{x}} \left[ u\_i \sum\_{j=1}^N \gamma\_{ij} \frac{\partial}{\partial \mathbf{x}} (\mathcal{K} \* u\_j) \right],\tag{1}$$

where *<sup>T</sup>* <sup>≥</sup> 0, *di* <sup>&</sup>gt; 0, *<sup>γ</sup>ij* <sup>∈</sup> <sup>R</sup>, and one can assume, without loss of generality, that the units are dimensionless and *d*<sup>1</sup> = 1. The case where *T* = 0 and *γii* = 0 for all *i* was studied in [12], and the reader is referred there for details of the non-dimensionalisation process. The dynamics take place on a unit line segment, [0, 1], with periodic boundary conditions

$$
u\_i(0, t) = 
u\_i(1, t),\tag{2}$$

so locations, *x*, live in the quotient space [0, 1]/{0, 1}. In Equation (1), K(*x*) is an integrable function, symmetric about *x* = 0 on the domain [0, 1]/{0, 1}, and K ∗ *uj* is the following spatial convolution

$$
\mathcal{K} \* \mu\_j(\mathbf{x}, t) = \int\_0^1 \mathcal{K}(\mathbf{x} - y) u\_j(y, t) \, \mathrm{d}y. \tag{3}
$$

The inter-population taxis term, in the right-hand summand of Equation (1), can come about through different biological mechanisms. The simplest is for the organisms in population *i* to sense directly the gradient of population density, *<sup>∂</sup> <sup>∂</sup><sup>x</sup>* (K ∗ *uj*), a mechanism that may be valid for very small individuals such as single-celled organisms or swarming insects. However, for larger organisms, this population gradient is more likely to be observed indirectly. For example, the deposition of marks in the environment from population *j* (e.g., through scenting; [29,30]) may indicate the population density, *uj*, across space. Similarly, memory of past interactions with individuals from population *j* can act as a proxy for sensing the population density gradient [31,32]. In [12], the authors showed how all of these biological mechanism can be modelled via the same taxis term, given in Equation (1).

The linear pattern formation properties of Equation (1) are analysed by perturbing the system about the constant steady-state solution, *ui*(*x*, *t*) = 1 for all *i*, *x*, *t*. Specifically, let **w**(*x*, *t*)=(*u*<sup>1</sup> − 1, ... , *uN* <sup>−</sup> <sup>1</sup>) = (*u*(0) <sup>1</sup> , ... , *<sup>u</sup>*(0) *<sup>N</sup>* ) exp(*σ<sup>t</sup>* <sup>+</sup> <sup>i</sup>*κx*), where *<sup>u</sup>*(0) <sup>1</sup> , ... , *<sup>u</sup>*(0) *<sup>N</sup>* <sup>∈</sup> <sup>R</sup>, *<sup>κ</sup>* <sup>∈</sup> <sup>R</sup>≥<sup>0</sup> and *<sup>σ</sup>* <sup>∈</sup> <sup>C</sup> are constants, and denotes matrix transpose. By neglecting non-linear terms, Equation (1) becomes

$$(T\sigma^2 + \sigma)\mathbf{w} = \kappa^2 M(\kappa)\mathbf{w},\tag{4}$$

where *M*(*κ*)=[*Mij*(*κ*)]*i*,*<sup>j</sup>* is an *N* × *N* matrix with

$$M\_{i\hat{j}}(\kappa) = \begin{cases} -d\_i + \gamma\_{ii}\hat{\mathcal{K}}(\kappa), & \text{if } i = j, \\ \gamma\_{i\hat{j}}\hat{\mathcal{K}}(\kappa), & \text{otherwise,} \end{cases} \tag{5}$$

and <sup>K</sup><sup>ˆ</sup> (*κ*) is the Fourier transform of <sup>K</sup>(*x*) on [0, 1]/{0, 1}.

Let *λ*1(*κ*), ... , *λN*(*κ*) be the eigenvalues of *M*(*κ*) (which are not necessarily distinct). If *T* = 0 then *σ* = *κ*2*λi*(*κ*) gives a solution to Equation (4) for some non-trivial vector **w**. From the perspective of pattern formation, there are three regimes of interest. These are well-documented [26] but it is valuable to summarise them briefly, for the purposes of introducing the key concepts and nomenclature used throughout this work:


The first of these states that the linear perturbation, **w**(*x*, *t*), will decay back to the homogeneous steady state, the second that **w**(*x*, *t*) will grow at small times for certain wavenumbers, *κ*, in a non-oscillatory fashion, and the third that **w**(*x*, *t*) will grow and oscillate at small times. These regimes give an indication for the pattern formation properties of the system. In the Stable region, the expectation is that spatial patterns do not form. In the Turing instability region, stationary patterns are predicted to form, which are fixed over time. If there is a Turing-Hopf instability, patterns that are in perpetual flux are expected to emerge. However, it is important to note that these are merely predictions, arising from a linearisation of the system, and that it is possible for the non-linear terms to cause different pattern formation properties asymptotically.

The main aim of this work is to examine how the demarcation into the above three regimes changes as *T* is increased from 0. When *T* > 0, for each *i* ∈ {1, ... , *N*}, there are up to two values of *σ* such that *Tσ*<sup>2</sup> + *σ* = *κ*2*λi*(*κ*). Thus I define

$$
\sigma\_i^{\pm}(\mathbf{x}, T) = \frac{-1 \pm \sqrt{1 + 4T\mathbf{x}^2 \lambda\_i(\mathbf{x})}}{2T},
\tag{6}
$$

for each *i* ∈ {1, ... , *N*}, and note that *σ*<sup>±</sup> *<sup>i</sup>* (*κ*, *T*) solves Equation (4) for each *i*. The main results from this work are given in the following three theorems.

**Theorem 1.** *In the case where the system is linearly stable for T* = 0 *(i.e., Re*(*λi*(*κ*)) < 0 *for all i* ∈ {1, ... , *N*}*, κ* > 0*) one of two situations can occur:*


**Proof.** For part (1), if *<sup>λ</sup>i*(*κ*) <sup>∈</sup> <sup>R</sup> for all *<sup>i</sup>*, *<sup>κ</sup>* then, since Re(*λi*(*κ*)) <sup>&</sup>lt; 0, the inequality 1 <sup>+</sup> <sup>4</sup>*Tκ*2*λi*(*κ*) <sup>&</sup>lt; <sup>1</sup> holds, so Re(*σ*± *<sup>i</sup>* (*κ*, *T*)) < 0 for all *i*, *κ*, *T*.

For part (2), let *<sup>i</sup>* ∈ {1, ... , *<sup>N</sup>*} and *<sup>κ</sup>* <sup>&</sup>gt; 0 such that *<sup>λ</sup>i*(*κ*) <sup>∈</sup>/ <sup>R</sup>. Assume, without loss of generality, that Im(*λi*(*κ*)) > 0 (otherwise, pick the complex conjugate of *λi*(*κ*)). Then if Re( <sup>1</sup> <sup>+</sup> <sup>4</sup>*Tκ*2*λi*(*κ*)) <sup>&</sup>gt; 1, the inequality Re(*σ*<sup>+</sup> *<sup>i</sup>* (*κ*, *T*)) > 0 holds, so that there is a Turing-Hopf instability.

I now show that if *T* is arbitrarily large, the criterion Re( <sup>1</sup> <sup>+</sup> <sup>4</sup>*Tκ*2*λi*(*κ*)) <sup>&</sup>gt; 1 is always satisfied. Here,

$$\operatorname{Re}(\sqrt{1+4T\kappa^2\lambda\_i(\kappa)}) \approx \operatorname{Re}(\sqrt{4T\kappa^2\lambda\_i(\kappa)}).\tag{7}$$

Now, arg( 4*Tκ*2*λi*(*κ*)) =arg(*λi*(*κ*))/2. Since Im(*λi*(*κ*)) <sup>&</sup>gt; 0 and Re(*λi*(*κ*)) <sup>&</sup>lt; 0, the following holds

$$
\pi/2 < \arg(\lambda\_i(\kappa)) < \pi. \tag{8}
$$

Therefore

$$0 < \arg(\sqrt{4T\kappa^2\lambda\_i(\kappa)}) < \pi/2,$$

so that Re(4*Tκ*2*λi*(*κ*)) <sup>&</sup>gt; 0. Hence Re( <sup>1</sup> <sup>+</sup> <sup>4</sup>*Tκ*2*λi*(*κ*)) <sup>&</sup>gt; 0 whenever *<sup>T</sup>* is sufficiently large. Furthermore, Re( <sup>1</sup> <sup>+</sup> <sup>4</sup>*Tκ*2*λi*(*κ*)) <sup>→</sup> <sup>∞</sup> as *<sup>T</sup>* <sup>→</sup> <sup>∞</sup>, so there exists some *<sup>T</sup>i*,*<sup>κ</sup>* <sup>∗</sup> such that Re( <sup>1</sup> <sup>+</sup> <sup>4</sup>*Tκ*2*λi*(*κ*)) <sup>&</sup>gt; 1 for all *<sup>T</sup>* <sup>&</sup>gt; *<sup>T</sup>i*,*<sup>κ</sup>* <sup>∗</sup> . There may be more than one *j* ∈ {1, ... , *N*} and *κ* > 0 such that *<sup>λ</sup>j*(*κ*) <sup>∈</sup>/ <sup>R</sup>. Thus let *<sup>T</sup>*<sup>∗</sup> be the minimum of *<sup>T</sup>j*,*<sup>κ</sup>* <sup>∗</sup> over such *j* and all *κ*. Then *T*<sup>∗</sup> satisfies the requirements of the theorem.

**Theorem 2.** *Consider the Turing instability case for <sup>T</sup>* <sup>=</sup> <sup>0</sup> *(i.e., argmaxλi*(*κ*)[*Re*(*λi*(*κ*)] <sup>∈</sup> <sup>R</sup>>0*). For a given κ, let λi*(*κ*) *be the dominant eigenvalue of M*(*κ*)*. Then one of two situations can occur:*


 $\text{Re}(\sqrt{1+4\kappa^2 T \lambda\_j(\kappa)}) > \sqrt{1+4\kappa^2 T \lambda\_i(\kappa)}$  for some j. Then there is a Turing-Hopf instability for all  $T > T\_\*$ .

**Proof.** For part (1), let *κ*, *T* > 0. If Re( <sup>1</sup> <sup>+</sup> <sup>4</sup>*κ*2*Tλj*(*κ*)) <sup>≤</sup> <sup>1</sup> <sup>+</sup> <sup>4</sup>*κ*2*Tλi*(*κ*) for all *<sup>j</sup>* then *<sup>σ</sup>*<sup>+</sup> *<sup>i</sup>* (*κ*, *T*) is the dominant eigenvalue. Since *<sup>λ</sup>i*(*κ*) <sup>∈</sup> <sup>R</sup>>0, <sup>1</sup> <sup>+</sup> <sup>4</sup>*κ*2*Tλi*(*κ*) <sup>∈</sup> <sup>R</sup>><sup>1</sup> so *<sup>σ</sup>*<sup>+</sup> *<sup>i</sup>* (*κ*, *<sup>T</sup>*) <sup>∈</sup> <sup>R</sup>><sup>0</sup> so there is a Turing instability for these values of *T* and *κ*.

For part (2), suppose there exist *j*, *T*, *κ* with Re( <sup>1</sup> <sup>+</sup> <sup>4</sup>*κ*2*Tλj*(*κ*)) <sup>&</sup>gt; <sup>1</sup> <sup>+</sup> <sup>4</sup>*κ*2*Tλi*(*κ*) and suppose Re( 1 + 4*κ*2*Tλj*(*κ*)) ≥ Re( <sup>1</sup> <sup>+</sup> <sup>4</sup>*κ*2*Tλk*(*κ*)) for all *<sup>k</sup>*. Then *<sup>σ</sup>*<sup>+</sup> *<sup>j</sup>* (*κ*, *T*) is the dominant eigenvalue. It is necessary to show that *σ*<sup>+</sup> *<sup>j</sup>* (*κ*, *<sup>T</sup>*) is not real. Therefore, for a contradiction, suppose 1 + 4*κ*2*Tλj*(*κ*) ∈ <sup>R</sup>. Then *<sup>λ</sup>j*(*κ*) <sup>∈</sup> <sup>R</sup>. However, *<sup>λ</sup>j*(*κ*) <sup>≤</sup> *<sup>λ</sup>i*(*κ*), since *<sup>λ</sup>i*(*κ*) is the dominant eigenvalue of *<sup>M</sup>*(*κ*), so <sup>1</sup> <sup>+</sup> <sup>4</sup>*κ*2*Tλj*(*κ*) <sup>≤</sup> <sup>1</sup> <sup>+</sup> <sup>4</sup>*κ*2*Tλi*(*κ*), which contradicts the assumption. Hence *<sup>σ</sup>*<sup>+</sup> *<sup>j</sup>* (*κ*, *<sup>T</sup>*) <sup>∈</sup>/ <sup>R</sup> so there is a Turing-Hopf instability for these values of *T* and *κ*.

*Corollary.* If there is some *j* such that Re( *<sup>λ</sup>j*(*κ*)) > *λi*(*κ*) then there is a Turing-Hopf instability at wavenumber *κ* for sufficiently large *T*.

**Theorem 3.** *Consider the case where there is a Turing-Hopf instability for T* = 0*. Then there is a Turing-Hopf instability for all T* > 0*.*

**Proof.** Let *i* be such that *λi*(*κ*) is the dominant eigenvalue of *M*(*κ*) for some *κ* where there is a Turing-Hopf instability for *T* = 0. Assume, without loss of generality, that Im(*λi*(*κ*)) > 0 (otherwise choose the complex conjugate of *λi*(*κ*)). Since Re(*λi*(*κ*)) > 0, the inequality Re(1 + 4*Tκ*2*λi*(*κ*)) > 1 holds, so Re(<sup>1</sup> <sup>+</sup> <sup>4</sup>*Tκ*2*λi*(*κ*)) <sup>&</sup>gt; 1. Since 1 <sup>+</sup> <sup>4</sup>*Tκ*<sup>2</sup> *λi*(*κ*) has positive real and imaginary part, Re(1 + 4*Tκ*2*λi*(*κ*)) < Re( <sup>1</sup> <sup>+</sup> <sup>4</sup>*Tκ*2*λi*(*κ*)). The latter follows from the following general calculation for *r* > 0 and 0 < *θ* < *π*/2

$$
\sqrt{\operatorname{Re}(r\mathbf{e}^{i\theta})} = \sqrt{r}\sqrt{\cos(\theta)} < \sqrt{r}\sqrt{\frac{\cos(\theta) + 1}{2}} = \sqrt{r}\cos\left(\frac{\theta}{2}\right) = \operatorname{Re}(\sqrt{r\mathbf{e}^{i\theta}}).\tag{9}
$$

The inequality in (9) requires 0 < cos(*θ*) < 1, which follows from 0 < *θ* < *π*/2.

It follows that Re( <sup>1</sup> + <sup>4</sup>*Tκλi*(*κ*)) > 1, and thus from Equation (6) that Re(*σ*<sup>+</sup> *<sup>i</sup>* (*κ*, *T*)) > 0 for all *T* > 0. To show that this is within the Turing-Hopf instability region, it is necessary to check that there does not exist *j* such that both *σ*<sup>+</sup> *<sup>j</sup>* (*κ*, *<sup>T</sup>*) <sup>∈</sup> <sup>R</sup>><sup>0</sup> and *<sup>σ</sup>*<sup>+</sup> *<sup>j</sup>* (*κ*, *<sup>T</sup>*) <sup>&</sup>gt; Re(*σ*<sup>+</sup> *<sup>i</sup>* (*κ*, *T*)), since otherwise this is the Turing instability region. For a contradiction, suppose such a *σ*<sup>+</sup> *<sup>j</sup>* (*κ*, *T*) exists. Then

$$\sqrt{1 + 4T\kappa^2 \lambda\_j(\kappa)} > \text{Re}(\sqrt{1 + 4T\kappa^2 \lambda\_i(\kappa)}) > \sqrt{\text{Re}(1 + 4T\kappa^2 \lambda\_i(\kappa))}\tag{10}$$

so 1 + 4*Tκ*2*λj*(*κ*) > Re(1 + 4*Tκ*2*λi*(*κ*)) so *λj*(*κ*) > Re(*λi*(*κ*)), which contradicts the fact that *λi*(*κ*) is the dominant eigenvalue for the *T* = 0 case.

#### **3. The Case of Two Interacting Populations (***N* **= 2)**

In this section, I look in detail at the specific example where *N* = 2 and *γii* = 0 for all *i* ∈ {1, . . . , *N*}, to show how persistent movement can drive qualitative changes in pattern formation properties even in this simple situation. In the case *T* = 0 (studied by [12]), the following holds

$$\lambda\_{\pm}(\kappa) = \frac{-(1+d\_2) \pm \sqrt{(1-d\_2)^2 + 4\gamma\_{12}\gamma\_{21}\hat{C}^2(\kappa)}}{2},\tag{11}$$

and it follows that the system for *T* = 0 is linearly stable to wave perturbations whenever *γ*12*γ*<sup>21</sup> < *d*2. When *<sup>γ</sup>*12*γ*<sup>21</sup> <sup>&</sup>gt; *<sup>d</sup>*2, the system can have a Turing instability, and always does in the limit <sup>K</sup><sup>ˆ</sup> (*κ*) <sup>→</sup> 1, where the spatial averaging given by K is arbitrarily small. The system never exhibits a Turing-Hopf instability for *T* = 0 [12].

In the Turing instability case for *T* = 0, both the eigenvalues are real, so the conditions of Theorem 2, part 2, are met. Thus there is a Turing instability for all *T* > 0. The case where *γ*12*γ*<sup>21</sup> < *d*2, however, is more interesting. Here, Theorem 1 states that a Turing-Hopf instability will occur for sufficiently large *<sup>T</sup>* as long as *<sup>λ</sup>*±(*κ*) <sup>∈</sup>/ <sup>R</sup>, that is

$$
\hat{\lambda}(1-d\_2)^2 + 4\gamma\_{12}\gamma\_{21}\hat{\mathcal{K}}^2(\mathbf{x}) < 0. \tag{12}
$$

The first thing to notice about the condition given by (12) is that *γ*<sup>12</sup> and *γ*<sup>21</sup> need to be of opposite signs. This means that one population tends to advect up the density gradient of the other, while the second population retreats down the gradient of the first. This is termed a 'pursuit-and-avoid' situation in [12]. In this situation, for *T* = 0, no patterns form: the populations each settle to a steady-state distribution whereby they are uniformly distributed on the terrain. However, if *T*<sup>∗</sup> is defined to be the minimum positive real number such that there exists *κ* > 0 with Re( <sup>1</sup> + <sup>4</sup>*Tκ*2*λ*+(*κ*)) > 1, then Theorem 1 implies there is a Turing-Hopf instability whenever *T* > *T*∗.

To understand how *<sup>T</sup>*<sup>∗</sup> depends on the parameters, I fix *<sup>d</sup>*<sup>2</sup> <sup>=</sup> 1 and set *<sup>γ</sup>* = +√−*<sup>γ</sup>*12*γ*21. (Notice that *<sup>γ</sup>* <sup>∈</sup> <sup>R</sup>><sup>0</sup> since *<sup>γ</sup>*<sup>12</sup> and *<sup>γ</sup>*<sup>21</sup> are of opposite signs.) It is also necessary to pick a particular functional form for K(*x*). This is a slightly delicate matter as it is not necessarily the case that the pattern formation problem is well-posed for an arbitrary K(*x*). For example, if K(*x*) = *δ*(*x*) is used, where *δ*(*x*) is the Dirac delta function, then, by Equation (6), the dominant eigenvalue is

$$
\sigma\_i^+ (\kappa, T) = \frac{-1 + \sqrt{1 - 4T\kappa^2 + 4T\text{i}\kappa^2\gamma}}{2T}. \tag{13}
$$

If *T* is fixed and *κ* is arbitrarily large, then Re(*σ*<sup>+</sup> *<sup>i</sup>* (*κ*, *T*)) ≈ *κ*Re( √−<sup>1</sup> <sup>+</sup> <sup>i</sup>*γ*/ <sup>√</sup>*T*) <sup>→</sup> <sup>∞</sup> as *<sup>κ</sup>* <sup>→</sup> <sup>∞</sup>. Hence patterns can form for arbitrarily high wavenumbers, meaning the pattern formation problem is ill-posed.

However, suppose instead that

$$\mathcal{K}(\mathbf{x}) = \begin{cases} \frac{\exp\left(-\frac{\mathbf{r}^2}{\rho^2}\right)}{\sigma\sqrt{\pi}\text{erf}\left(\frac{1}{2\pi}\right)}, & \text{if } -1/2 < \mathbf{x} < 1/2\\ 0, & \text{otherwise} \end{cases} \tag{14}$$

for some *σ* > 0, so that

$$\hat{\mathcal{K}}(\kappa) = \frac{\text{erf}\left(\frac{1}{2\tau} - \frac{\text{i}\kappa\tau}{2}\right) + \text{erf}\left(\frac{1}{2\tau} + \frac{\text{i}\kappa\tau}{2}\right)}{2\text{erf}\left(\frac{1}{2\tau}\right)} \exp\left(\frac{-\kappa^2 \sigma^2}{4}\right). \tag{15}$$

The case of interest is where *σ* is small, in which case the following approximation can be made

$$
\hat{\mathcal{K}}(\kappa) \approx \exp\left(\frac{-\kappa^2 \sigma^2}{4}\right). \tag{16}
$$

Then Equation (6) gives

$$\sigma\_i^+\left(\mathbf{x}, T\right) = \frac{-1 + \sqrt{1 - 4T\kappa^2 + 4T\text{i}\gamma\kappa^2 \exp\left(-\kappa^2 \sigma^2 / 4\right)}}{2T}. \tag{17}$$

Now, notice that if *T* is fixed and *κ* is arbitrarily large, then

$$\operatorname{Re}(\sigma\_i^+(\kappa, T)) \approx \kappa \text{Re}\left(\frac{\sqrt{-1 + \operatorname{i}\gamma \exp(-\kappa^2 \sigma^2 / 4)}}{\sqrt{T}}\right) \to 0,\tag{18}$$

as *κ* → ∞. Therefore, with the Gaussian averaging kernel from Equation (14), the pattern formation problem is well-posed, in the sense that patterns cannot form at arbitrarily high wavenumbers.

By Theorem 1, there is some *T*<sup>∗</sup> such that there is a Turing-Hopf instability for all *T* > *T*∗. Furthermore, this Theorem states that *T*<sup>∗</sup> is the minimum *T* such that the following holds for some *κ*

$$f(\kappa) := \operatorname{Re}[\sqrt{1 - 4T\kappa^2 + 4T\text{i}\gamma\kappa^2 \exp(-\kappa^2\sigma^2/4)}] > 1. \tag{19}$$

In Figure 1a, *f*(*κ*) is plotted against *κ* for *γ* = 0.2, *σ* = 0.05, and varying values of *T*. This reveals that *T* not only affects whether patterns form, but the range of wavenumbers for which patterns may form. In this example, there is a Turing-Hopf bifurcation for value of *T* = *T*<sup>∗</sup> somewhere between *T* = 0.05 and *T* = 0.1.

The precise values of *T*<sup>∗</sup> for a range of *γ*- and *σ*-values are plotted in Figure 1b,c. The *γ* parameter encodes the strength of attraction/avoidance. Figure 1b,c shows that *T*<sup>∗</sup> decays for increasing *γ* and grows with the width of spatial averaging, *σ*.

**Figure 1. Critical value of** *T* **for pattern formation.** In Panel (**a**), Equation (19) is plotted for *γ* = 0.2, *σ* = 0.05, and various values of *T*. Where *f*(*κ*) > 1, there is a Turing-Hopf instability. In Panels (**b**,**c**), the value, *T*∗, above which there is a Turing-Hopf instability for some *κ* and below which there is not, is plotted for various values of *σ* and *γ*.

To understand the qualitative properties of the patterns that emerge from increasing *T*<sup>∗</sup> past the Turing-Hopf bifurcation point, the system in Equations (1)–(3) is solved numerically, with *N* = 2, *<sup>d</sup>*<sup>2</sup> <sup>=</sup> 1, *<sup>γ</sup>* = +√−*<sup>γ</sup>*12*γ*21, and the Gaussian spatial averaging from Equation (14). In Figure 2, the results of these numerics are shown for certain values of *γ*, *σ* and *T*. Here, rather than decaying to the homogeneous steady state (*T* = 0), for higher *T* the populations move across the landscape while maintaining a non-constant population distribution.

**Figure 2. Example numerics for** *N* = 2**.** Numerical solutions of Equation (1) for *N* = 2, *d*<sup>2</sup> = 1, *γ*<sup>12</sup> = 1, *γ*<sup>21</sup> = −1, and K(*x*) as in Equation (14) with *σ* = 0.1. Both examples have random continuous initial conditions, but forced to be symmetric about *x* = 0.5 to satisfy periodic boundary conditions. In Panel (**a**), *T* = 0, and the initial perturbation of the constant steady state decays to the solution *u*1(*x*, ∞) = 1. In Panel (**b**), *T* = 1. Here, the population *u*<sup>1</sup> moves across the landscape, not settling to a constant distribution.

#### **4. Discussion**

This study examines a telegrapher-taxis system of *N* populations to demonstrate how directional correlation can drive pattern formation in systems of between-population taxis. General criteria are given for switches in pattern formation regime driven by directional correlation and the *N* = 2 case is examined in detail. These results demonstrate the importance of considering directional correlation when seeking to understanding the spatio-temporal population distributions of moving organisms. As our ability to capture data on animal movement is becoming increasingly sophisticated, better understanding of how the details of animal movement can affect population dynamics is becoming an ever-more pertinent question [6], with implications for a diverse range of ecological areas such as connectivity dynamics [33] and conservation [34].

The model of [12], on which the present model is based, is closely related to aggregation and chemotaxis models inspired by cell biology [35–38]. Although many such models only examine a single population, there are various examples of two-population models [28,39], and some that incorporate an arbitrary number of populations [40,41]. In the latter examples, the diffusion-taxis equations are coupled via interactions with a diffusive chemical, different to the model studied here. While cells lack the momentum of much larger organisms, directional correlation is known to be a factor in the movement of cells in certain circumstances [42–44]. Therefore the results presented here suggest that it is worth exploring how directional correlation may affect the pattern formation properties of such aggregation and chemotaxis models.

Here, ecosystems are modelled on a constrained timescale whereby births and deaths are negligible. However, it would be valuable to extend the model presented here to incorporate such effects, via competition and/or predation terms (i.e., kinetics), and so explore the effect of directional persistence over longer timescales. Adding directional persistence to models that incorporate kinetics is non-trivial, though, and does not simply involve adding a second temporal derivative to a reaction-diffusion-taxis model [45].

Nonetheless, building such a system of reaction-telegrapher-taxis equations would enable connection with cross-diffusion models (e.g., [16,17,46,47]). These are generalisations of reactiondiffusion systems, whereby terms are added that allow for motion driven by the presence of foreign populations. These terms include, but are not limited to, the taxis terms observed in the model from Equation (1). However, incorporating directional correlation via a telegrapher's term into a cross-diffusive setting is much rarer (but see [48]) and the pattern formation properties of such systems are not well-explored. Given the results presented here, it may be interesting to extend cross-diffusive models in this way, to show how directional correlation may affect pattern formation in such models.

In summary, the simple but general results shown here demonstrate that directional correlation of individuals' movement can have a great effect on the spatio-temporal distribution of species. While I have demonstrated this in a model of ecosystems relevant over relatively short timescales, where births and deaths are minimal, it highlights a general principle that is little-studied and may have much wider implications.

**Acknowledgments:** JRP acknowledges support from the Natural Environment Research Council (NERC) (grant number NE/R001669/1).

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
