The effect of landscape fragmentation on Turing-pattern formation

1.   Introduction

Spatial patterns are ubiquitous in ecological systems and are widely investigated using reaction-diffusion equations that account for the random motion and population dynamics of the species involved. When the landscape itself is heterogeneous, e.g., in the distribution of abiotic conditions, it seems plausible that the distribution of the species that respond to these conditions is also heterogeneous, following the abiotic conditions. Turing's surprising discovery was that even in a perfectly homogeneous landscape, spatially patterned species distributions can arise from the interplay between species interaction and movement in space ^[1]. In fact, Turing gave conditions under which a spatially homogeneous steady state of two interacting species is stable in the absence of diffusion but becomes unstable in the presence of diffusion. This effect is often called diffusion-driven instability (DDI), and the resulting patterns are known as Turing patterns. The crucial ingredient in Turing's theory is the sign pattern of the interaction matrix (Jacobian matrix) of the two species at the homogeneous steady state: only if the sign pattern is of activator–inhibitor, i.e.,

or of positive feedback type ^[2,3], can DDI arise. It only does arise if the inhibitor species has a much higher diffusion rate than the activator.

Although Turing's study described chemical reactions, his ideas were soon applied to biological systems, in particular to developmental biology and animal-coat patterns ^[2]. Relatively few applications to ecology exist, because some of the conditions seem quite restrictive for ecological systems. An order of magnitude difference in dispersal rates between species is not easy to satisfy. And the Rosenzweig–MacArthur model, the most commonly used predator–prey model, cannot not exhibit DDI because its Jacobian matrix cannot have one of the required sign structures ^[4]. Nonetheless, examples exist, some of which we mention below. More recently, the search for Turing patterns in ecological systems has seen a revival ^[5,6]. In this work, we investigate how spatial heterogeneity in landscape characteristics (e.g., abiotic) can alter the conditions for pattern formation in an ecological predator–prey system. We find mechanisms that can lead to DDI and patterns in heterogeneous landscape where none would exist in a homogeneous landscape.

Segel and Jackson were the first to find DDI in ecological dynamics ^[7]. They chose the Lotka-Volterra predator–prey model with an additional Allee effect for the prey and a self-limitation term for the predator. For a certain range of parameters, the sign structure of the Jacobian matrix at the positive steady state for this model is of activator-inhibitor type with the prey as an activator. As a result, this model can give rise to Turing patterns if the prey disperse sufficiently less than the predator. Turing's mechanism was also invoked by Levin and Segel to explain the emergence of patchy population distributions of oceanic plankton ^[8]. If biological interactions between phytoplankton and herbivores reduce the efficiency of herbivory as phytoplankton density increases and if differential dispersal rates favour higher herbivore mortality, such patterns arise.

While the Rosenzweig–MacArthur model does not lead to Turing patterns, several other predator–prey models do, such as the Beddington–DeAngelis model ^[9], the ratio-dependent model ^[10] and the May model ^[11]. More precisely, if diffusion is added to these models, they can give rise to DDI and Turing patterns if the prey disperse sufficiently less than the predator ^[12,13,14]. Alternatively, certain variants of the Rosenzweig–MacArthur model can exhibit DDI and Turing patterns. Fasani and Rinaldi included various aspects of predator behavior into this model and gave conditions for when DDI can result ^[15]. Specifically, they chose small perturbations that would lead to a required sign pattern, in this case the positive feedback structure ^[2,3]. Consequently, DDI and Turing-pattern formation are possible when the predator disperses sufficiently less than the prey.

By way of example, we make these considerations more concrete based on the May model ^[11], which we also use later in our examples. In this model, the prey species $u = u(t)$ grows logistically and is preyed upon with a Holling type II functional response by a predator species $v = v(t)$. The predator species grows logistically with carrying capacity proportional to the prey density. The nondimensional equations read

Here, $b$ denotes the half-saturation constant in the functional response ^[16,17], $s$ is the intrinsic growth rate of the predator, and $qu$ its carrying capacity. All three parameters are assumed positive.

Eq (1.2) has a unique positive steady state, given by

The Jacobian matrix at this steady state is

Since the entries in the second column of this matrix are negative, the predator inhibits its own and the prey's growth at steady state. Since the (2, 1) entry is positive, the prey activates the predator's growth at steady state. The (1, 1) entry can change sign. If it is positive, the prey also activates its own growth, and we have the required activator–inhibitor sign pattern with the prey as the activator and the predator as the inhibitor. If it is negative, DDI is impossible ^[2,3].

Let us choose $q$ as our parameter of interest and fix the other two. Then the (1, 1) entry is positive if

For DDI, we also require that the positive steady state of the nonspatial model be stable. Since the determinant of the Jacobian matrix in Eq (1.4) is positive, the steady state of the May model is stable if the trace is negative. The resulting calculations are tedious but not difficult. One finds that if $s$ is small enough, then the trace is positive for some range $q_{H, 1} < q < q_{H, 2}$, where the thresholds $q_{H, i}$ can be calculated explicitly. Otherwise, the trace is negative and the steady state is stable. More details can be found in ^[18,19].

To consider spatial effects in the May model, we denote by $u = u(x, t)$ and $v = v(x, t)$ the densities of prey and predator in a one-dimensional spatial domain and we add diffusion terms to each of the variables in Eq (1.2). We denote the diffusion rates of the two species by $D^u$ and $D^v$, respectively. The steady state $(u^*, v^*)$ of the nonspatial model above becomes a spatially constant steady state of the spatial model. This state is stable with respect to spatially constant perturbations but can be unstable with respect to spatial perturbations of certain wavelengths if the Jacobian matrix has the sign patterns indicated above and if $D^v\gg D^u$. Specifically, there exists a threshold $D_c > 1$, such that DDI occurs for $D^v/D^u > D_c$ ^[2,3]; see also below.

The study of DDI and Turing patterns in spatially heterogeneous landscapes is much less developed and general insights are relatively rare. In part, this is because of the myriad of possibilities in which spatial heterogeneity can present itself in the model. In the simplest case, a heterogeneous landscape consists of only two patches with an abrupt change at the interface between them. Then the model parameters are piecewise constant, i.e., assume one value on each patch.

Benson et al. ^[20,21] considered such a scenario and studied Turing-pattern formation in embryological and ecological models with a general two-species reaction-diffusion model in one space dimension, where the domain is composed of two parts. They required that population densities and fluxes are continuous at the patch interface. Their analysis provides an understanding of how environmental heterogeneities can modulate the pattern-forming capabilities of reaction-diffusion systems. While Benson et al. concentrated on diffusion coefficients changing in space, Page et al. focused on the case of kinetic parameters varying in space ^[22]. They showed that a step-function heterogeneity in a kinetic parameter of a reaction–diffusion system can lead to spatial pattern formation outside the classical Turing parameter regime for patterning. They found that the resulting patterns are spatially localised around the parameter discontinuity because the discontinuity in the steady state values acts like a local perturbation and can induce pattern formation if the homogeneous steady state on either side of the boundary is unstable. In ^[23], Page et al. extended their model to smoothly varying kinetic parameters, either spatially monotone or periodic. More recently, Kozák and coauthors expanded the analysis of pattern formation with a step-function heterogeneity in the kinetic parameters and found different frequency and amplitude patterns on the two sides of the discontinuity ^[24].

Sheffer et al. studied the interaction between patterns that arise in response to landscape heterogeneity and patterns that arise from DDI ^[25]. They used theoretical and empirical approaches for the formation of patterned vegetation in semi-arid regions.They showed that both mechanisms significantly affect the pattern-formation process, with their relative contribution depending on water availability. They explored the effect of spatial scales on pattern formation by comparing the length scale of the spatial extent of landscape feature to the length scale on which biological feedbacks operate.

Cobbold et al. studied how extrinsic factors, environmental variation and intrinsic interaction can lead to spatial patterns in a predator–prey model in a heterogeneous landscape ^[18]. Their landscape consists of discrete patches of two different types that alternate in space. Population dynamics follow ordinary differential equations on each patch (e.g., the May model above) and movement is modeled by discrete, not necessarily symmetric, diffusion between patches. With a finite number of patches, they found a large number of patterns, often stably coexisting, and complex bifurcation diagrams through a numerical bifurcation analysis. With infinitely many patches, they derived the dispersion relation of the spatially homogeneous steady state and the corresponding pattern formation conditions.

We use a reaction-diffusion system to study these questions in a continuous space and time setting. Hence, our study is in some sense a continuous-space version of the model by Cobbold et al. ^[18]. Working with continuous space, however, makes our work not only closer to Turing's original ideas, it also allows us to explore several additional features, such as the influence of patch sizes, and use additional techniques, such as homogenization. Specifically, we consider periodically alternating patches of two types on the real line. On each patch, we consider a system of reaction-diffusion equations of predator and prey. Between adjacent patches, we have conditions that match population densities and fluxes across the interface. While population fluxes are continuous across an interface, their densities need not be if they represent certain movement behavior at the interface ^[26,27]. Our methods are completely different from those used in ^[18] as we employ averaging theory to derive a homogenized model under the assumption that the landscape period is small ^[28,29]. We derive the DDI conditions for this homogenized model. We use numerical evaluation to analyze and visualize the DDI conditions. We also simulate the homogenized model to test the prediction of the DDI conditions. Since parameter space is fairly large in our model, we carefully select three instructive scenarios to study how the interaction of spatial heterogeneity, movement behavior and population dynamics can lead to pattern formation in ranges of parameter space where we would not expect them (Sections 3–5). We note that some reaction-diffusion systems support patterned states that do not arise via bifurcations from a homogeneous state; we do not study such patterns.

2.   The model set-up and method

We present a general predator–prey model in a one-dimensional heterogeneous landscape that contains infinitely many patches: patches of Type 1, or "good patches", and patches of Type 2, or "bad patches", which are periodically alternating. We denote by $L_1$ and $L_2$ the length of good and bad patches, respectively, and by $L = L_1+L_2$ the spatial period. We also denote $l_1 = \frac{L_1}{L}$ as the fraction of good patches and $l_2 = \frac{L_2}{L}$ as the fraction of bad patches, so that $l_1+l_2 = 1$. Without loss of generality, we choose a good patch to be located at $(-L_1, 0)$ and other good patches $L$-periodic from thereon. Accordingly, bad patches are located at $(0, L_2)$ and $L$-periodic from thereon. The population densities for prey and predator on patch $i$ are denoted by $u_i(x, t)$ and $v_i(x, t)$, respectively. Then we describe the population dynamics and individual movement in our predator–prey model for each species on each patch as follows:

The equations for $u_1$ and $v_1$ hold on all good patches, i.e., for $x\in (-L_1, 0)+L\mathbb{Z}$, and analogously for $u_2$ and $v_2$ on bad patches. We denote the diffusion coefficients on patch type $i$ for species $u$ and $v$ by $D_i^u, D_i^v$, respectively. The interaction functions on patch type $i$ are given by $f_i$ and $g_i$. We have two types of interfaces between adjacent patches: At $x = 0$ we have a good patch on the left and a bad patch on the right; at $x = L_2$, the situation is reversed. Since the landscape is periodic, so are the interface locations. We need to impose matching conditions for population densities and fluxes at these interfaces. We choose the formulation by Ovaskainen and Cornell ^[26], discussed in detail by Maciel and Lutscher ^[27]. We denote by $p^u_i$ and $p^v_i$ the probabilities for prey and predator, respectively, that an individual at an interface moves to the patch of Type $i$ (with $p^{u, v}_1 = 1-p^{u, v}_2$). Then the matching conditions read:

where $k_u = \frac{p^u_1}{p^u_2}\frac{D^u_2}{D^u_1}$ and $k_v = \frac{p^v_1}{p^v_2}\frac{D^v_2}{D^v_1}$ are dimensionless parameters. Since they are in general different from unity, the population densities are discontinuous across an interface (first and third line in Eq (2.2)). Continuity of the fluxes, however, is guaranteed by the matching conditions in the second and fourth line in Eq (2.2).

The above model is quite difficult to study. Mathematically, the discontinuity at the interfaces requires careful consideration for the existence of solutions of the nonlinear problem ^[30] and eigenvalues for the linearized problem ^[31]. Biologically, the number of parameters prohibits a complete qualitative analysis of the model. We therefore simplify the task somewhat by assuming that dispersal is relatively large compared to landscape period. This assumption allows us to apply homogenization to Eq (2.1) and then follow the steps of deriving DDI conditions in a homogeneous landscape to determine instability conditions for the homogenized equations.

The homogenized model is given by

where $x\in \mathbb{R}$. The derivation follows Yurk and Cobbold ^[28] and is quite lengthy and technical; see also ^[19,29] for details and alternative formulations. In this model, the locally averaged population densities for prey and predator are denoted by $U(x, t)$ and $V(x, t)$, respectively. The homogenized diffusion coefficients for prey and predator are given by

and

Here, $r_i^u$ and $r_i^v$ are the residence indices ^[32] for prey and predator on patches of Type $i$, given by

and

In homogenized Eq (2.3), $\hat{f}(U, V)$ and $\hat{g}(U, V)$ are the homogenized reaction terms for prey and predator, whose expressions are

and

Here, $\langle r^u \rangle$ and $\langle r^v \rangle$ denote the average residence indices for prey and predator, given by

and

respectively.

We now follow the usual steps of a linear stability analysis for the homogenized Eq (2.3). We assume that the system has a spatially constant positive steady state, which we denote by $(U^*, V^*)$. We speak of diffusion-driven instability (DDI) if (i) the homogenized steady state is linearly stable with respect to spatially constant perturbations, and (ii) it is unstable to certain spatially varying perturbations; see ^[2,3]. The DDI conditions for the general homogenized Eq (2.3) are given as follows:

The first two of these conditions ensure stability with respect to homogeneous perturbations; the last two guarantee instability to certain spatial perturbations. In the Appendix, we give the expressions of these DDI conditions in terms of the patch-level quantities.

We had mentioned in the introduction that the spatial Rosenzweig–MacArthur predator–prey model cannot support Turing patterns. The reason is that the equation for the predator is linear in the predator density, so that the (2, 2) entry in the Jacobian matrix is zero and the required activator–inhibitor sign pattern is impossible (see Introduction). The following lemma generalizes this observation and shows that spatial heterogeneity alone cannot generate Turing patterns in a spatial Rosenzweig–MacArthur model of this type.

Lemma 1. If function $g_i(u_i, v_i)$ is linear in $v_i$ on both patches, then it is impossible to get Turing-pattern formation in the homogenized model.

Proof. Since $g_i(u_i, v_i)$ is linear in $v_i$ on both patches, in homogenized Eq (2.3), $\hat{g}(U, V)$ is linear in $V$. Hence, we can write $\hat g(U, V) = G(U)V$. At steady state then, we necessarily have $G(U^*) = 0$. The $(2, 2)$ entry of the Jacobian matrix at this steady state is $\partial_V \hat g(U^*, V^*) = G(U^*) = 0$. Hence, we do not have the required activator–inhibitor sign structure from Eq (1.1) in the model to support the existence of DDI.

Condition (DDI 4) is a quadratic equation in the ratio $\hat{D}^v/\hat{D}^u$. It can therefore be solved explicitly to get the condition that DDI can only occur if the ratio of the diffusion coefficients exceeds a critical ratio, namely,

At first sight, this expression is exactly the same as for the homogeneous problem. In particular, it appears that one can choose population dynamic parameters in such a way that the sign pattern requirement of the Jacobian matrix is satisfied and then choose diffusion constants so that the instability conditions are satisfied. However, the homogenized quantities on the right-hand side of Eq (2.13) already contain all model parameters, including the patch-level diffusion coefficients, as can be seen in the explicit expressions above. Hence, we cannot choose population dynamics and movement parameter separately. We will return to this issue later.

At the onset of DDI, when the inequality in Eq (2.13) is an equality, perturbations of a certain critical wave number will begin to grow. This critical wave number, $k_c^2$, is given by

and the corresponding critical wavelength, $w_c$, is

Since in our simulations, we set $\hat{D}_u = 1$, we have $\hat{D}_v = \hat{D}_c$ and can evaluate the critical wavelength with the above formula.

For our particular application, we use the May model as the reaction term on good patches; see Introduction and ^[11,33]. On bad patches, we replace the logistic growth of the prey by a simple linear death term. All other terms remain unchanged. In particular, predation still occurs in bad patches. Then the model is given by Eqs (2.1) and (2.2), where

Here, parameter $m$ denotes the prey's mortality rate on bad patches. As before, $b$ denotes the half-saturation constant of the Holling type II functional response ^[16,17], $s$ is the low-density growth rate of the predator, and $q$ is the proportionality factor between prey density and predator carrying capacity.

Remark 1. The May model is not defined when the prey density is zero. Since the positive steady state for the $v$-equation is $v = qu$, one can set $g(0, v) = 0$, and thereby allow $(0, 0)$ to be a biologically reasonable steady state: in the absence of prey, the predator cannot survive. However, we will only be interested in the positive coexistence steady state. Hence, our results will not be affected by singularity in the predator equation at zero.

If the landscape was homogeneous with reaction terms $f_1$ and $g_1$, we could choose parameter values so that the DDI conditions were satisfied and Turing patterns would form. If the landscape was homogeneous with reaction terms $f_2$ and $g_2$, the only spatially homogeneous steady state of the model would be the trivial state $(0, 0)$. Since DDI implies that the resulting pattern varies around the steady state and since the system preserves positivity, no DDI-induced patterns could possibly form there. In fact, for this case, any patterned steady states can be excluded since $f_2(u, v)\leq -mu$, so that $u\to 0$ irrespective of the value of $v$.

The homogenized reaction terms for the heterogeneous May model are explicitly given by

and

A spatially constant positive steady state, $(U^*, V^*) = (\langle r^u\rangle\widetilde U, \langle r^v\rangle\widetilde V)$, is given by

and

While the equation for $\widetilde V$ closely resembles the corresponding Eq (1.3), the equation for $\widetilde U$ is a cubic, except in the special case where $r_1^u = r_2^u$. The coefficients of the cubic depend on too many parameters to make meaningful statements about the roots in general. What can be said is that if $l_1 r_1^u-l_2 r_2^u m < 0$, then the sequence of coefficients of the cubic has no sign change and therefore there is no positive real root by Descartes' rule of signs. When the reverse inequality holds, then there is either one or three sign changes, which means that there exists at least one positive real root. The inequality is equivalent to the prey species not being able to grow at low density, see Eq (3.1). We solved for the roots numerically and checked that in all cases that we simulated there was at most one real root and hence only one positive steady state. We used this steady state to evaluate the DDI conditions. We simulate the reaction-diffusion system in Eq (2.3) to test the prediction of the DDI conditions. We give details on the numerical aspects in Section 6, where we also compare the solution of the homogenized and the nonhomogenized model.

Our model is too complex to determine all possible cases of DDI and pattern formation. Instead, we study several carefully chosen scenarios that illustrate several fundamental mechanisms by which landscape heterogeneity can enhance or reduce the likelihood of DDI in our predator–prey system.

Scenario 1 We choose parameters on patches of Type 1 such that all DDI conditions are satisfied on a homogeneous landscape of Type 1. As we mentioned above, no patterns can form on a homogeneous landscape of Type 2. We then expect that patterns emerge or disappear as we change the proportion of the two types of patches in the landscape.

Scenario 2 We choose parameters such that the DDI conditions are not satisfied on patches of Type 1. Since no pattern formation is possible on patches of Type 2, there is no obvious reason to believe that patterns could form on a landscape with both patch types mixed. However, we shall show that patterns are possible on some mixed landscapes. We consider two subcases, depending on which pattern formation conditions are violated on patches of Type 1.

(a) Here, we assume that the Jacobian matrix has an activator–inhibitor sign pattern but the ratio of the diffusion coefficients is below the critical ratio in patches of Type 1.

(b) Here, we assume that both diagonal entries of the Jacobian matrix are negative on patches of Type 1 so that there is no activator–inhibitor sign pattern.

In each scenario, we consider three different aspects of heterogeneity: (i) only population dynamics vary between patches, diffusion rates are constant in the landscape, and there is no preference for a particular type of patch; (ii) population dynamics and movement behaviour vary between patches but there is no preference for a particular type of patch; (iii) population dynamics vary between patches and there is preference for a patch type, but movement behaviour does not vary.

For the convenience of the reader, we summarize all parameters of the local and the homogenized model in Table 1.

3.   Analysis of Scenario 1

We choose population dynamics parameters such that on patches of Type 1 (good patches), we have an activator–inhibitor sign-pattern, whereas on patches of Type 2 (bad patches), we have no positive coexistence state. Then, we can choose diffusion coefficients in patches of Type 1, such that Turing patterns would form on an infinite homogeneous landscape ($l_2 = 0$) of this type. No Turing patterns are possible on an infinite homogeneous landscape of Type 2 ($l_2 = 1$). By continuity, we expect that patterns will form for small enough $l_2$ and will not form for large enough $l_2$ in homogenized Eq (2.3), where $\hat{f}(U, V)$ and $\hat{g}(U, V)$ are defined in Eqs (2.17) and (2.18), respectively. These values will depend on the other model parameters. We shall see that, in some cases, there is a unique threshold value of $l_2$ that separates the two cases, whereas in other cases, there are intermediate regions of pattern formation.

For this section, we choose the population dynamics parameters to be $b = 0.1$, $s = 0.2$ and $q = 0.8$, so that an activator–inhibitor sign pattern arises on patches of Type 1. On patches of Type 2, we choose $m = 0.6$. In the following section, we set the movement parameters equal in both patches; in later sections, we consider spatially varying movement and patch preferences.

3.1.   Patchiness in terms of population dynamics

In this section, only population dynamics but not movement behaviour change between patch types. For the population dynamics parameters in patches of Type 1, we calculate the critical diffusion ratio for pattern formation on an infinite landscape of this type to be $D_{c}\approx28.6$. We choose diffusion rates for prey and predator on patches of Type 1 such that their ratio exceeds this critical ratio. Since the movement behaviour does not change between patches, we choose the same values for diffusion coefficients in patches of Type 2. Specifically, for simulations we set $D^v_1 = D^v_2 = 40$ and $D^u_1 = D^u_2 = 1$. We assume that there is no patch preference for either species, so that $p^u_1 = p^v_1 = 0.5$.

A minimum requirement for pattern formation is that there will be a coexistence state. Explicitly evaluating the steady-state expressions in Eqs (2.19) and (2.20) is tedious if not impossible. A weaker, necessary condition is that the prey can persist in the system. The prey can persist if they can grow at low density. To evaluate this condition, we linearize the first equation of Eq (2.3) at the trivial equilibrium and derive conditions for which $\hat f_U$ is positive. We find that $\hat f_U > 0$ if and only if

This value serves as an upper bound for the threshold that we are looking for.

To find the actual range of values for $l_2$ where pattern formation is possible, we evaluate conditions DDI 3 and DDI 4 numerically as functions of $l_2$ in the range where DDI 1 and DDI 2 are satisfied. We plot the left-hand side of DDI 3 in Eq (2.12) as a function of $l_2$ (red curve in Figure 1). For this reason, DDI 3 is satisfied wherever the curve is positive. For DDI 4, we plot the critical value of Eq (2.13) of the homogenized diffusion coefficients, $\hat{D}_c$, as a function of $l_2$ (black solid in Figure 1). We also plot the actual ratio of the homogenized diffusion coefficients, $\frac{\hat{D}^v}{\hat{D}^u}$, as a function of $l_2$ (black dashed in Figure 1). Since the movement rates are identical between patches, this value is constant with respect to $l_2$. Then DDI 4 is satisfied when the black solid curve is below the black dashed line. The threshold value of $l_2$ for pattern formation is the smaller of the two values; in this case, it is the value where DDI 4 turns into an equality.

We check this result by simulating the homogenized model in Eq (2.3). We fix all biological parameters as above. According to Figure 1, pattern formation can occur for $l_2 < 0.43$. We choose $l_2 = 0.3$ for our simulations. With these parameter values, the steady state is given by $(U^*, V^*)\approx(0.1198, 0.0958)$ and the critical value of Eq (2.13) of the homogenized diffusion coefficient is $\hat{D}_c\approx 12$. Figure 2 shows a periodic spatial pattern as the steady-state profile of our system. From it, we calculate $w_c\approx 24$ as the critical wavelength.

Now we investigate how the characteristics of the patterns depend on the size of bad patches, $l_2$. More precisely, we compare how the critical wave number and wavelength that arise from DDI 4 (when all other conditions are satisfied) depend on $l_2$. Since we have chosen the diffusion coefficients independent of patch type and the patch preferences to equal 1/2 (i.e., there is no patch preference), neither the diffusion coefficients nor the patch preferences affect the averaged reaction terms $\hat{f}$ and $\hat{g}$. Also, the averaged diffusion coefficients equal the actual diffusion coefficients. Hence, the diffusion coefficients are independent parameters in the DDI conditions, just as they were in ^[2,3], and we can find the critical values. The situation will change in later sections.

We plot the critical wave number and wavelength as the red and blue curves, respectively, in Figure 3. For easier comparison, we include the critical values for a homogeneous landscape of Type 1 as the black dashed and solid lines. We observe that the critical wavelength is a concave up function of $l_2$. Compared to a homogeneous Type 1 landscape, the critical wavelength first decreases as the size of Type 2 patches increases. Then there is a minimal value, and when bad patches become even larger, the critical wavelength of the pattern increases. It appears that patterns are lost when the critical wavelength approaches infinity.

3.2.   Patchiness in terms of movement rates

In the previous section, movement behaviour was identical between the two patch types and there was no patch preference. Therefore, the homogenized diffusion coefficients were independent of the size of the two habitat types. There is ample evidence that individuals adjust their movement behaviour to their environment and do show preference for certain habitat types ^[34]. In this section, we explore how differences in movement rate affect pattern formation, in particular through Conditions DDI 3 and DDI 4. In the next section, we study the effect of patch preference. Similar considerations in a discrete-patch model showed that spatially varying movement rates can significantly affect parameter ranges for DDI ^[18].

For the population dynamics parameters in patches of Type 1, the critical diffusion ratio for pattern formation on an infinite landscape of this type is $D_{c}\approx 28.6$. We choose diffusion coefficients for both species in patches of Type 1 to be the same as Section 3.1. In our first scenario, we vary only the diffusion coefficient of species $v$ on patches of Type 2. We set $D^v_1 = 40$ and $D^u_1 = D^u_2 = 1$ as above and evaluate DDI conditions in the two-parameter plane of $D^v_2$ and $l_2$. Our second scenario is the same, except that we vary the diffusion coefficient of species $u$ on patches of Type 2. For both scenarios, we assume that there is no patch preference for either species, so that $p^u_1 = p^v_1 = 0.5$.

The two plots in Figure 4 show the regions in parameter space where the four DDI conditions are or are not satisfied. As in the previous section, there is a threshold value of $l_2$, above which the prey cannot persist. This threshold is indicated as the 'Extinction Boundary' (pink color). In the left plot, this threshold is a vertical line, because the critical value in Eq (3.1) does not depend on $D^v_2$. In the right plot, it is a curve, since the threshold does depend on $D^u_2$. To the right of the extinction boundary, no patterns can form, because the prey cannot even persist. To the left of the boundary, all pattern-formation conditions are satisfied in the red area. Patterns cease to exist because DDI 4 is not satisfied (blue area). In the white area, Conditions DDI 3 and 4 are both not satisfied, whereas in the green area, only DDI 3 is violated. Conditions DDI 1 and 2 are satisfied throughout the region where the prey can persist in the left plot, but there is a tiny region (orange) in the right plot, where DDI 1 is violated. Hence, there are (at least) two unstable modes in this case: the zero mode of a spatially constant perturbation with a pair of complex conjugate eigenvalues with positive real part, and a nonzero mode with positive real part for the usual diffusion-driven instability. We also see that when $D^v_2$ is large enough or small enough, there is a unique threshold value of $l_2$ that separates pattern formation from no pattern formation. For some intermediate values of $l_2$, however, there are two distinct intervals with respect to $l_2$ for which patterns can form.

We simulated the non-spatial system in the orange region of parameter space and found that, indeed, the coexistence state was unstable and temporally oscillating solutions emerged (Figure 5(a)). These solutions correspond to spatially constant, time-periodic solutions of the reaction-diffusion system in Eq (2.3), because of the no-flux boundary conditions. We also simulated the reaction-diffusion system for the same parameter values with the initial condition a small random perturbation of the spatially constant coexistence state. Those simulations show the emergence of spatial patterns (Figure 5(b)).

Finally, we take a closer look at the region where pattern formation occurs by including information about the critical wavelength, similar to Figure 3. In Figure 6 and similar figures to follow, the white curve delineates the region where all DDI conditions are satisfied (blue to yellow colours) from where they are not (uniform blue colour). Where the DDI conditions are satisfied, colours indicate the level sets of the critical wavelength with lighter colours indicating longer wavelengths (see side bar). The critical wavelength was computed as follows: for a given set of parameters, the homogeneous steady state was calculated and the Jacobian matrix was evaluated there. The critical wavelength is then computed from the critical wave number, where Condition DDI 4 becomes an equality. We emphasize that this value generally depends on the diffusion coefficients and the patch preferences through the averaging process. In particular, even if the Jacobian matrix has the required sign pattern, the choice of movement parameters may lead to DDI 4 not being satisfied. This is in contrast to the spatially homogeneous theory, where the Jacobian matrix is independent of the diffusion coefficients. In that case, as long as the Jacobian matrix has an appropriate sign pattern, one can always choose movement parameters so that Condition DDI 4 is satisfied.

We see from Figure 6 that the pattern in Figure 3 is quite common: along a horizontal transect, the critical wavelength is a concave up function of $l_2$: it decreases for small $l_2$ and increases for larger $l_2$. The largest wavelengths occur as one approaches the boundary where pattern formation is possible. For the left plot in Figure 6, we note that there is an interval for $D^v_2$ where pattern formation occurs for very small values of $l_2$ and again for intermediate values, but not in between. For the right plot in Figure 6, we note that the small region where the homogeneous steady state is unstable is shown with the white boundary.

From Figures 4 and 6, we see that it is possible to get patterns even for small values of $D^v_2$ and for larger values of $D^u_2$. As long as patches of Type 2 are small enough, the homogenized diffusion coefficient is still large enough so that the DDI conditions are satisfied.

3.3.   Patchiness in terms of patch preferences

In this section, we study the second aspect of movement behaviour in a patchy landscape: patch preference. We keep the diffusion coefficients constant in space. While it seems more reasonable to assume that organisms who can adjust their movement behaviour to local habitat conditions would adjust movement rate and patch preference, we consider each aspect separately to disentangle the different effects.

Given the population dynamics parameters in patches of Type 1, the critical diffusion ratio for pattern formation on an infinite landscape of this type remains the same as in the previous two sections, $D_{c}\approx 28.6$. We choose diffusion rates for prey and predator on patches of Type 1 such that their ratio exceeds this critical ratio. Since the movement behaviour does not change between patches, we choose the same values for the diffusion coefficients in both patches. Specifically, for simulations, we set $D^v_1 = D^v_2 = 40$, $D^u_1 = D^u_2 = 1$. In our first scenario, we vary only the patch preference of the prey species ($u$) and set $p^v_1 = 0.5$ as above. In the second scenario, we vary patch preference for the predator ($v$) but not for the prey ($p^u_1 = 0.5$). We evaluate the DDI conditions and illustrate our results in two-parameter planes of the preference of Type 1 patches of prey (see Figure 7(a) with $p_1^u = 1-p_2^u$) and predator (see Figure 7(b) with $p^v_1 = 1-p^v_2)$.

Several conclusions that can be drawn from this figure are similar to those in the preceding figure. For example, for a fixed value of $p^u_1$ or $p_1^v$, the critical wavelength is typically a concave-up function of $l_2$, and patterns disappear as the wavelength becomes large. At the boundary (white curve), condition DDI 4 becomes an equality and then fails to hold so that no patterns form. Also, pattern formation is possible for quite a wide range of parameter values. Interestingly, the pattern formation boundary is almost a horizontal line in plot (a): the critical value of $p_1^u$ is nearly constant for a range of $l_2$ between 0.1 and 0.5. In plot (b), on the other hand, the boundary is nearly vertical near $l_2\approx 0.45$ where $p^v_1$ can range from around 0.3 to around 0.7 with no visible change in the pattern formation boundary.

In the middle of the turquoise region of Figure 7(a), we see a small white region. This occurs as before where condition DDI 1 is violated and the non-spatial system has a periodic orbit. The situation here is the same as in the corresponding situation in the previous section.

4.   Analysis of Scenario 2a

In Scenario 1, the two patch types were chosen such that pattern formation occurred in a homogeneous landscape of Type 1 but not of Type 2. By shifting parameter $l_2$, which controls the relative size of the patch types, we could obtain patterns or have them disappear. In Scenario 2, we choose both patch types such that no patterns form on the corresponding homogeneous landscapes. Yet we will find that patterns can form for intermediate ratios of the two patch types in a heterogeneous landscape. We split this scenario into two cases:

(a) we assume that the Jacobian matrix on patches of Type 1 has the correct sign pattern (activator–inhibitor), but the ratio of the diffusion rates is too small to generate patterns;

(b) we assume that the sign pattern of the Jacobian matrix does not allow pattern formation on either patch type in isolation.

In this section, we choose parameters for population dynamics such that on patches of Type 1 (good patches), we have an activator–inhibitor sign-pattern and Conditions DDI 1 and DDI 2 are satisfied, whereas on patches of Type 2 (bad patches), we have no positive coexistence state. Then we calculate the critical diffusion ratio on patches of Type 1, where Conditions DDI 3 and DDI 4 are satisfied, and choose a ratio smaller than that so that no patterns can form on an infinite homogeneous Type 1 landscape. Obviously, pattern formation is impossible on a homogeneous Type 2 landscape. We conduct the same numerical experiments as in the preceding section.

4.1.   Patchiness in terms of population dynamics

We begin with the case where only population dynamics but not movement behaviour change between patch types. With the same population-dynamic parameters as in the preceding section, Conditions DDI 1 and DDI 2 are satisfied on Type 1 patches, and the critical diffusion ratio is $D_{c}\approx 28.6$ as before. We choose diffusion coefficients with a ratio smaller than that. Specifically, for simulations, we set $D^v_1 = D^v_2 = 15$ and $D^u_1 = D^u_2 = 1$, which gives a ratio of just over half of the critical ratio. We assume that there is no patch preference for either species, so that $p^u_1 = p^v_1 = 0.5$. As in Section 3, an upper bound for the threshold of $l_2$, above which the prey cannot persist, is given in Eq (3.1).

To see whether pattern formation is possible for the homogenized model, we evaluate Conditions DDI 3 and DDI 4 numerically as functions of $l_2$ in the range where DDI 1 and DDI 2 are satisfied. While condition DDI 3 is satisfied for all $l_2$ below some threshold (red curve in Figure 8), DDI 4 is satisfied only for intermediate values of $l_2$ but not as $l_2$ approaches zero (black curve in Figure 8). This behaviour, of course, reflects our choice of condition DDI 4 not being satisfied on a homogeneous Type 1 landscape. We compare and contrast this with Figure 1, where Condition DDI 4 holds for all small enough $l_2$ because Condition DDI 4 holds on a homogeneous Type 1 landscape there.

As before, we check whether and which patterns actually form by simulating homogenized system in Eq (2.3), where $\hat{f}(U, V)$ and $\hat{g}(U, V)$ are given by Eqs (2.17) and (2.18), respectively. We use the Crank-Nicolson scheme. We choose a specific value for $l_2$ in the range where the preceding analysis predicts pattern formation (here $l_2 = 0.1$), and we fix the other parameters as mentioned above. With these parameter values, the steady state is given by $(U^*, V^*)\approx(0.279, 0.223)$, and the critical value of Eq (2.13) for the homogenized diffusion coefficient is $\hat{D}_c\approx11.86$. Figure 9 shows a spatially periodic pattern as the steady-state profile of our system. From it, we calculate $w_c\approx21.84$ as the critical wavelength.

4.2.   Patchiness in terms of movement rates

We continue our investigation of Scenario 2a by choosing the movement rates to vary as we did in Scenario 1, Section 3.2. While keeping population dynamics parameters as in the first part of Scenario 2a, we now choose $D^v_1 = 15$ and $D^u_1 = D^u_2 = 1$ and evaluate DDI conditions in the two-parameter plane of $D^v_2$ and $l_2$; see Figure 10(a). We also vary the diffusion coefficient of species $u$ on patches of Type 2 see Figure 10(b). For both scenarios, we assume that there is no patch preference for either species, so that $p^u_1 = p^v_1 = 0.5$.

In Figure 10, the white curve delineates the region where all DDI conditions are satisfied (blue to yellow colours) from where they are not all satisfied (uniform blue colour). Where the DDI conditions are satisfied, the colours indicate the level sets of the wavelength of emerging patterns with lighter colours indicating longer wavelengths (see side bar). As in the previous section, there is a threshold of $l_2$, the 'Extinction Boundary' (pink), above which the prey cannot persist.

Figure 10 shares many similarities with Figure 6 but also shows some crucial differences. The most important is that pattern formation is impossible for very small and large values of $l_2$ by our choice of Scenario 2a. Nonetheless, we see that pattern formation is possible for a large range of intermediate parameter values. Patterns are more likely to form when $D^v_2$ is large or $D^u_2$ is small, since then the ratio of the homogenized diffusion coefficients is larger, and we know that this ratio typically needs to be large for patterns to form.

Similar to Figure 6, Figure 10(a) also shows that for fixed values of $D^v_2$, the critical wavelength is a concave up function of $l_2$. Patterns cease to exist when the wavelength approaches infinity. The situation is similar in Figure 10(b) where the critical wavelength is concave up with respect to $l_2$ for any fixed $D^u_2$. In that Figure, we also see again a small spot where Condition DDI 1 is violated in the turquoise region. As in the scenarios before, there is a periodic orbit in the non-spatial system.

4.3.   Patchiness in terms of patch preferences

We conclude our investigation of Scenario 2a by varying the patch preferences of the two species while keeping other parameters fixed, similar to what we did in Scenario 1, Section 3.3. Parameters other than patch preferences are set exactly as in Section 4.1. As in Figure 7, we vary the patch preference of the predator and prey species separately and evaluate the corresponding DDI conditions. We present and illustrate our results as contours of critical wave lengths in the corresponding two-parameter planes of $p^u_1 = 1-p^u_2$ or $p^v_1 = 1-p^v_2$ and $l_2$; see Figure 11.

The plots in Figure 11 should be compared to those in Figure 7. They show a number of similarities, but their biggest difference is that pattern formation is not possible for $l_2 = 0$ by the setting of Scenario 2. Hence, there is a region without pattern formation near $l_2 = 0$. Interestingly, when $p^u_1$ is very small, pattern formation is possible for very small values of $l_2$. There seems to be a trade-off: when $l_2$ is small, bad patches are small, but when the preference for these bad patches is high, their effect appears larger. As before, for fixed values of patch preferences, the critical wavelength is a concave-up function of $l_2$.

5.   Analysis of Scenario 2b

In contrast to the last two sections, we now choose population dynamics parameters such that on patches of Type 1 (good patches), we do not have an activator–inhibitor sign-pattern. It is therefore impossible to get Turing patterns on a homogeneous Type 1 landscape. More precisely, by choosing parameter $q$ to be less than $q_A$, the (1, 1) entry in the Jacobian matrix is negative. This means that Condition DDI 3 cannot be satisfied for any values of the diffusion coefficients, even though Conditions DDI 1 and DDI 2 can be. Patches of Type 2 (bad patches) remain the same as before, where we have no positive coexistence state, so that no Turing patterns are possible on an infinite homogeneous landscape on Type 2.

For numerical simulations, we choose the population dynamics parameters $b = 0.1$, $s = 0.2$ and $q = 0.5 < q_A\approx 0.672$. This ensures that the Jacobian matrix on Type 1 patches does not have an activator–inhibitor pattern. On patches of Type 2, we choose $m = 0.6$ as before. As in the previous scenarios, we first set the movement parameters equal in both patches; later, we consider spatially varying movement and patch preferences.

5.1.   Patchiness in terms of population dynamics

As in the previous scenarios, we begin our investigation of Scenario 2b with the case where only population dynamics but not movement behaviour changes between patch types. We set $D^v_1 = D^v_2 = 40$, and $D^u_1 = D^u_2 = 1$. We assume that there is no patch preference for either species, so that $p^u_1 = p^v_2 = 0.5$. As before, Eq (3.1) defines a threshold of $l_2$, above which the prey cannot persist.

We begin with the analogue of Figures 1 and 8; i.e., we visualize Conditions DDI 3 and DDI 4 as a function of $l_2$. We checked that Conditions DDI 1 and DDI 2 are satisfied in the range $l_2\in[0, 0.6]$. The plot in Figure 12 shows that Condition DDI 4 is satisfied for $l_2$ near zero, but DDI 3 is not. This is in contrast to Scenario 1, where both were satisfied near zero, and Scenario 2a, where DDI 3 was satisfied near zero but DDI 4 was not. As $l_2$ is increased, DDI 4 is violated before DDI 3 holds. Then there is an intermediate range of $l_2$ where both conditions hold (approximately for $l_2\in[0.29, 0.44]$) before both are violated for larger values of $l_2$.

We can explain the emergence of an activator–inhibitor sign pattern for intermediate values of $l_2$ mathematically and biologically. Mathematically, we note that the prey density at the coexistence state decreases as $l_2$ increases. Furthermore, the derivative of $\hat f$ with respect to $U$, the (1, 1)-entry of the Jacobian matrix, is a decreasing function of $U$. Hence, there can be a range of $l_2$ where $U$ is small enough to make this entry positive and give the required sign pattern. Biologically speaking, the productivity of a population (its rate of change) is typically the highest at intermediate population densities. At low density, there are not enough individuals around; at high density, there are too many. When there are bad patches into which some individuals are moving, the steady state density is decreased, and productivity is increased. If productivity is high enough, we have an activator.

As in the last two scenarios, we illustrate pattern formation by simulating homogenized Eq (2.3), using the Crank-Nicolson scheme. We fix all parameters as above and choose $l_2 = 0.4$. Then the steady state is given by $(U^*, V^*)\approx(0.1296, 0.0648)$, and the critical value of Eq (2.13) of the homogenized diffusion coefficient is $\hat{D}_c\approx 28.84$. Figure 13 shows a periodic spatial pattern as the steady-state profile of our system. From it, we calculate $w_c\approx32.53$ as the critical wavelength.

5.2.   Patchiness in terms of movement rates

In the previous section, movement behaviour was identical between the two patch types, so that the homogenized diffusion coefficients were independent of the sizes of the two habitat types. Similar to Sections 3.2 and 4.2, we now explore how differences in movement rate affect pattern formation, in particular through Conditions DDI 3 and DDI 4. As before, we vary $D^v_2$ and $D^u_2$ as well as $l_2$ and visualize the critical wavelength in the corresponding parameter planes. Unless otherwise noted, all parameters are as in Section 5.1.

In Figure 14, the colour coding, the extinction boundary and the white curve that delineates the regions where pattern formation is possible are exactly the same as in previous figures, such as Figure 3. No pattern formation is possible near $l_2 = 0$. Figure 14(a) shows the intermediate range of $l_2$ where pattern formation is possible provided $D^v_2$ is large enough. Similarly, Figure 14(b) shows the intermediate range of $l_2$ where patterns can form when $D^u_2$ is small enough. We note that pattern formation is possible very close to the extinction boundary here. As we have seen in all previous plots of this kind, the largest wavelengths occur as one approaches the boundary where pattern formation is possible.

5.3.   Patchiness in terms of patch preferences

We conclude our investigation of Scenario 2b by allowing patch preferences to vary while keeping movement rates constant in space, as we did in Sections 3.3 and 4.3. Population dynamics parameters and diffusion coefficients are as in Section 5.1. We evaluate the DDI conditions and illustrate our results in two-parameter planes of the preference of Type 1 patches of prey (see Figure 15(a) with $p_1^u = 1-p_2^u$) and predator (see Figure 15(b) with $p^v_1 = 1-p^v_2$).

Several conclusions that can be drawn from this figure are similar to those in the corresponding preceding figures. For example, at the boundary (white curve), Condition DDI 4 becomes an equality and then fails to hold so that no patterns form. Pattern formation is not possible for quite a wide range of parameter values, especially in plot (b). When the preference of the prey for patches of Type 1 is large, no patterns can form. When the preference is low, in particular when prey prefer patches of Type 2, patterns can form for intermediate values of $l_2$. The biological explanation behind this mathematical observation is again in the question of release from intraspecific competition. When prey have a preference for Type 1 patches, not many leave for Type 2 patches; hence, there is no relief from this competition. When the preference is for Type 2 patches and those patches are large enough, the release from competition is strong enough so that an activator–inhibitor sign pattern exists and patterns can form. If, however, the Type 2 patches are too long, the overall population density will be too small for pattern formation.

In plot (a), as the preference of prey for patches of Type 1 ($p^u_1$) increases, we get Turing patterns for the larger values and a wider range of $l_2$. Moreover, in this plot, we can see that when $p^u_1$ exceeds 0.5, $l_2$ needs to be large enough for patterns to form. The reason is again that a large enough bad patch will provide release from intraspecific competition in patches of Type 1. Finally, in contrast to plots (a) in Figures 7 and 11, in this plot, we do not have any hole in the region where we get patterns. In plot (b), we observe that it is possible to get patterns when $p^v_1$ varies approximately between 0.3 and 0.7, but only when Type 2 patches are long enough.

6.   Numerical aspects

In this section, we describe in more detail how we solved the homogenized and the nonhomogenized reaction-diffusion equations and we give some examples that the former approximates the latter quite well.

To simulate the homogeneous or homogenized model, we chose two independent implementations of numerical reaction-diffusion solvers. One was $\mathsf{Matlab's}$ built-in $\mathtt{pdepe}$ solver, which we used with the standard error tolerances. The other was the finite-difference Crank-Nicolson scheme ^[35], which we implemented also in $\mathsf{Matlab}$. We used the same spatial and temporal discretization as well as the final time for both schemes, namely $\mathrm{d}x = 0.05$, $\mathrm{d}t = 0.1$ and $T = 3000$ unless otherwise noted. With these parameters, we found the simulations to be stable and solutions of the time-dependent problem to reach the time-independent steady state. The simulations of the two different methods agreed. To approximate the infinite spatial domain, on which the analysis rests, we chose a relatively large bounded numerical domain and implemented no-flux (Neumann) boundary conditions. The formulas obtained analytically guided our choice of domain size in that we ensured that the numerical domain was large enough to contain several periods of the expected pattern.

Finally, we implemented a fractional step method Strang-splitting to solve the nonhomogenized reaction-diffusion system. To iterate each time step, we used half a time step of diffusion, implemented as a Crank-Nicholson scheme (see above), followed by a full time step of reaction using a fourth order Runge–Kutta scheme, followed by half a time step of diffusion, implemented as before. To ensure that the whole scheme remains second-order accurate, we used second-order forward or backward difference to approximate the derivatives in the interface conditions. Since each patch is very small, we used a smaller spatial discretization than in the homogeneous case, which then also required a smaller time step ($\mathrm{d}x = 0.01, \mathrm{d}t = 8 \times 10^{-4}$). Since the pattern was already well established at $T = 1000$, we ended the simulations there. Numerical simulations for a single population have already demonstrated that the homogenization provides a very good approximation to the nonhomogenized model, e.g., Figure 2 in ^[29]. We note that the requirement of very small time steps in the nonhomogenized model leads to a much greater demand on CPU time for simulation than in the homogenized model at roughly equal accuracy. In this case, the time requirement was about $10^5$ times that of the homogenized model.

We show one representative simulation outcome in Figure 16. This figure is based on Scenario 2a. As we had done there, we chose the predator to have no patch preference and equal diffusion rates in both patches. Consequently, the predator density is continuous at all interfaces. Hence the simulation result from the nonhomogenized model (black) and the homogenized model (red) closely match (Figure 16(a)). We chose the prey to have preference for patch type 2, so that the homogenized model predicts pattern formation (see Figure 7(a)). Now the prey density is discontinuous at each interface, and the homogenized model is expected to approximate some average spatial density in the nonhomogenized model. This is indeed the case. The solid red curve in Figure 16(b) shows the solution of the homogenized model, whereas the black curve shows the solution of the nonhomogenized model. The two dashed red curves are multiples of the solid red curve, given by $u*r_1^u/\langle r^u\rangle$ (bottom curve) and $u*r_2^u/\langle r^u\rangle$ (top curve). Hence, these two dashed curves are what one obtains when recovering the nonhomogenized solution from the homogenized one. They show excellent agreement.

7.   Conclusions and discussions

Uncovering mechanisms that lead to observed spatial patterns in ecological systems is still a fundamental challenge in ecology. The pioneering work on pattern formation was conducted by Alan Turing ^[1], who discovered pattern formation through DDI by deriving conditions under which a spatially homogeneous steady state that is stable in the absence of diffusion becomes unstable in the presence of diffusion. This mechanism determines the spatial pattern that evolves ^[2,3]. Although Turing's study was originally proposed in chemistry, his idea spread to biology; e.g., developmental biology and animal-coat patterns ^[2]. Turing patterns are generated by the interaction of two substances, an activator and an inhibitor, with different diffusion rates: the inhibitor must diffuse faster than the activator. We specialize the general Turing pattern idea to ecological population dynamics and predator–prey interactions.

A number of authors have explored the application of DDI to population dynamics in ecology ^[7,9,15], although these have been largely confined to studying pattern formation in a homogeneous landscape. In this paper, we examined the possibility for Turing-pattern formation in a heterogenous landscape. Since we assumed that the period of the underlying landscape heterogeniety is small, we could study pattern formation in a patchy landscape by applying the homogenization technique to derive a homogenized model. From this homogenized model we could apply the standard approach ^[2,3] to derive DDI conditions for this simplified model. Our main results (Figures 6, 10 and 15) show that even when pattern formation is impossible in a homogeneous landscape, heterogeneity could generate pattern formation. By varying the relative size ($l_2$) of the bad patches in the landscape we found that we could control whether we obtained patterns or, by contrast, have them disappear (Figures 4 and 7).

Two requirements of diffusion driven instability have limited the application of Turing patten formation in ecology, one being the ecological requirement of an Allee effect in prey growth rate or self-damping in the predator, the second requirement being that the predator diffuses significantly further than the prey. Together, these are quite restrictive and not necessarily common place in predator–prey systems. Much work has been conducted into finding ways to relax these requirements so that Turing's theory can be applied to an ecological setting, and here we argue that heterogeneity can provide at least a partial solution to this problem. In Scenario 2a we relaxed the requirement of high predator diffusion (Section 4). While we still required the predator to diffuse more rapidly than the prey, we found we could allow these rates of diffusion to be much closer in value and still obtain pattern formation. Preference for a particular patch type and the relative sizes of patches in which movement is fast or slow act to speed up or slow down movement at the landscape scale, thereby enabling the system to satisfy the Turing pattern formation conditions at the level of the landscape (homogenized model), but not requiring these conditions to be met locally (patch model). Finally, in Scenario 2b, we showed how the need for an Allee effect or self-damping could be relaxed. We found spatial heterogeneity could turn a system without an activator–inhibitor sign structure at the fine spatial scale into a system that has one at landscape spatial scale, enabling Turing patterns to emerge.

Our study is closest to the one by Cobbold et al. ^[18], who also consider pattern formation on a heterogeneous landscape composed of patches. They also find that spatial heterogeneity can facilitate pattern formation. They assume that movement within a patch is fast, but movement between patches is localized and confined to neighbouring patches only, which is in contrast to our assumption, that movement is rapid and individuals encounter many patches over their lifetime. Despite this difference in modelling assumptions, several of our results are similar to those of Cobbold et al. For example, they found that no patterns form when prey movement was heavily biased toward good patches. Similar results can be found in our work, as seen in Figures 11(a) and 15(a). Cobbold et al. also compared the range of pattern formation in a homogeneous landscape with a heterogeneous landscape. They showed that the range of parameter $q$ that leads to pattern formation in a homogeneous landscape is much smaller than in a heterogeneous landscape. We reached a similar conclusion, in Scenario 2b when we compared the effects of parameter $q$ in the two different landscapes. The consistency of these results across both fine grain (our work) and coarse grain (Cobbold et al.) spatial heterogeneity highlights the robustness of many of our findings.

Further studies that have considered pattern formation on heterogeneous domains have focused on only the diffusion coefficients depending on space ^[20,21], which has the benefit that there is a spatially constant steady state around which one can linearize. These studies chose diffusion coefficients such that pattern formation is possible on one side of a two-patch domain but not on the other. Under appropriate conditions, they found that patterns can emerge on the entire domain through DDI. In the models by Page et al. ^[22,23], only the population dynamics change in space but not the diffusion terms. They showed that, under some conditions, the mismatch of the steady state at the interface can generate patterns that are centred around this discontinuity and then decays away from it to the boundaries of the domain. The patterns found in a similar set-up by Kozák ^[24] extend on either or both sides of the discontinuity but are still of very small wavelength compared to the patch size. These spatially localized patterns are in contrast to the patterns we find, which propagate across the entire domain and are generated by Turing instability rather than the discontinuity at the patch interface. Underlying the work of Benson et al., Page et al. and Kozák et al. is the assumption that the length of patches is larger than the length of the patterns and patterns emerge as variation in the densities in one patch. Instead, we have assumed that the pattern wavelength is much larger than the length of patches. There is currently no theory in place to study the case when the pattern wavelength and patch lengths are comparable.

Homogenization has provided a useful method for studying Turing instabilities in spatially heterogeneous systems with rapidly varying heterogeniety. Krause et al. ^[36] recently used WKBJ asymptotics to derive a local version of the classical Turing conditions which allows one to formally study the case of slowly varying spatial heterogeneiety. In contrast to our study, Krause et al. find that patterns are localised in space, similar to Page et al. Also in contrast to our study, the local Turing conditions require that the classical Turing conditions hold locally, which is not required in the landscapes we consider. This raises the question: at what scale of heterogeneity do our pattern formation conditions break down? We know from the work of Page, Benson and Krause (see above) that large scale heterogeneity leads to localised patterns and these occur in a more restricted range of parameter space, but what happens at intermediate scales of heterogeneity still remains an open question.

Preliminary work suggests that our findings are likely to be robust to the choice of predator–prey model. We chose the May model on good patches and the no-growth model for prey on bad patches. If we replace the May model by the Beddington–DeAngelis functional response, we also observe that it is possible to get Turing patterns for some range of the size of bad patches (plots not shown). Hence, an interesting task is to replace the May model on Type 1 patches with different functional responses such as the ratio-dependent model ^[10] and the modified Rosenzweig–MacArthur model that was introduced by Fasani and Rinaldi ^[15]. Then we can study how heterogeneity may affect Turing patterns for each of these models in terms of how the characteristics of the patterns such as the wavelength etc. depend on the size of bad patches, $l_2$. Another interesting question is to explore spatial patterns when we no longer have strict bad and good types of patches; e.g., a heterogeneous landscape consisting of two different type of patches, with the same functional form of the model, say the May model, in both patches but different values for the model parameters. Arrangement of patches is another modelling assumptions that should be tested for robustness.

Our work here has relied on the technique of homogenisation which requires the size of the good and bad patches to be very small. Yurk and Cobbold have shown that homogenization can also give a good approximation for heterogeneous reaction-diffusion equations even for moderate size patches ^[28]. Moreover the theory of homogenization has also been extended to two dimensional systems ^[37] and to advective systems ^[38] which paves the way for our approach to be extended to more complex settings in the future.

Acknowledgements

CC was supported by a Leverhulme Research Fellowship RF-2018- 577\textbackslash 9. FL gratefully acknowledges funding from the Natural Sciences and Engineering Research Council of Canada under a Discovery Grant (RGPIN-2016-04795) and a Discovery Accelerator Supplement (RGPAS-492878-2016).

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix

We express the DDI conditions for the homogenized model in Eq (2.3) in terms of the patchy growth functions, $f_i$ and $g_i$, in order to evaluate them at the steady state in the patch-level model. It may sometimes be easier to evaluate the DDI conditions on the patchy level rather than on the homogenized level.

We derive these conditions on the patch-level by applying the chain rule for each of the following partial derivatives to obtain

and

We have similar expressions for partial derivatives of function $g_i$.

As a result, Conditions $\mathrm{DDI}\; 1$–$\mathrm{DDI}\; 4$ in terms of $f_i$ and $g_i$ are given by

where all these expressions are evaluated at the steady state $(U^*\frac{r^u_i}{\langle r^u \rangle}, V^*\frac{r^v_i}{\langle r^v \rangle})$.