<figure><img src="image_1.png" /></figure>

HAL
open science

Dynamics of a Diffusive Two-Prey-One-Predator Model
with Nonlocal Intra-Specific Competition for Both the
Prey Species

Kalyan Manna, Vitaly Volpert, Malay Banerjee

► To cite this version:

Kalyan Manna, Vitaly Volpert, Malay Banerjee. Dynamics of a Diffusive Two-Prey-One-Predator Model with Nonlocal Intra-Specific Competition for Both the Prey Species. Mathematics , 2020, 8 (1), pp.1-28. 10.3390/math8010101 . hal-03027297

HAL Id: hal-03027297

https://hal.science/hal-03027297v1

Submitted on 12 Mar 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Article

# Dynamics of a Diffusive Two-Prey-One-Predator Model with Nonlocal Intra-Specific Competition for Both the Prey Species

Kalyan Manna<fnref n="1" />, Vitaly Volpert<fnref n="2" />,<fnref n="3" />,<fnref n="4" /> and Malay Banerjee <fnref n="1" />,*

Received: 13 December 2019; Accepted: 2 January 2020; Published: 7 January 2020

<figure><img src="image_4.png" /></figure>

check for
updates

**Abstract:** Investigation of interacting populations is an active area of research, and various modeling approaches have been adopted to describe their dynamics. Mathematical models of such interactions using differential equations are capable to mimic the stationary and oscillating (regular or irregular) population distributions. Recently, some researchers have paid their attention to explain the consequences of transient dynamics of population density (especially the long transients) and able to capture such behaviors with simple models. Existence of multiple stationary patches and settlement to a stable distribution after a long quasi-stable transient dynamics can be explained by spatiotemporal models with nonlocal interaction terms. However, the studies of such interesting phenomena for three interacting species are not abundant in literature. Motivated by these facts here we have considered a three species prey-predator model where the predator is generalist in nature as it survives on two prey species. Nonlocalities are introduced in the intra-specific competition terms for the two prey species in order to model the accessibility of nearby resources. Using linear analysis, we have derived the Turing instability conditions for both the spatiotemporal models with and without nonlocal interactions. Validation of such conditions indicates the possibility of existence of stationary spatially heterogeneous distributions for all the three species. Existence of long transient dynamics has been presented under certain parametric domain. Exhaustive numerical simulations reveal various scenarios of stabilization of population distribution due to the presence of nonlocal intra-specific competition for the two prey species. Chaotic oscillation exhibited by the temporal model is significantly suppressed when the populations are allowed to move over their habitat and prey species can access the nearby resources.

**Keywords:** prey-predator; diffusion; nonlocal interaction; Turing instability; spatiotemporal pattern

## 1. Introduction

Generally, the ecological systems and interactions among different species in nature are too complex to be understood with ease. In order to comprehend and make predictions about the dynamics of a complex ecological system, mathematical models and quantitative mathematical methods act as useful tools. Since Paine's discovery on predation-mediated species diversity [1], several theoretical studies have been conducted in order to understand the role played by predation on the dynamics of temporal food web models (especially on two-prey-one-predator temporal models).

<footnote n="1">¹ Department of Mathematics and Statistics, Indian Institute of Technology Kanpur, Kanpur 208016, India; kalyan274667@gmail.com</footnote>

<footnote n="2">² Institut Camille Jordan, UMR 5208 CNRS, University Lyon 1, 69622 Villeurbanne, France; volpert@math.univ-lyon1.fr</footnote>

<footnote n="3">³ Peoples Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya St, Moscow 117198, Russia</footnote>

<footnote n="4">⁴ INRIA Team Dracula, INRIA Lyon La Doua, 69603 Villeurbanne, France</footnote>

<footnote>* Correspondence: malayb@iitk.ac.in</footnote>
In this direction, Fujii [2] considered a two-prey-one-predator model with the competition between the two species of prey proposed in [3] and investigated for predation-mediated stabilization of the system. He showed the existence of a globally stable limit cycle around the unstable coexistence equilibrium point, whereas the corresponding two-species system without predator is unstable. On the other hand, Vance [4] argued that this coexistence of three species is not a result of predation alone but it occurs due to the combined effect of resource partitioning for the two prey species and predation. Working on the same model with Lotka-Volterra type linear functional responses, Takeuchi and Adachi [5] showed that chaotic coexistence is also possible apart from the already established periodic one. It was demonstrated that chaotic coexistence bifurcates from the periodic one when one prey's competitive ability is greater than the other. A criterion for permanent coexistence was established in [6] for the two-prey-one-predator model with intra-specific competition for predator.

Moreover, Matsuda et al. [7] investigated the effect of switching property of predator on the stability of the three-species model and concluded that heterogeneity in intrinsic growth rates for prey is an important component in order to achieve the stabilization. They further showed that the tendency of switching response toward the prey with lesser intrinsic growth rate enhances the stabilizing effect. Abrams and Matsuda [8] investigated the consequences of the adaptive anti-predator behavior of prey in a two-prey-one-predator system. In [9], the authors showed that the occurrence of chaotic dynamics for two-prey-one-predator food chain is independent of the particular choice of the model. Abrams showed that the amplitude of population cycles increases through enrichment in a model with saturating functional response and this enlargement of cycles favors the endangered species in [10]. In [11], the authors defined a functional response taking into account both the prey densities and showed that the switching effect of predator can potentially cause the extinction of predator species. Lee and Kajiwara [12] proposed a two-prey-one-predator model where one prey species consumes the remaining part of carcass of other prey species left by predator. They showed that this type of consumption can eventually lead to chaos. In case of non-interacting preys, some sufficient conditions for the stable coexistence equilibrium and periodic solution were provided in [13]. Recently, Groll et al. [14] examined a three-species microbial model and found the existence of chaotic regime between two regimes of stable limit cycles. They associated the reason for this type of chaotic regime to the competition between two distinct limit cycles.

Prey-predator patterns for ecological communities are ubiquitous in natural habitat such as the patchy spatial distribution in plankton ecosystem and the chaotic distribution for vole populations in northern Fennoscandia [15,16]. Generally, the spatial component acts as a crucial ingredient in shaping the spatial structure of an ecological community. Over the past few decades, many researchers have paid attention in order to understand the emergence of spatiotemporal patterns in ecological communities and their possible mechanisms since the seminal work of Turing [17]. Theoretical investigations on spatiotemporal pattern formations are generally accomplished by analyzing the associated reaction-diffusion equations. Till date, numerous works on the pattern formation for two-species models are present in literature, and some of these works can be found in [18–20] and the references therein. As a result, a wide variety of stationary patterns such as hot spots, cold spots, labyrinthine and mixture of spots and stripes, and dynamic patterns such as traveling waves, periodic traveling waves and spatiotemporal chaos etc. have been reported [18–20]. Also, Turing and Hopf-Turing bifurcations have been identified as two common mechanisms for the emergence of these different types of patterns [18–20].

However, the number of works on pattern formations for three-species food chain models is somewhat limited as compared to the two-species cases. For instance, White and Gilligan provided sufficient conditions to induce diffusion-driven instability for a generic three-species reaction-diffusion system in [21]. Petrovskii et al. [22] considered a three-species competitive spatial model of Lotka-Volterra type, and showed the existence of traveling waves and spatiotemporal chaos. Chen and Peng [23] proved the existence of stationary patterns for a competitor-competitor-mutualist model with cross-diffusion by establishing the existence of non-constant positive steady states. Also,
they demonstrated that the corresponding model without cross-diffusion fails to produce any stationary pattern. In [24], the authors examined a spatially extended one-prey-two-predator system with defense switching mechanism for prey and collaborative exploitation of prey by predators. Their study also held cross-diffusion responsible for the emergence of stationary patterns. On the other hand, Wang [25] found stationary patterns for three-species models with both the prey-dependent and ratio-dependent functional responses as a result of self-diffusion. The author further showed the existence of stationary patterns for these models in presence of the cross-diffusion as well in [26]. Tian [27] presented both the theoretical and numerical results on cross-diffusion induced Turing patterns for Holling type II and Leslie-Gower type three-species food chain model. Similar types of results were presented for a predator-prey-mutualist system in [28]. Rao [29] investigated the complex spatiotemporal dynamics for an one-prey-two-predator system with ratio-dependent functional responses and showed the transition of different types of stationary patterns to chaotic wave patterns depending on the magnitudes of the maximum ingestion rate or the mortality rate of intermediate predator. In [30], the authors examined for the possible pattern formations in a three-species food chain model with the Holling type II and a modified Leslie-Gower type functional responses over a circular domain. Some other studies in this direction can be found in [31–36]. Very recently, Mukherjee et al. [37] employed weakly nonlinear analysis for a three-species model in order to obtain amplitude equations which determine the thresholds for emergence and stability of Turing patterns in the vicinity of the Turing bifurcation boundary.

The majority of the existing works on pattern formation in ecological communities have been accomplished by considering the local intra- and inter-specific interactions, that is, an individual located at a specific spatial point $\hat{x}$ can only interact with individuals present at that point. However, spatial mobility of a species gives rise to the scenario that the resource consumptions at the nearby locations and hence interactions should depend on a weighted average of the population density in a neighborhood of $\hat{x}$ instead of the pointwise values of the population density [38,39]. Nonlocal interaction can be modeled by incorporating a convolution integral with a specified non-negative kernel function which takes care of the weighting aspect [38,39]. The incorporation of nonlocality into a model results in a system of integro-partial differential equations. In the available literature, several different types of kernel functions have been considered to model nonlocal interactions. These different types of kernel functions include top-hat kernel, Gaussian kernel, Laplacian kernel, triangular kernel, parabolic kernel, etc. [40–42]. Among these kernel functions, the simplest one is the symmetric top-hat kernel which is basically a step function, and incorporation of it in the modeling approach can lead to very rich and complex spatiotemporal pattern formations [40,42].

Over the last decade, considerable amount of efforts have been made in order to understand the effects of nonlocal interactions on the spatiotemporal dynamics possessed by ecological communities. For instance, Gourley considered a nonlocal Fisher equation and established the existence of traveling wavefront solutions connecting two spatially homogeneous steady states for sufficiently weak nonlocality in [43]. Merchant and Nagata showed the destabilizing effect of nonlocal prey competition and reported various complex spatiotemporal patterns such as stationary periodic in space patterns, wave trains, and irregular spatiotemporal oscillations in [41]. In a subsequent study [44], they studied the correlation between the properties of selected wave trains and the standard deviation of nonlocality, and suggested that highly nonlocal systems are more likely to produce spatiotemporal chaos. Segal et al. [42] investigated pattern formations in a nonlocal model of competing populations and showed that the population distribution can eventually evolve into localized island structure. Further, they described analytically the structure of these islands by using step-function kernel. In [38], the authors studied a model for two competing species with asymmetric nonlocal coupling and found the existence of different types of spatiotemporal solutions such as propagating traveling waves of islands, colony formation, and modulated traveling waves. Deng and Wu [45] analyzed the global stability property of a nonlocal single species population model. In [46], the authors investigated the
dynamics of a prey–predator model with nonlocal consumption of prey, and showed the existence of multiple stationary states, simple and periodic traveling waves. They further showed that the incorporation of nonlocal interactions in the classical Rosenzweig–MacArthur model can produce Turing patterns, while the local counterpart of it is unable to do so [47]. Some other interesting works on the impacts of nonlocal interactions on spatiotemporal dynamics for two-species models can be found in [40,48–54].

On the other hand, there exist very few studies on the effect of nonlocal interactions for three-species food chain models. Recently, Autry et al. [39] considered a three-species food chain model with ratio-dependent functional response and two types of nonlocality. They addressed the issue of biological pest control and showed that biological control cannot be achieved for highly diffusive pest species, however, robust partial control can be achieved if the super-predator is sufficiently diffusive and behaves nonlocally. In [55], the authors investigated the spatiotemporal pattern formation in a nonlocal intraguild predation model where the nonlocal interaction is incorporated in the growth of the shared resource. They identified the transition of stationary Turing patterns to spatiotemporal chaotic patterns through dynamic oscillatory patterns depending upon the magnitudes of the extent of nonlocality.

In this article, our aim is to provide rich spatiotemporal dynamics possessed by models of three interacting species (in particular, two preys and their common predator). Mainly, we are interested to examine the effects of random dispersal of all the three species and nonlocal intra-specific competition for two prey species on the population distributions. Apart from these, we are also interested to investigate some of the key issues in ecology such as long transient phenomenon and habitat size dependence of resulting patterns. In this regard, we investigate systematically a temporal model, a reaction-diffusion model and a nonlocal model. These systematic investigations can be useful in comparing the resulting dynamics and their appropriate causes as these three models constitute a system of nested models. The rest of this article is organized in the following fashion. In Section 2, we introduce a three-species (two-prey-one-predator) temporal model and briefly discuss about the dynamical behaviors. Then in Section 3, we extend the model by incorporating spatial components and present a wide variety of spatiotemporal dynamics possessed by it. The spatiotemporal model is further extended by considering the nonlocal intra-specific competitions for both the prey species in Section 4 with the manifestation of possible impacts of the extent of nonlocality on the spatiotemporal dynamics. Finally, we end this article with a brief discussion in Section 5.

## 2. Temporal Model

In this section, we introduce a three-species (two-prey-one-predator) temporal model for prey–predator interactions presented in [5] which will serve as a base model for our study. The concerned temporal model is given by the following system of three ordinary differential equations:

$$ \frac{du_1}{dt} = u_1(b_1 - u_1 - \alpha u_2 - \epsilon v), \quad (1) $$

$$ \frac{du_2}{dt} = u_2(b_2 - \beta u_1 - u_2 - \mu v), \quad (2) $$

$$ \frac{dv}{dt} = v(-b_3 + d\epsilon u_1 + d\mu u_2), \quad (3) $$

where $u_1(t)$, $u_2(t)$ and $v(t)$ represent the densities of two prey species and one predator species at time $t$, respectively. The parameters $b_i$ ($i = 1, 2, 3$) involved in this system denote the intrinsic rates of increase or decrease in density for the considered three species. The rates of competitive effects between the two prey species are represented by the parameters $\alpha$ and $\beta$, while the parameters $\epsilon$ and $\mu$ account for the diminishing rates in density of both the prey species due to predation by the considered predator. The parameter $d$ denotes the transformation rate of prey biomass into predator biomass.
Note that all the parameters involved in this prey–predator system (1)–(3) are positive constants. Also, the model (1)–(3) is subjected to the non-negative initial conditions.

The possible equilibrium points for the temporal model (1)–(3) can be explicitly derived from the system of algebraic equations arising by equating the reaction parts to zero. The system (1)–(3) admits the following six trivial and semi-trivial equilibrium points: $E_{000} = (0,0,0)$, $E_{+00} = (b_1,0,0)$, $E_{0+0} = (0,b_2,0)$, $E_{++0} = (\frac{b_1-\alpha b_2}{1-\alpha\beta}, \frac{b_2-\beta b_1}{1-\alpha\beta},0)$, $E_{+0+} = (\frac{b_3}{d\epsilon},0, \frac{db_1\epsilon-b_3}{d\epsilon^2})$ and $E_{0++} = (0, \frac{b_3}{d\mu}, \frac{db_2\mu-b_3}{d\mu^2})$. Also, this system admits a unique coexistence equilibrium point $E_{+++} = (u_1^*, u_2^*, v^*)$, where

$$ 
\begin{aligned}
u_1^* &= \frac{b_3\epsilon - db_2\epsilon\mu - \alpha b_3\mu + db_1\mu^2}{d(\epsilon^2 + \mu^2 - (\alpha + \beta)\epsilon\mu)}, \\
u_2^* &= \frac{b_3\mu - db_1\epsilon\mu - \beta b_3\epsilon + db_2\epsilon^2}{d(\epsilon^2 + \mu^2 - (\alpha + \beta)\epsilon\mu)}, \\
v^* &= \frac{b_3(\alpha\beta - 1) + d\mu(b_2 - \beta b_1) + d\epsilon(b_1 - \alpha b_2)}{d(\epsilon^2 + \mu^2 - (\alpha + \beta)\epsilon\mu)}.
\end{aligned}
$$

The feasibility conditions for the coexistence equilibrium point $E_{+++}$ are given by (i) $b_3\epsilon + db_1\mu^2 > \mu(db_2\epsilon + \alpha b_3)$, $b_3\mu + db_2\epsilon^2 > \epsilon(db_1\mu + \beta b_3)$, $b_3\alpha\beta + d(b_1\epsilon + b_2\mu) > b_3 + d(b_2\alpha\epsilon + b_1\beta\mu)$ and $\epsilon^2 + \mu^2 > (\alpha + \beta)\epsilon\mu$, or (ii) $b_3\epsilon + db_1\mu^2 < \mu(db_2\epsilon + \alpha b_3)$, $b_3\mu + db_2\epsilon^2 < \epsilon(db_1\mu + \beta b_3)$, $b_3\alpha\beta + d(b_1\epsilon + b_2\mu) < b_3 + d(b_2\alpha\epsilon + b_1\beta\mu)$ and $\epsilon^2 + \mu^2 < (\alpha + \beta)\epsilon\mu$.

In this study, we are mainly concerned with the dynamical behaviors of the system around the coexistence equilibrium point since the existence of this equilibrium point indicates the species diversity. For detailed theoretical and numerical investigations on the dynamics of this system in the vicinity of the coexistence equilibrium point, one can go through the references [5,56]. However, here we briefly describe the theoretical results on local stability and Hopf bifurcation of the coexistence equilibrium point $E_{+++}$. Linearizing the above model (1)–(3) around the equilibrium point $E_{+++}$, we obtain the following characteristic equation

$$ \lambda^3 + a_2\lambda^2 + a_1\lambda + a_0 = 0, \quad (4) $$

where

$$ a_2 = u_1^* + u_2^*, \quad (5) $$

$$ a_1 = (1 - \alpha\beta)u_1^*u_2^* + d(\epsilon^2 u_1^* + \mu^2 u_2^*)v^*, \quad (6) $$

$$ a_0 = d\{(\epsilon^2 + \mu^2) - (\alpha + \beta)\epsilon\mu\}u_1^*u_2^*v^*. \quad (7) $$

Therefore, the equilibrium point $E_{+++}$ is locally asymptotically stable provided the following Routh–Hurwitz criteria are satisfied:

$$ a_2 > 0, \quad a_0 > 0, \quad a_2a_1 - a_0 > 0. \quad (8) $$

The condition $a_2 > 0$ is automatically satisfied for feasible $E_{+++}$. However, we need to have the parametric restriction $\epsilon^2 + \mu^2 > (\alpha + \beta)\epsilon\mu$ in order to satisfy the condition $a_0 > 0$. Therefore, the parametric restriction $\epsilon^2 + \mu^2 > (\alpha + \beta)\epsilon\mu$ acts as a necessary condition for local stability. Furthermore, the model (1)–(3) undergoes Hopf bifurcation at $\epsilon = \epsilon_*$ when the following conditions are satisfied:

$$ (a_2a_1 - a_0)|_{\epsilon=\epsilon_*} = 0, \quad \frac{d}{d\epsilon}(a_2a_1 - a_0)|_{\epsilon=\epsilon_*} \neq 0, \quad (9) $$

with the parametric restriction $\epsilon^2 + \mu^2 > (\alpha + \beta)\epsilon\mu$.

Now, we encapsulate numerically the dynamical behaviors around the coexistence equilibrium point $E_{+++}$ in the form of a bifurcation diagram with respect to the parameter accounting for the
diminishing rate in first prey density due to predation (ε). For this purpose, we consider a fixed set of parameter values

$$b_1 = 1.0, \alpha = 1.0, b_2 = 1.0, \beta = 1.5, \mu = 1.0, b_3 = 1.0, d = 0.5, \qquad (10)$$

and vary the parameter ε. The resulting dynamical behaviors are encapsulated in Figure 1. From this figure, we can observe that the stable coexistence equilibrium point E$_{+++}$ loses its stability through the occurrence of Hopf bifurcation approximately at ε = 5.486 and stable limit cycles appear for higher ε-values than this threshold value. The size of the limit cycle gradually increases as the ε-value moves away in positive direction from this threshold. For sufficiently large values of ε (for example, ε = 8.0), the system admits two-periodic solutions. Further increments in the value of ε result in chaotic dynamics (for example, ε = 10.0). This figure also clearly illustrates that the system takes the period-doubling route to reach chaos finally.

<figure><img src="image_5.png" /><figcaption>Figure 1. Single parameter bifurcation diagram for the density of first prey population (u₁) with respect to the parameter ε for the temporal model (1)–(3). Other parameter values are mentioned in (10).</figcaption></figure>

### 3. Spatiotemporal Model

In this section, we extend the temporal model (1)–(3) presented in the preceding section by incorporating the random dispersal of all the three species. The corresponding spatiotemporal model is given by

$$\frac{\partial u_1}{\partial t} = d_1 \frac{\partial^2 u_1}{\partial x^2} + u_1 (b_1 - u_1 - \alpha u_2 - \epsilon v), \qquad (11)$$

$$\frac{\partial u_2}{\partial t} = d_2 \frac{\partial^2 u_2}{\partial x^2} + u_2 (b_2 - \beta u_1 - u_2 - \mu v), \qquad (12)$$

$$\frac{\partial v}{\partial t} = d_3 \frac{\partial^2 v}{\partial x^2} + v(-b_3 + d\epsilon u_1 + d\mu u_2), \qquad (13)$$

subjected to non-negative initial conditions and periodic boundary conditions. Here, u₁(x, t), u₂(x, t) and v(x, t) respectively denote the densities of the three species at spatial position x and time t. The bounded spatial domain is given by [-L, L], where L(> 0) ∈ ℝ. The positive parameters dᵢ (i = 1, 2, 3) represent the diffusion coefficients for two prey and one predator species, respectively.
Note that the equilibrium points corresponding to the temporal model (1)–(3) also serve as the spatially homogeneous steady states for the model (11)–(13).

Turing instability occurs when the stable spatially homogeneous steady state becomes unstable due to small amplitude heterogeneous perturbations around the spatially homogeneous steady state. Introducing small amplitude spatiotemporal perturbation around the spatially homogeneous coexistence steady state $E_{+++} = (u_1^*, u_2^*, v^*)$ and then linearizing we find the following Jacobian matrix:

$$ \mathcal{J}(k^2) = \begin{bmatrix} -u_1^* - d_1 k^2 & -\alpha u_1^* & -\epsilon u_1^* \\ -\beta u_2^* & -u_2^* - d_2 k^2 & -\mu u_2^* \\ d\epsilon v^* & d\mu v^* & -d_3 k^2 \end{bmatrix}, \quad (14) $$

where $k$ represents the wavenumber. Therefore, the characteristic equation is given by

$$ \lambda^3 + A(k^2)\lambda^2 + B(k^2)\lambda + C(k^2) = 0, \quad (15) $$

where

$$ A(k^2) = (d_1 + d_2 + d_3)k^2 + u_1^* + u_2^*, $$

$$ B(k^2) = (d_1d_2 + d_2d_3 + d_3d_1)k^4 + \{d_1u_2^* + d_2u_1^* + d_3(u_1^* + u_2^*)\}k^2 + (1 - \alpha\beta)u_1^*u_2^* + d(\epsilon^2u_1^* + \mu^2u_2^*)v^*, $$

$$ C(k^2) = d_1d_2d_3k^6 + d_3(d_1u_2^* + d_2u_1^*)k^4 + \{d_1d\mu^2u_2^*v^* + d_2d\epsilon^2u_1^*v^* + d_3(1 - \alpha\beta)u_1^*u_2^*\}k^2 + \{d(\epsilon^2 + \mu^2) - d(\alpha + \beta)\epsilon\mu\}u_1^*u_2^*v^*. $$

Now, Re($\lambda$) is less than zero provided the Routh-Hurwitz criteria are satisfied:

$$ A(k^2) > 0, C(k^2) > 0, A(k^2)B(k^2) - C(k^2) > 0. \quad (16) $$

When all the conditions presented in (16) are satisfied, then the spatially homogeneous steady state is stable. Diffusion-driven instability occurs only when one eigenvalue is zero while other two eigenvalues still have negative real parts [21,57]. If $\lambda_1, \lambda_2$ and $\lambda_3$ denote the roots of the characteristic Equation (15), then from the properties of the roots of a cubic equation we have the following relations:

$$ \lambda_1 + \lambda_2 + \lambda_3 = -A(k^2), \quad (17) $$

$$ \lambda_1\lambda_2 + \lambda_2\lambda_3 + \lambda_3\lambda_1 = B(k^2), \quad (18) $$

$$ \lambda_1\lambda_2\lambda_3 = -C(k^2), \quad (19) $$

$$ -(\lambda_1 + \lambda_2)(\lambda_2 + \lambda_3)(\lambda_3 + \lambda_1) = A(k^2)B(k^2) - C(k^2). \quad (20) $$

Below the Turing bifurcation threshold, the spatially homogeneous steady state is stable and accordingly all the eigenvalues have negative real parts. At the Turing bifurcation threshold, exactly one eigenvalue becomes zero at the critical wavenumber $k_T$, whereas other two eigenvalues still have negative real parts [58]. Since at the critical wavenumber $k = k_T$ we have one of the roots is equal to zero, without any loss of generality we assume

$$ \lambda_1 |_{k^2=k_T^2} = 0, \quad \text{Re}(\lambda_2) |_{k^2=k_T^2} < 0 \quad \text{and} \quad \text{Re}(\lambda_3) |_{k^2=k_T^2} < 0. \quad (21) $$

Hence, at the critical wavenumber $k = k_T$ we obtain $C(k_T^2) = 0$. Further due to the above conditions (21) for Turing instability, we get $A(k_T^2) > 0, B(k_T^2) > 0$ and $A(k_T^2)B(k_T^2) - C(k_T^2) = A(k_T^2)B(k_T^2) > 0$. Thus, the system remains stable if $C(k^2) > 0$ holds for all $k$ and it becomes Turing unstable if $C(k^2) < 0$ holds for at least one $k$. Also, $C(k^2) < 0$ holds for a range of $k$-values around $k_T$ beyond the Turing bifurcation threshold. The expression of $C(k^2)$ can be rewritten as

$$ C(k^2) \equiv C_3(k^2)^3 + C_2(k^2)^2 + C_1(k^2) + C_0, \quad (22) $$
where $C_3 > 0$ since the diffusion coefficients are positive constants and $C_0 > 0$ from the local asymptotic stability condition of $E_{+++}$. Also, the coefficient $C_2 > 0$ as the expression of it is given by $C_2 = d_3(d_1u_2^* + d_2u_1^*)$. The minimum for $C(k^2)$ occurs at $k = k_T$ where

$$k_T^2 = \frac{-C_2 + \sqrt{C_2^2 - 3C_1C_3}}{3C_3}, \qquad (23)$$

as we have $\frac{dC(k^2)}{d(k^2)} = 0$ and $\frac{d^2C(k^2)}{d(k^2)^2} > 0$ at $k = k_T$. Now from the expression (23), we can notice that $k_T^2$ is positive if $C_1 < 0$ or $C_2 < 0$ and $C_2^2 > 3C_1C_3$. Since $C_2 > 0$ holds for our considered system, then the condition for the positivity of $k_T^2$ reduces to $C_1 < 0$. Therefore, the Turing bifurcation boundary is given by

$$2C_2^3 - 9C_1C_2C_3 - 2(C_2^2 - 3C_1C_3)^{\frac{3}{2}} + 27C_0C_3^2 = 0. \qquad (24)$$

This is an implicit expression for the Turing bifurcation curve whenever other relevant conditions are satisfied.

## Numerical Results

In this subsection, we present some numerical results in order to validate the analytical predictions and to manifest the spatiotemporal dynamics which are beyond the scope of linear analysis. For this purpose, we considered a set of parameter values

$$b_1 = 1.0, \quad \alpha = 1.0, \quad b_2 = 1.0, \quad \beta = 1.5, \quad \mu = 1.0, \quad d = 0.5, \quad b_3 = 1.0, \quad d_1 = 0.1, \quad d_2 = 0.2, \qquad (25)$$

and varied the parameters $\epsilon$ and $d_3$. For this set of parameter values, the temporal-Hopf bifurcation threshold is approximately given by $\epsilon = 5.486$. First, we present the temporal-Hopf and Turing bifurcation curves for the local spatiotemporal model (11)–(13) in $(\epsilon, d_3)$-parameter space in Figure 2. From this figure, we can notice that these two bifurcation curves divide the two-dimensional parameter space into four different regions, such as stable region (left to the temporal-Hopf curve and below the Turing curve), Turing region (left to the temporal-Hopf curve and above the Turing curve), Hopf-Turing region (right to the temporal-Hopf curve and above the Turing curve) and pure Hopf region (right to the temporal-Hopf curve and below the Turing curve).

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Plots of temporal-Hopf and Turing bifurcation curves in $(\epsilon, d_3)$-parameter space for the local spatiotemporal model (11)–(13). Here, red vertical straight line represents the temporal-Hopf bifurcation threshold and the blue curve represents the Turing bifurcation curve. Other parameter values are mentioned in (25). The red vertical straight line has been drawn by solving the equations presented in (9) simultaneously and the blue curve has been drawn from the Equation (24).</figcaption></figure>
Before proceeding further, here we briefly describe about the numerical discretization of the system (11)–(13). We have used the forward Euler scheme and central difference scheme to discretize the reaction and diffusion parts, respectively. The periodic boundary conditions have been employed at the boundary of the computational domain. For all the patterns presented in this subsection, we have used the following pulse-type initial conditions:

$$ 
\begin{align} 
u_1(x, 0) &= \begin{cases} u_1^* + \eta_1, & |x| < 10, \\ u_1^*, & 10 \le |x| \le L, \end{cases} \nonumber \\ 
u_2(x, 0) &= \begin{cases} u_2^* + \eta_2, & |x| < 10, \\ u_2^*, & 10 \le |x| \le L, \end{cases} \tag{26} \\ 
v(x, 0) &= \begin{cases} v^* + \eta_3, & |x| < 10, \\ v^*, & 10 \le |x| \le L, \end{cases} \nonumber 
\end{align} 
$$

where $|\eta_r| \ll 1$, for $r = 1, 2, 3$.

Now, we consider a point $(\epsilon, d_3) = (5.4, 9.0)$, which lies in the Turing region. For this choice of $\epsilon$, we obtain the coexistence steady state $E_{+++} = (0.2641, 0.5738, 0.03)$ which is locally asymptotically stable in the absence of diffusion. Figure 3 shows the stationary Turing pattern produced by $u_1$-component and the corresponding temporal evolution of the average densities of all considered species. We would like to mention here that the patterns for the second prey and predator species are qualitatively similar with the pattern for the first prey species presented in Figure 3a. However, the regions with high first prey density correspond to the same for predator but to low second prey density. These scenarios for the density correlation among the three considered species are presented in Figure 4. This type of density correlation among the three considered species has also been observed for other patterns presented in this section. Also, note that we have considered the one-dimensional spatial domain $[-100, 100]$ for the clarity of understanding in Figure 4, but all the patterns presented in this section are obtained for the spatial domain $[-200, 200]$.

<figure><img src="image_7.png" /><figcaption>Figure 3. (a) Plot of stationary Turing pattern exhibited by $u_1$, (b) time evolution of the average densities of three species. $\epsilon = 5.4$ and $d_3 = 9.0$, other parameter values are mentioned in (25).</figcaption></figure>

Now if we choose $\epsilon = 6.0$, then $E_{+++} = (0.2273, 0.6364, 0.0227)$ is unstable, due to temporal-Hopf bifurcation, in the absence of diffusion terms. Hence, we can observe homogeneous in space and periodic in time pattern for $d_3 = 3.0$ in Figure 5a, periodic in both space and time pattern for $d_3 = 3.75$ in Figure 5b, and periodic in space but stationary in time pattern for $d_3 = 9.0$ in Figure 5c. Phase diagrams presented in Figure 6 confirm the temporal periodicity of the patterns presented in Figure 5a,b, respectively. In Figure 7, we present the bifurcation diagram for spatially averaged
density of first prey species in order to illustrate the transition of different types of patterns depending on the value of the diffusion coefficient $d_3$ for $\epsilon = 6.0$. In order to prepare this bifurcation diagram, we have used small pulse type perturbation about the mid-point around the spatially homogeneous steady state for each $d_3$-value and then plotted the maxima and minima of spatially averaged densities. From this figure, we can observe that periodic in time and homogeneous in space pattern persists up to $d_3 = 3.67$. The periodic in both space and time pattern exists for $d_3 \in [3.68, 3.91]$ and then stationary pattern emerges for $d_3 > 3.91$. Interestingly, the number of patches with high first prey density gradually decreases with the increment in the value of $d_3$.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Correlation of densities among two prey and one predator species at a particular instant of time when the patterns are stationary in time. Left panel (a) shows the density correlation between first prey and predator species and right panel (b) shows the same between second prey and predator species. Numerical simulation has been performed for the same parameter values as in the Figure 3 over the one-dimensional spatial domain [-100, 100].</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Plots of different types of patterns exhibited by $u_1$ for $\epsilon = 6.0$. Three different patterns are obtained for (a) $d_3 = 3.0$; (b) $d_3 = 3.75$; (c) $d_3 = 9.0$. Other parameter values are mentioned in (25).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Phase diagrams of spatial average densities of the associated three species corresponding to (a) Figure 5a, (b) Figure 5b.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Single parameter bifurcation diagram for the spatially averaged density of first prey population ($<u_1>$) with respect to the diffusion coefficient $d_3$. Here, $\epsilon = 6.0$ and other parameter values are mentioned in (25).</figcaption></figure>

If we further increase the magnitude of $\epsilon$ to 8.0, then the system (11)–(13) admits the spatially homogeneous coexistence steady state $E_{+++} = (0.1556, 0.7556, 0.0111)$ which is temporally unstable. Then the choice $d_3 = 1.0$ leads to the spatiotemporal chaotic pattern (see Figure 8a). A close look at the chaotic pattern reveals that the distribution of the species is symmetric about the midpoint in space and this is due to our consideration of the symmetric pulse type initial conditions. The chaotic nature has been cross verified by checking the sensitivity to initial conditions as described in [59], but we do not present the results here for the sake of brevity. Interestingly, spatiotemporal chaotic pattern disappears and eventually settles to a stationary pattern for a higher value of $d_3$ (for example, $d_3 = 10.0$). This transition from spatiotemporal chaotic to stationary pattern takes place via a pattern where patches move periodically for intermediate values of $d_3$ (for example, $d_3 = 9.5$). The resulting patterns are presented in Figure 8. The phase diagrams presented in Figure 9 confirm the chaotic nature of the pattern presented in Figure 8a and the periodic nature of the pattern presented in Figure 8b.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Plots of different types of patterns exhibited by $u_1$ for $\epsilon = 8.0$. Three different patterns are obtained for (a) $d_3 = 1.0$; (b) $d_3 = 9.5$; (c) $d_3 = 10.0$. Other parameter values are mentioned in (25).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Phase diagrams of spatial average densities of the associated three species corresponding to (a) Figure 8a, (b) Figure 8b.</figcaption></figure>

All patterns exhibited in Figures 3a, 5 and 8 are stable in the sense that their nature will not alter with further increment in time. However, the considered system (11)–(13) can produce long transient patterns in order to reach the final stable patterns. Such an instance of long transient pattern can be found for $\epsilon = 8.0$ and $d_3 = 9.0$ along with other parameter values mentioned above. Figure 10 exhibits the phase diagrams for the spatially averaged densities of the considered three species for $t \in [18,000, 20,000]$ and $t \in [58,000, 60,000]$, respectively. Since the phase diagram presented in Figure 10a is restricted to an annular region, we can term the nature of the corresponding spatiotemporal distribution as quasi-periodic which persists for approximately $t = 50,000$. On the other hand, the phase diagram presented in Figure 10b clearly shows the periodic nature of the
corresponding spatiotemporal distribution and this pattern is stable for the considered parameter values. With the increment in time, the width of the annular region gradually decreases and finally the region settles into a limit cycle. However, this long time period for the existence of quasi-stable and quasi-periodic pattern can be few or even hundreds of generations for certain species. Therefore, this type of existence of long transient pattern and then finally convergence into a stable pattern can provide an alternative explanation for ecological regime shifts apart from the usual concept of exogenously sparked regime shifts.

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Phase diagrams of spatial average densities of the associated three species for $\epsilon = 8.0$ and $d_3 = 9.0$ with (a) $t \in [18,000,20,000]$, (b) $t \in [58,000,60,000]$. Other parameter values are mentioned in (25).</figcaption></figure>

Now, we want to investigate the influence of random dispersal of the three species on the dynamics presented in Figure 1 for the temporal model. For this purpose, we considered the same parameter values as for Figure 1 with $d_1 = 0.1$, $d_2 = 0.2$ and $d_3 = 1.0$, and varied the value of the parameter $\epsilon$. By employing spatially homogeneous initial population distributions slightly perturbed from the spatially homogeneous coexistence steady state, we found the bifurcation diagram which was more or less similar to the Figure 1 and chose not to display it here for the sake of brevity. However, consideration of pulse type initial population distributions led to different dynamical behaviors for a certain range of $\epsilon$-values and the resulting dynamics are encapsulated in Figure 11. Comparing Figures 1 and 11, we can observe that the dynamics for both the models were same up to $\epsilon = 8.5$ approximately. For higher values of $\epsilon$ (that is, $\epsilon \in (8.5, 10.0]$), we observed the difference in chaotic dynamics. For spatiotemporal model (11)–(13), the amplitudes of chaotic oscillations were much less than that of the temporal case. We can attribute the reason for this type of phenomenon to the emergence of spatial heterogeneity in the population distributions as we have found spatially heterogeneous chaotic distributions in this case.

Further, we examined the stability of this bifurcation diagram by employing both the forward and backward continuation techniques. In order to construct a bifurcation diagram using forward continuation technique, we employed the following algorithm:

* First, consider a pulse-type initial condition for the starting value of the bifurcation parameter and run the simulation for sufficient time period such that a stable state is reached.

* Make a tiny increment in the value of the bifurcation parameter.

* Then, use the stable state obtained in the previous step as an initial condition for the present value of the bifurcation parameter and run the simulation till a stable state is reached.

* Continue the previous two steps till the last value of the bifurcation parameter is reached.
In a similar manner, one can efficiently construct a bifurcation diagram using backward continuation technique by reversing the process encapsulated in the above algorithm. Forward continuation technique results in more or less similar bifurcation diagram as the Figure 11 and accordingly we do not include it here. On the other hand, we observed significant changes in the bifurcation diagram using backward continuation technique and the resulting dynamics are encapsulated in Figure 12. From this diagram, we can observe that the period-doubling route to chaos has been destroyed and spatially heterogeneous chaotic regime survives for larger range of $\epsilon$-values. Overall, bifurcation diagrams presented in Figures 1, 11 and 12 indicate the existence of alternative stable states for the model (11)–(13). Note that these states indicate the coexistence of all the three species and can be both the spatially homogeneous and heterogeneous population distributions.

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Single parameter bifurcation diagram for the spatially averaged density of first prey population ($<u_1>$) with respect to the parameter $\epsilon$ for the spatiotemporal model (11)–(13). Here, we have used pulse type initial population distributions for each value of $\epsilon$. Other parameter values are mentioned in (25) with $d_3 = 1.0$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Single parameter bifurcation diagram for the spatially averaged density of first prey population ($<u_1>$) with respect to the parameter $\epsilon$ for the spatiotemporal model (11)–(13). Here, backward continuation technique has been employed. Other parameter values are mentioned in (25) with $d_3 = 1.0$.</figcaption></figure>
Another interesting finding is that the spatial population distributions depended on the spatial
domain size when the parameter values were taken far from the temporal-Hopf bifurcation threshold.
The patterns presented in Figure 13 are numerically computed over the one-dimensional spatial domain
[-50, 50] with the same parameter values used for Figure 8a,c. Figure 13a suggests that the smaller
domain size can produce a two-periodic pattern while the larger domain size produces a chaotic one.
Also, Figure 13c indicates that the smaller domain size can give rise to periodic in both the space and
time population distribution while the larger domain size gives a stationary heterogeneous distribution.
The two-periodic and periodic nature of the patterns presented in Figure 13a,c can be confirmed by the
phase portraits of the spatially averaged densities in Figure 13b,d, respectively. However, we did not
find any effect of the domain size on the resulting patterns presented in Figures 3, 5 and 8b. Therefore,
our numerical findings suggest that not only certain parameters associated with a model play crucial
roles in the emergence of different types of spatial population distributions but also the spatial domain
size can potentially play a key role in shaping the population distributions.

<figure><img src="image_17.png" /><figcaption>**Figure 13.** (a) Plot of spatial pattern exhibited by $u_1$, (b) corresponding phase portrait of the average densities of three species for $\epsilon = 8.0$, $d_3 = 1.0$. (c) Plot of spatial pattern exhibited by $u_1$, (d) corresponding phase portrait of the average densities of three species for $\epsilon = 8.0$, $d_3 = 10.0$. Spatial domain size is [-50, 50] and other parameter values are mentioned in (25).</figcaption></figure>
## 4. Nonlocal Model

In this section, we further extend the spatiotemporal model (11)–(13) by incorporating nonlocal intra-specific competition for both the prey species. Accordingly, the model (11)–(13) becomes

$$ \frac{\partial u_1}{\partial t} = d_1 \frac{\partial^2 u_1}{\partial x^2} + u_1 (b_1 - J(u_1) - \alpha u_2 - \epsilon v), \quad (27) $$

$$ \frac{\partial u_2}{\partial t} = d_2 \frac{\partial^2 u_2}{\partial x^2} + u_2 (b_2 - \beta u_1 - J(u_2) - \mu v), \quad (28) $$

$$ \frac{\partial v}{\partial t} = d_3 \frac{\partial^2 v}{\partial x^2} + v(-b_3 + d\epsilon u_1 + d\mu u_2), \quad (29) $$

where

$$ J(u_r) = \int_{-\infty}^{\infty} u_r(y,t) \phi_r(x-y) dy, \quad r = 1,2, $$

which take care of the nonlocality in intra-specific interactions for both the prey species arising primarily as a result of nonlocal consumption of resources. However, there exist other mechanisms such as cannibalism, space sharing, fights, etc., which can eventually lead to nonlocal intra-specific competition. The nonlocal system (27)–(29) is subjected to the same initial and boundary conditions mentioned before for the local spatiotemporal system (11)–(13). Here, the functions $\phi_r$ ($r = 1,2$) represent the non-negative kernel functions. Further, in order to capture the unbiased movement of the two prey populations we assume the kernel functions $\phi_r$ ($r = 1,2$) as even functions. Throughout this paper, we restrict ourselves only to the top-hat kernel functions defined as follows:

$$ \phi_r(x) = \begin{cases} \frac{1}{2\delta_r}, & |x| \le \delta_r, \\ 0, & \text{otherwise,} \end{cases} \quad (30) $$

for $r = 1,2$. The newly introduced parameters $\delta_r$ ($r = 1,2$) represent the extent of the nonlocality by controlling the width of the kernel functions $\phi_r$. Also, one can easily verify that the spatially homogeneous steady states of the nonlocal system (27)–(29) are same as that of the corresponding local system (11)–(13) and accordingly we use the same notations for them.

In order to linearize the nonlocal system (27)–(29) about the spatially homogeneous coexistence steady state $E_{+++} = (u_1^*, u_2^*, v^*)$, we introduce the perturbation $u_1 = u_1^* + \epsilon \tilde{u}_1 e^{\lambda t + ikx}$, $u_2 = u_2^* + \epsilon \tilde{u}_2 e^{\lambda t + ikx}$ and $v = v^* + \epsilon \tilde{v} e^{\lambda t + ikx}$ with $|\epsilon| \ll 1$. Then the Jacobian matrix corresponding to the linearized system is given by

$$ \begin{aligned} \tilde{J}(k, \delta_1, \delta_2) &= \begin{bmatrix} -u_1^* \int_{-\infty}^{\infty} \phi_1(y) e^{-iky} dy - d_1 k^2 & -\alpha u_1^* & -\epsilon u_1^* \\ -\beta u_2^* & -u_2^* \int_{-\infty}^{\infty} \phi_2(y) e^{-iky} dy - d_2 k^2 & -\mu u_2^* \\ d\epsilon v^* & d\mu v^* & -d_3 k^2 \end{bmatrix} \\ &= \begin{bmatrix} -u_1^* \frac{\sin(k\delta_1)}{k\delta_1} - d_1 k^2 & -\alpha u_1^* & -\epsilon u_1^* \\ -\beta u_2^* & -u_2^* \frac{\sin(k\delta_2)}{k\delta_2} - d_2 k^2 & -\mu u_2^* \\ d\epsilon v^* & d\mu v^* & -d_3 k^2 \end{bmatrix}. \end{aligned} \quad (31) $$

Therefore, the characteristic equation is given by

$$ \lambda^3 + \tilde{A}(k, \delta_1, \delta_2)\lambda^2 + \tilde{B}(k, \delta_1, \delta_2)\lambda + \tilde{C}(k, \delta_1, \delta_2) = 0, \quad (32) $$
where

$$ \tilde{A}(k, \delta_1, \delta_2) = (d_1 + d_2 + d_3)k^2 + u_1^* \frac{\sin(k\delta_1)}{k\delta_1} + u_2^* \frac{\sin(k\delta_2)}{k\delta_2}, \quad (33) $$

$$ \begin{aligned} \tilde{B}(k, \delta_1, \delta_2) &= d_3 k^2 \left\{ (d_1 + d_2)k^2 + u_1^* \frac{\sin(k\delta_1)}{k\delta_1} + u_2^* \frac{\sin(k\delta_2)}{k\delta_2} \right\} + \left( d_1 k^2 + u_1^* \frac{\sin(k\delta_1)}{k\delta_1} \right) \left( d_2 k^2 + u_2^* \frac{\sin(k\delta_2)}{k\delta_2} \right) \\ &\quad + d(\epsilon^2 u_1^* + \mu^2 u_2^*) v^* - \alpha \beta u_1^* u_2^* \\ &= (d_1 d_2 + d_2 d_3 + d_3 d_1)k^4 + \left\{ d_1 u_2^* \frac{\sin(k\delta_2)}{k\delta_2} + d_2 u_1^* \frac{\sin(k\delta_1)}{k\delta_1} + d_3 \left( u_1^* \frac{\sin(k\delta_1)}{k\delta_1} + u_2^* \frac{\sin(k\delta_2)}{k\delta_2} \right) \right\} k^2 \\ &\quad + \left( \frac{\sin(k\delta_1) \sin(k\delta_2)}{k^2 \delta_1 \delta_2} - \alpha \beta \right) u_1^* u_2^* + d(\epsilon^2 u_1^* + \mu^2 u_2^*) v^*, \end{aligned} \quad (34) $$

$$ \begin{aligned} \tilde{C}(k, \delta_1, \delta_2) &= d_3 k^2 \left( d_1 k^2 + u_1^* \frac{\sin(k\delta_1)}{k\delta_1} \right) \left( d_2 k^2 + u_2^* \frac{\sin(k\delta_2)}{k\delta_2} \right) + d\epsilon^2 u_1^* v^* \left( d_2 k^2 + u_2^* \frac{\sin(k\delta_2)}{k\delta_2} \right) + \\ &d\mu^2 u_2^* v^* \left( d_1 k^2 + u_1^* \frac{\sin(k\delta_1)}{k\delta_1} \right) - d_3 \alpha \beta u_1^* u_2^* k^2 - d(\alpha + \beta) \epsilon \mu u_1^* u_2^* v^* \\ &= d_1 d_2 d_3 k^6 + d_3 \left\{ d_1 u_2^* \frac{\sin(k\delta_2)}{k\delta_2} + d_2 u_1^* \frac{\sin(k\delta_1)}{k\delta_1} \right\} k^4 + \\ &\left\{ d_1 d\mu^2 u_2^* v^* + d_2 d\epsilon^2 u_1^* v^* + d_3 \left( \frac{\sin(k\delta_1) \sin(k\delta_2)}{k^2 \delta_1 \delta_2} - \alpha \beta \right) u_1^* u_2^* \right\} k^2 \\ &+ d \left\{ \left( \epsilon^2 \frac{\sin(k\delta_2)}{k\delta_2} + \mu^2 \frac{\sin(k\delta_1)}{k\delta_1} \right) - (\alpha + \beta) \epsilon \mu \right\} u_1^* u_2^* v^*. \end{aligned} \quad (35) $$

Using the Routh-Hurwitz criteria, we can say that the spatially homogeneous coexistence steady state $E_{+++} = (u_1^*, u_2^*, v^*)$ will be asymptotically stable if and only if

$$ \tilde{A}(k, \delta_1, \delta_2) > 0, \quad \tilde{C}(k, \delta_1, \delta_2) > 0, \quad \tilde{A}(k, \delta_1, \delta_2) \tilde{B}(k, \delta_1, \delta_2) - \tilde{C}(k, \delta_1, \delta_2) > 0. \quad (36) $$

In order to attain Turing instability threshold we need to have one zero eigenvalue and other two eigenvalues with negative real parts. Following the process described in the previous section, the conditions for achieving Turing instability threshold imply the following conditions:

$$ \tilde{A}(k, \delta_1, \delta_2) > 0, \quad \tilde{B}(k, \delta_1, \delta_2) > 0, \quad \tilde{C}(k, \delta_1, \delta_2) = 0. \quad (37) $$

Note that the expressions of $\tilde{A}$, $\tilde{B}$ and $\tilde{C}$ are transcendental in nature of their arguments $k$, $\delta_1$ and $\delta_2$. The associated Turing bifurcation threshold can be obtained by solving the following two equations simultaneously

$$ \tilde{C}(k, \delta_1, \delta_2) = 0, \quad \frac{\partial}{\partial k} \tilde{C}(k, \delta_1, \delta_2) = 0. \quad (38) $$

Solving the first equation of (38) for $d_3$, we obtain

$$ d_3 = - \frac{d \left\{ (d_1 \mu^2 u_2^* + d_2 \epsilon^2 u_1^*) k^2 + \left( \epsilon^2 \frac{\sin(k\delta_2)}{k\delta_2} + \mu^2 \frac{\sin(k\delta_1)}{k\delta_1} - (\alpha + \beta) \epsilon \mu \right) u_1^* u_2^* \right\} v^*}{d_1 d_2 k^6 + \left\{ d_1 u_2^* \frac{\sin(k\delta_2)}{k\delta_2} + d_2 u_1^* \frac{\sin(k\delta_1)}{k\delta_1} \right\} k^3 + \left\{ \frac{\sin(k\delta_1) \sin(k\delta_2)}{k^2 \delta_1 \delta_2} - \alpha \beta \right\} u_1^* u_2^* k^2}. \quad (39) $$

Now, differentiating $\tilde{C}(k, \delta_1, \delta_2)$ with respect to $k$, we get

$$ \begin{aligned} \frac{\partial \tilde{C}}{\partial k} = d_3 & \left[ 6d_1 d_2 k^5 + \left\{ d_1 u_2^* \cos(k\delta_2) + d_2 u_1^* \cos(k\delta_1) \right\} k^3 + 3 \left\{ d_1 u_2^* \frac{\sin(k\delta_2)}{k\delta_2} + d_2 u_1^* \frac{\sin(k\delta_1)}{k\delta_1} \right\} k^2 - \right. \\ & \left. 2\alpha \beta u_1^* u_2^* k + \left\{ \frac{\cos(k\delta_1) \sin(k\delta_2)}{k\delta_2} + \frac{\sin(k\delta_1) \cos(k\delta_2)}{k\delta_1} \right\} u_1^* u_2^* \right] + d \left[ 2(d_1 \mu^2 u_2^* + d_2 \epsilon^2 u_1^*) k + \right. \\ & \left. \left\{ \epsilon^2 \left( \frac{\cos(k\delta_2)}{k} - \frac{\sin(k\delta_2)}{k^2 \delta_2} \right) + \mu^2 \left( \frac{\cos(k\delta_1)}{k} - \frac{\sin(k\delta_1)}{k^2 \delta_1} \right) \right\} u_1^* u_2^* \right] v^*. \end{aligned} \quad (40) $$

Substituting the expression for $d_3$ in the second equation of (38), we obtain

$$ \mathcal{F}(k, \delta_1, \delta_2) = 0, \quad (41) $$
where the expression for $\mathcal{F}$, which is transcendental in nature, is given in Appendix A. Solving the Equation (41) for $k$, we obtain the critical wavenumber $k_T$. But the transcendental nature of the Equation (41) prevents us from finding an analytical solution for $k_T$ and hence, we find the critical value numerically. Further, substituting the value of $k_T$ into the Equation (39) we obtain the critical value of the diffusion coefficient of predator $d_3^T$ in order to induce Turing instability.

## Numerical Results

In this subsection, we first investigate the influence of the diffusion coefficients on the positioning of the Turing curve. As the expressions for Turing instability conditions for both the local and nonlocal models involve the diffusion coefficients (that is, $d_1$, $d_2$ and $d_3$), one can expect that variations in their magnitudes can change the corresponding positions of the Turing curves. For this purpose, we present two diagrams in $\epsilon - d_3$ parameter space for $d_2 = 0.2$ (Figure 14a) and $d_2 = 0.1$ (Figure 14b). The other parameter values used in preparation of these diagrams were $b_1 = 1.0$, $\alpha = 1.0$, $b_2 = 1.0$, $\beta = 1.5$, $\mu = 1.0$, $d = 0.5$, $b_3 = 1.0$ and $d_1 = 0.1$. In these two figures, the black vertical straight line (at $\epsilon = 5.486$) represents the temporal-Hopf bifurcation threshold. Since the expression for the temporal-Hopf bifurcation threshold did not involve any of the diffusion coefficients, we had the exact same position for this line in both the figures for the considered two different values of $d_2$. For the both of these figures, the blue curve represents the Turing curve for the local model (i.e., $\delta_1 = \delta_2 = 0$), magenta curve represents the Turing curve for the nonlocal model where the nonlocality is incorporated only in the first prey species (i.e., $\delta_2 = 0$) with $\delta_1 = 1.0$, green curve denotes the Turing curve for the nonlocal model where the nonlocality is incorporated only in the second prey species (i.e., $\delta_1 = 0$) with $\delta_2 = 1.0$, and red curve denotes the Turing curve for the nonlocal model (27)–(29) with $\delta_1 = \delta_2 = 1.0$. From these two diagrams (i.e., Figure 14a,b), we can observe that the lower value of the parameter $d_2$ enhanced the Turing instability region significantly. Also, these figures indicate that the introduction of nonlocal intra-specific competition for both the prey species led to the enlargement of the Turing instability region. For the chosen set of parameter values, we can easily notice that incorporation of nonlocality for the second prey species had greater impact on the enlargement of the Turing instability region than that of the first prey species.

<figure><img src="image_18.png" /><figcaption>**Figure 14.** Plots of the temporal-Hopf and Turing bifurcation curves in $(\epsilon, d_3)$-parameter space for (a) $d_2 = 0.2$, and (b) $d_2 = 0.1$. Here, black vertical straight line corresponds to the temporal-Hopf threshold, blue curve represents the Turing curve for $\delta_1 = \delta_2 = 0$, magenta curve represents the Turing curve for $\delta_1 = 1.0$ and $\delta_2 = 0$, green curve represents the Turing curve for $\delta_1 = 0$ and $\delta_2 = 1.0$, and red curve represents the Turing curve for $\delta_1 = \delta_2 = 1.0$. The black vertical straight line has been drawn by solving the equations presented in (9) simultaneously and the red curve has been drawn from the Equation (39). The magenta and green curves have been drawn from Equations (A4) and (A5), respectively. Brief derivations of Equations (A4) and (A5) are provided in Appendix B.</figcaption></figure>
For the nonlocal system (27)–(29), we used the central difference approximation for diffusion part and forward Euler scheme for the reaction part. The nonlocal interaction terms have been numerically approximated by trapezoidal rule. For all the patterns presented in this subsection, we employed the pulse-type initial conditions (26) and the periodic boundary conditions. Now, we are interested to investigate numerically the dependence of the number of patches with high population density for the stationary patterns on the extent of nonlocality inside the Turing domain. For this purpose, we considered the parameter values $d_2 = 0.2$, $d_3 = 9.0$, $\epsilon = 5.0$ and the remaining parameter values as the same as above. Also, we took equal extents of nonlocality for both the prey species (that is, $\delta_1 = \delta_2 = \delta$). We present the corresponding numerical results for the first prey species in Figure 15 over one-dimensional spatial domain [-50, 50]. We employed forward continuation technique in order to prepare this figure. This figure indicates the existence of different regimes with different number of patches depending upon the value of $\delta$. Interestingly, the number of patches initially increased with the increment in $\delta$ and then decreased.

<figure><img src="image_19.png" /><figcaption>**Figure 15.** Color plot of density of the first prey species ($u_1$) using the forward continuation technique. Here, $\epsilon = 5.0$, $d_3 = 9.0$ and other parameter values are mentioned in (25).</figcaption></figure>

Further, we are interested to explore the effect of the extent of nonlocality on the dynamic patterns exhibited by the corresponding local model. In order to do so, we first considered the parameter values for the homogeneous in space and periodic in time pattern presented in Figure 5a, and varied the extent of nonlocality ($\delta$). The resulting patterns for the first prey species are presented in Figure 16. The resulting patterns suggest that small values of $\delta$ retained the spatiotemporal distribution of the populations (e.g., Figure 16a for $\delta = 0.5$), a slight increment in $\delta$ altered the spatiotemporal distribution by giving rise to periodic in both the space and time population distribution (e.g., Figure 16b for $\delta = 0.7$), and further increments in $\delta$ resulted in the emergence of the stationary periodic in space population distribution (e.g., Figure 16c for $\delta = 0.8$). The dynamic nature of the patterns presented in Figure 16a,b can be easily understood from the phase portraits illustrated in Figure 17.

Now, we considered the parameter values for the chaotic pattern presented in Figure 8a, and varied the value of $\delta$ to observe the effect of it on the chaotic dynamics of the local model.
The corresponding numerical simulations suggest that the smaller values of $\delta$ were unable to alter the chaotic nature of the population distribution, however, interestingly destroyed the symmetry observed for the local model (e.g., Figure 18a for $\delta = 1.0$). A gradual increment in $\delta$ resulted in a chaotic pattern where long patches broke into smaller patches (e.g., Figure 18b for $\delta = 2.0$), and finally settled into a stationary periodic in space pattern for higher values of $\delta$ (e.g., Figure 18c for $\delta = 2.5$). Chaotic nature of the patterns presented in Figure 18a,b can be easily observed from the phase portraits illustrated in Figure 19. The chaotic nature of these patterns have also been cross verified by checking the sensitivity to initial conditions as described in [59], but we do not present the results here for the sake of brevity. For the intermediate periodic in both the space and time population distributions presented in Figures 5b and 8b, we found that the small extents of nonlocality for both the prey species could potentially stabilize the population distributions by giving rise to the stationary periodic in space distributions and we chose not to display these results here for the sake of brevity. Overall, the above observations indicate the stabilizing effect of the extents of nonlocality for all the three considered species.

<figure><img src="image_20.png" /><figcaption>**Figure 16.** Plots of different types of patterns exhibited by $u_1$ representing the effect of the extent of nonlocality for $\epsilon = 6.0$ and $d_3 = 3.0$. Three different patterns are obtained for (a) $\delta = 0.5$; (b) $\delta = 0.7$; (c) $\delta = 0.8$. Other parameter values are the same as that of Figure 5a.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 17.** Phase diagrams of spatial average densities of the associated three species corresponding to (a) Figure 16a, (b) Figure 16b.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 18.** Plots of different types of patterns exhibited by $u_1$ representing the effect of the extent of nonlocality for $\epsilon = 8.0$ and $d_3 = 1.0$. Three different patterns are obtained for (a) $\delta = 1.0$; (b) $\delta = 2.0$; (c) $\delta = 2.5$. Other parameter values are the same as that of Figure 8a.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 19.** Phase diagrams of spatial average densities of the associated three species corresponding to (a) Figure 18a, (b) Figure 18b.</figcaption></figure>

Finally, we want to explore the effects of the extent of nonlocality on the dynamics presented in Figures 1 and 11. For this purpose, we considered the same parameter values as it was used to prepare the bifurcation diagram presented in Figure 11 with $\delta_1 = \delta_2 = 2.0$. The resulting dynamics are encapsulated in Figure 20. In this figure, the red colored curve corresponds to the first prey component from the spatially homogeneous coexistence steady state and blue color stands for the simulation results of the nonlocal model (27)–(29). Here, we can observe that the average densities of the first prey species stayed almost below the steady state values throughout the region in presence of the nonlocal interactions. Also, note that spatial heterogeneity in population distribution appeared throughout the region as opposed to that for both the temporal and local spatiotemporal cases. This figure reconfirms that the spatially heterogeneous population distribution can suppress the spatially averaged density for the first prey species. The same scenario happens for the second prey species, whereas, the spatial heterogeneity enhances the spatially averaged density for the predator population. For the sake of brevity, we restrict ourselves from displaying the corresponding bifurcation diagrams for both the second prey and predator species.

<figure><img src="image_24.png" /><figcaption>**Figure 20.** Single parameter bifurcation diagram for the spatially averaged density of first prey population ($< u_1 >$) with respect to the parameter $\epsilon$ for the nonlocal model (27)–(29). Here, $\delta_1 = \delta_2 = 2.0$, $d_3 = 1.0$ and other parameter values are mentioned in (25).</figcaption></figure>
**5. Discussion**

The number of studies on spatiotemporal pattern formation for three-species models is comparatively less in existing literature and there exist even fewer works for models with nonlocal interactions. In this article, we have investigated for the possible impacts of random dispersal and nonlocal intra-specific interactions on the dynamics of a two-prey-one-predator system. For the temporal version of the model, we have mentioned the period-doubling route to chaos by means of numerical simulations. For the corresponding spatiotemporal models with both the local and nonlocal intra-specific interactions for the two considered prey species, we have derived the conditions for Turing instability through linear analysis. Further, we have carried out extensive numerical simulations in order to validate our theoretical predictions as well as to explore rich complex spatiotemporal dynamical behaviors which are beyond the scope of linear analysis. Overall, we have presented a comparative study on the dynamics possessed by the temporal, local, and nonlocal spatiotemporal models.

Local spatiotemporal model (11)–(13) along with the periodic boundary conditions possesses a wide variety of stationary and dynamic patterns. For the parameter values well inside the Turing domain, the model produces stationary periodic in space patterns for all the three considered species. However, an interesting phenomenon has been observed for these stationary patterns where the regions with high first prey density correspond to the same for predator but to low second prey density. This phenomenon arises due to the fact that the chosen values of the parameter $\epsilon$ are much large than the parameter $\mu$ and this makes the first prey species more favorable to predator than the second prey species. Apart from this, numerical simulations show the stabilizing effect of diffusion on the spatial population distributions by making the periodic in time homogeneous in space or even chaotic population distributions to stationary periodic in space one for higher values of diffusion coefficient of predator species ($d_3$). Both these transitions occur through the periodic in both space and time population distribution for intermediate values of $d_3$. We would like to mention here that the choice of other diffusion coefficients as bifurcation parameter will not reveal any additional feature as the main determining factor is the location of the parameter value with respect to the temporal-Hopf and Turing bifurcation curves.

Traditionally, the spatiotemporal dynamics of an ecological system have been investigated in the asymptotic sense. Sometimes it can take a very long time period for an ecological system to eventually settle down to a spatiotemporal state and this long time period can be as long as a few generations for certain species. Therefore, from the perspective of ecological time scales it is very important to investigate the quasi-stable transients for an ecological system and to identify the key factors for its appearance [60]. A long transient exhibits a quasi-stable dynamics for a few generations and then the ecological system experiences an abrupt transition to another regime [60]. Hence, it can provide an alternative explanation for ecological regime shifts apart from the usual concept of bifurcation [60]. However, bifurcation of dynamical behaviors for a system arises due to directional change in parameter values (for example, Figures 7 and 11) and this directional change is often considered to be mediated through external factors, such as climate change [60]. Therefore, bifurcation can provide useful explanation for exogenously triggered regime shifts whereas long transients can serve as alternative underlying mechanisms for regime shifts [60]. Hence, the bifurcation diagrams such as Figures 7 and 11 exhibit exogenously triggered regime shifts while the phase portraits presented in Figure 10 provide evidence for endogenously triggered regime shift for the local spatiotemporal model (11)–(13).

Further, we have investigated for the possible effects of random dispersal of all the considered species on the dynamics possessed by the corresponding temporal model. In this regard, we have observed that the local spatiotemporal model (11)–(13) with spatially homogeneous initial conditions preserves the temporal dynamics whereas the system with pulse type initial conditions produces smaller amplitude chaotic oscillations than the temporal chaos. This type of reduction in amplitude of chaotic oscillations arises due to the emergence of spatial heterogeneity in population distributions.
It is worthy to mention here that the spatially homogeneous perturbation over the entire domain is not ecologically realistic. Therefore, our results suggest that one should expect small amplitude chaotic oscillations for a spatially heterogeneous population as compared to a spatially homogeneous population. Also, we have shown that the backward continuation technique significantly enhances the spatially heterogeneous chaotic regime by destroying the period-doubling route to chaos. These bifurcation diagrams (that is, Figures 1, 11 and 12) vouch for the existence of alternative stable states for the spatiotemporal model without nonlocal interactions (11)–(13).

Another interesting dynamical aspect of the local spatiotemporal model (11)–(13) is that the nature of the patterns produced is not independent of the spatial domain size for certain parameter range, in particular, for parameter values in temporal-Hopf unstable region sufficiently far from the temporal-Hopf threshold. In this regard, we have shown that a chaotic pattern can become a two-periodic one in smaller spatial domain and even a stationary heterogeneous pattern can lose its stability through becoming a periodic one in both the space and time in smaller spatial domain. These spatial domain size dependence of spatiotemporal patterns have been reported in Figure 13. Therefore, our study suggests that one can get different spatiotemporal patterns depending on the spatial domain size. At this point, we would like to emphasize on the fact that this type of domain size dependence for the resulting spatiotemporal patterns is not unheard of at all. For instance, Barrio et al. [61] showed that both the domain size and shape can have significant roles in selection of different types of patterns. Also, Neville et al. [62] investigated how domain growth can influence pattern formation. Hence, our findings on domain size dependence of patterns are in accordance with literature and carry forward our understanding in this direction.

The dependence of results on the domain size can be explained as follows. If there exists a spatially periodic solution $u(x, t)$ with the wavenumber $k$ and the length of the interval $L_1 = \frac{2\pi n_1}{k}$, where $n_1 \in \mathbb{N}$, then this solution can be naturally extended by periodicity to the interval $L_2 = \frac{2\pi n_2}{k}$, where $n_2 > n_1$. Moreover, if we increase the length of the interval, solutions with other wavenumbers can appear due to their successive bifurcations from the spatially homogeneous stationary solutions. Furthermore, solutions with a fixed wavenumber persist under the variation of the length of the interval in some range of lengths. Hence, solutions with different wavenumbers can coexist for the same length, and their number increases with the increase of $L$. This multiplicity of solutions is conventionally observed in bifurcation problems. Their stability and basins of attraction depend on the length $L$. The interaction of different solutions can result in complex dynamics including chaotic behavior. The transition from two-periodic oscillations to chaotic oscillations for a larger length (cf. Figures 8 and 13) corresponds to this understanding of the dynamics of solutions. It is also possible that the initial condition remains in the basin of attraction of the same solution, such that we observe the same dynamics for a changing interval.

Finally, we have studied the effects of nonlocal interactions for both the prey species on the spatiotemporal dynamics. Our study suggests that the introduction of nonlocal interactions for both the prey species or any one of them can lead to the enlargement of the Turing instability region where nonlocality for the second prey species possess greater impact than that for the first prey species. In order to understand the impact of nonlocality on the stationary Turing patterns, we have employed forward continuation technique which shows the existence of multiple stationary regimes with different number of patches. The number of patches initially increases and then decreases with increasing extent of nonlocality. This scenario is encapsulated in Figure 15 and it is different from the result presented in [48] where the number of patches gradually increases. Another interesting feature of the dynamics of nonlocal model (27)–(29) is that it breaks the spatial symmetry possessed by local model (11)–(13) for chaotic population distribution. Our results also indicate that the nonlocality acts as a stabilizing factor as large extent of it can suppress the oscillatory behaviors by changing them to stationary patterns. Ecologically it can be interpreted as the access to the nearby resources allow the species to be confined within a patch rather than changing its position in a regular or irregular fashion. Further, our investigation suggests that the presence of spatial heterogeneity in
population distributions can diminish the spatially averaged density for both the prey species (for instance, see the Figure 20 for the first prey species) and enhances the same for the predator species.

This study leads to some other interesting and relevant problems to investigate. In this study, we have confined ourselves to only the top-hat kernel function in order to model the nonlocal interactions. However, it will be a fascinating future research prospect to investigate the spatiotemporal dynamics with other types of kernel functions such as Gaussian, Laplacian and triangular kernels etc. Another interesting problem will be to consider nonlocal consumption for predator species and examine its effect on dynamics. The present study has considered Lotka–Volterra type linear functional responses in order to capture the prey–predator interactions whereas the consideration of more suited functional responses such as Holling type-II and Beddington–DeAngelis functional responses will make the investigation ecologically more viable and challenging. These interesting as well as challenging issues will be taken care of in our future works.

**Author Contributions:** M.B. and V.V. proposed the formulation of the problem; K.M. performed the mathematical calculations and numerical simulations; M.B. analyzed and verified the results. All authors participated in the preparation of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** The first author gratefully acknowledges the financial support provided by Science and Engineering Research Board, Government of India for pursuing his post-doctoral research. V.V. was supported by the "RUDN University Program 5-100" and the French-Russian project PRC2307. M.B. was supported by SERB grant MTR/2018/000527.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A

The transcendental expression for $\mathcal{F}(k, \delta_1, \delta_2)$ introduced in Section 4 is given by

$$ \begin{aligned} \mathcal{F}(k, \delta_1, \delta_2) = & 4d_1d_2(d_1\mu^2u_2^* + d_2\epsilon^2u_1^*)k^9 + \{d_1^2\mu^2u_2^{*2}\cos(k\delta_2) + d_2^2\epsilon^2u_1^{*2}\cos(k\delta_1) - 6d_1d_2(\alpha + \beta)\epsilon\mu u_1^*u_2^*\}k^7 \\ & + \{8d_1d_2u_1^*u_2^*\left(\epsilon^2\frac{\sin(k\delta_2)}{\delta_2} + \mu^2\frac{\sin(k\delta_1)}{\delta_1}\right) + d_1^2\mu^2u_2^{*2}\frac{\sin(k\delta_2)}{\delta_2} + d_2^2\epsilon^2u_1^{*2}\frac{\sin(k\delta_1)}{\delta_1}\}k^6 \\ & - \{(\alpha + \beta)\epsilon\mu u_1^*u_2^*(d_1u_2^*\cos(k\delta_2) + d_2u_1^*\cos(k\delta_1))\}k^5 \\ & + u_1^*u_2^*\left\{2\epsilon^2\left(d_1u_2^*\frac{\sin(k\delta_2)\cos(k\delta_2)}{\delta_2} + d_2u_1^*\frac{\sin(k\delta_1)\cos(k\delta_2)}{\delta_1} + d_1u_2^*\frac{\sin(k\delta_1)\cos(k\delta_2)}{\delta_1}\right) + 2\mu^2\left(d_2u_1^*\frac{\sin(k\delta_1)\cos(k\delta_1)}{\delta_1} + d_1u_2^*\frac{\sin(k\delta_2)\cos(k\delta_2)}{\delta_2}\right)\right\}k^4 \\ & + u_1^*u_2^*\left\{4\left(d_1\epsilon^2u_2^*\frac{\sin^2(k\delta_2)}{\delta_2^2} + d_2\mu^2u_1^*\frac{\sin^2(k\delta_1)}{\delta_1^2}\right) + \alpha\beta u_1^*u_2^*(\epsilon^2\cos(k\delta_2) + \mu^2\cos(k\delta_1)) + 2(d_1\mu^2u_2^* + d_2\epsilon^2u_1^*)\frac{\sin(k\delta_1)\sin(k\delta_2)}{\delta_1\delta_2} + 2\alpha\beta(\alpha + \beta)\epsilon\mu u_1^*u_2^*\right\}k^3 \\ & - u_1^*u_2^*\left\{3\alpha\beta\left(\epsilon^2\frac{\sin(k\delta_2)}{\delta_2} + \mu^2\frac{\sin(k\delta_1)}{\delta_1}\right) + (\alpha + \beta)\epsilon\mu\left(\frac{\cos(k\delta_1)\sin(k\delta_2)}{\delta_2} + \frac{\sin(k\delta_1)\cos(k\delta_2)}{\delta_1}\right)\right\}k^2 \\ & + u_1^*u_2^*\left\{\epsilon^2\frac{\cos(k\delta_1)\sin^2(k\delta_2)}{\delta_2^2} + \mu^2\frac{\sin^2(k\delta_1)\cos(k\delta_2)}{\delta_1^2}\right\}k + u_1^*u_2^*\left\{\left(\epsilon^2\frac{\sin(k\delta_2)}{\delta_2} + \mu^2\frac{\sin(k\delta_1)}{\delta_1}\right)\frac{\sin(k\delta_1)\sin(k\delta_2)}{\delta_1\delta_2}\right\}. \end{aligned} \quad (A1) $$

## Appendix B

Here, we briefly describe about the calculation of critical value $d_3^T$ for inducing Turing instability with nonlocal intra-specific competition being present in only one prey species.

Case 1: $\delta_1 \neq 0, \delta_2 = 0$

Proceeding in a similar manner as described in Section 4, we have

$$ \tilde{C}(k, \delta_1) = d_1d_2d_3k^6 + d_3\left\{d_1u_2^* + d_2u_1^*\frac{\sin(k\delta_1)}{k\delta_1}\right\}k^4 + \left\{d(d_1\mu^2u_2^* + d_2\epsilon^2u_1^*)v^* + d_3\left(\frac{\sin(k\delta_1)}{k\delta_1} - \alpha\beta\right)u_1^*u_2^*\right\}k^2 \\ + d\left\{(\epsilon^2 + \mu^2\frac{\sin(k\delta_1)}{k\delta_1}) - (\alpha + \beta)\epsilon\mu\right\}u_1^*u_2^*v^*. \quad (A2) $$
Then, the Turing bifurcation threshold can be obtained by solving the following two equations simultaneously

$$ \tilde{C}(k, \delta_1) = 0, \quad \frac{\partial}{\partial k} \tilde{C}(k, \delta_1) = 0. \qquad (A3) $$

Solving the first equation of (A3) for $d_3$, we obtain

$$ d_3 = - \frac{d \left\{ (d_1 \mu^2 u_2^* + d_2 \epsilon^2 u_1^*) k^2 + \left( \epsilon^2 + \mu^2 \frac{\sin(k \delta_1)}{k \delta_1} - (\alpha + \beta) \epsilon \mu \right) u_1^* u_2^* \right\} v^*}{d_1 d_2 k^6 + \left\{ d_1 u_2^* + d_2 u_1^* \frac{\sin(k \delta_1)}{k \delta_1} \right\} k^4 + \left\{ \frac{\sin(k \delta_1)}{k \delta_1} - \alpha \beta \right\} u_1^* u_2^* k^2}. \qquad (A4) $$

Now, substituting the above expression for $d_3$ in the second equation of (A3) and solving for $k$, we obtain the critical wavenumber $k_T$. Further, substitution of $k = k_T$ into the Equation (A4) gives us the desired critical value $d_3^T$.

Case 2: $\delta_1 = 0, \delta_2 \neq 0$

Proceeding in a similar manner as described above, we obtain the following critical value $d_3^T$ to induce Turing instability in this case:

$$ d_3^T = - \frac{d \left\{ (d_1 \mu^2 u_2^* + d_2 \epsilon^2 u_1^*) k_T^2 + \left( \epsilon^2 \frac{\sin(k_T \delta_2)}{k_T \delta_2} + \mu^2 - (\alpha + \beta) \epsilon \mu \right) u_1^* u_2^* \right\} v^*}{d_1 d_2 k_T^6 + \left\{ d_1 u_2^* \frac{\sin(k_T \delta_2)}{k_T \delta_2} + d_2 u_1^* \right\} k_T^4 + \left\{ \frac{\sin(k_T \delta_2)}{k_T \delta_2} - \alpha \beta \right\} u_1^* u_2^* k_T^2}. \qquad (A5) $$

## References

1. Paine, R.T. Food web complexity and species diversity. Am. Nat. **1966**, 100, 65–75. [CrossRef]

2. Fujii, K. Complexity-stability relationship of two-prey-one-predator species system model: Local and global stability. J. Theor. Biol. **1977**, 69, 613–623. [CrossRef]

3. Parrish, J.D.; Saila, S.B. Interspecific competition, predation and species diversity. J. Theor. Biol. **1970**, 27, 207–220. [CrossRef]

4. Vance, R.R. Predation and resource partitioning in one predator-two prey model communities. Am. Nat. **1978**, 112, 797–813. [CrossRef]

5. Takeuchi, Y.; Adachi, N. Existence and bifurcation of stable equilibrium in two-prey, one-predator communities. Bull. Math. Biol. **1983**, 45, 877–900. [CrossRef]

6. Hutson, V.; Vickers, G.T. A criterion for permanent coexistence of species, with an application to a two-prey one-predator system. Math. Biosci. **1983**, 63, 253–269. [CrossRef]

7. Matsuda, H.; Kawasaki, K.; Shigesada, N.; Teramoto, E.; Ricciardi, L.M. Switching effect on the stability of the prey-predator system with three trophic levels. J. Theor. Biol. **1986**, 122, 251–262. [CrossRef]

8. Abrams, P.; Matsuda, H. Effects of adaptive predatory and anti-predator behaviour in a two-prey-one-predator system. Evol. Ecol. **1993**, 7, 312–326. [CrossRef]

9. Klebanoff, A.; Hastings, A. Chaos in one-predator, two-prey models: General results from bifurcation theory. Math. Biosci. **1994**, 122, 221–233. [CrossRef]

10. Abrams, P.A. Is predator-mediated coexistence possible in unstable systems? Ecology **1999**, 80, 608–621. [CrossRef]

11. van Leeuwen, E.; Jansen, V.A.A.; Bright, P.W. How population dynamics shape the functional response in a one-predator-two-prey system. Ecology **2007**, 88, 1571–1581. [CrossRef] [PubMed]

12. Lee, S.S.; Kajiwara, T. The effect of the remains of the carcass in a two-prey, one-predator model. Discret. Contin. Dyn. Syst. Ser. B **2008**, 9, 353–374.

13. Zaman, G.; Saker, S.H. Dynamics and control of a system of two non-interacting preys with common predator. Math. Methods Appl. Sci. **2011**, 34, 2259–2273. [CrossRef]
14. Groll, F.; Arndt, H.; Altland, A. Chaotic attractor in two-prey one-preator system originates from interplay of limit cycles. *Theor. Ecol.* **2016**, *10*, 147–154. [CrossRef]

15. Turchin, P.; Ellner, S.P. Living on the edge of chaos: Population dynamics of Fennoscandian voles. *Ecology* **2000**, *81*, 3099–3116. [CrossRef]

16. Levin, S.A.; Segel, L.A. Hypothesis for origin of planktonic patchiness. *Nature* **1976**, *259*, 659. [CrossRef]

17. Turing, A.M. The chemical basis of morphogenesis. *Philos. Trans. R. Soc. Lond. B Biol. Sci.* **1952**, *237*, 37–72.

18. Manna, K.; Banerjee, M. Stationary, non-stationary and invasive patterns for a prey-predator system with additive Allee effect in prey growth. *Ecol. Complex.* **2018**, *36*, 206–217. [CrossRef]

19. Baurmann, M.; Gross, T.; Feudel, U. Instabilities in spatially extended predator-prey systems: Spatio-temporal patterns in the neighborhood of Turing-Hopf bifurcations. *J. Theor. Biol.* **2007**, *245*, 220–229. [CrossRef]

20. Zhang, L.; Liu, J.; Banerjee, M. Hopf and steady state bifurcation analysis in a ratio-dependent predator-prey model. *Commun. Nonlinear Sci. Numer. Simul.* **2017**, *44*, 52–73. [CrossRef]

21. White, K.A.J.; Gilligan, C.A. Spatial heterogeneity in three-species, plant-parasite-hyperparasite, systems. *Phil. Trans. R. Soc. Lond. B* **1998**, *353*, 543–557. [CrossRef]

22. Petrovskii, S.; Kawasaki, K.; Takasu, F.; Shigesada, N. Diffusive waves, dynamical stabilization and spatio-temporal chaos in a community of three competitive species. *Jpn. J. Ind. Appl. Math.* **2001**, *18*, 459–481. [CrossRef]

23. Chen, W.; Peng, R. Stationary patterns created by cross-diffusion for the competitor-competitor-mutualist model. *J. Math. Anal. Appl.* **2004**, *291*, 550–564. [CrossRef]

24. Pang, P.Y.H.; Wang, M. Strategy and stationary pattern in a three-species predator-prey model. *J. Differ. Equ.* **2004**, *200*, 245–273. [CrossRef]

25. Wang, M. Stationary patterns for a prey-predator model with prey-dependent and ratio-dependent functional responses and diffusion. *Phys. D* **2004**, *196*, 172–192. [CrossRef]

26. Wang, M. Stationary patterns caused by cross-diffusion for a three-species prey-predator model. *Comput. Math. Appl.* **2006**, *52*, 707–720. [CrossRef]

27. Tian, C. Turing patterns created by cross-diffusion for a Holling II and Leslie-Gower type three species food chain model. *J. Math. Chem.* **2011**, *49*, 1128–1150. [CrossRef]

28. Tian, C.; Ling, Z.; Lin, Z. Turing pattern formation in a predator-prey-mutualist system. *Nonlinear Anal. Real World Appl.* **2011**, *12*, 3224–3237. [CrossRef]

29. Rao, F. Spatiotemporal complexity of a three-species ratio-dependent food chain model. *Nonlinear Dyn.* **2014**, *76*, 1661–1676. [CrossRef]

30. Abid, W.; Yafia, R.; Alaoui, M.A.A.; Bouhafa, H.; Abichou, A. Instability and pattern formation in three-species food chain model via Holling type II functional response on a circular domain. *Int. J. Bifurc. Chaos* **2015**, *25*, 1550092. [CrossRef]

31. Ikeda, T.; Mimura, M. An interfacial approach to regional segregation of two competing species mediated by a predator. *J. Math. Biol.* **1993**, *31*, 215–240. [CrossRef]

32. Ko, W.; Ryu, K. Analysis of diffusive two-competing-prey and one-predator systems with Beddington-DeAngelis functional response. *Nonlinear Anal. Theory Methods Appl.* **2009**, *71*, 4185–4202. [CrossRef]

33. Lv, Y.; Yuan, R.; Pei, Y. Turing pattern formation in a three species model with generalist predator and cross-difusion. *Nonlinear Anal. Theory, Methods Appl.* **2013**, *85*, 214–232. [CrossRef]

34. Xie, Z. Cross-diffusion induced Turing instability for a three species food chain model. *J. Math. Anal. Appl.* **2012**, *388*, 539–547. [CrossRef]

35. Peng, R.; Shi, J.; Wang, M. Stationary pattern of a ratio-dependent food chain model with diffusion. *SIAM J. Appl. Math.* **2007**, *67*, 1479–1503. [CrossRef]

36. Feng, W. Coexistence, stability, and limiting behavior in a one-predator-two-prey model. *J. Math. Anal. Appl.* **1993**, *179*, 592–609. [CrossRef]

37. Mukherjee, N.; Ghorai, S.; Banerjee, M. Detection of Turing patterns in a three species food chain model via amplitude equation. *Commun. Nonlinear Sci. Numer. Simul.* **2019**, *69*, 219–236. [CrossRef]

38. Tanzy, M.C.; Volpert, V.A.; Bayliss, A.; Nehrkorn, M.E. Stability and pattern formation for competing populations with asymmetric nonlocal coupling. *Math. Biosci.* **2013**, *246*, 14–26. [CrossRef]
39. Autry, E.A.; Bayliss, A.; Volpert, V.A. Biological control with nonlocal interactions. *Math. Biosci.* **2018**, *301*, 129–146. [CrossRef]

40. Pal, S.; Ghorai, S.; Banerjee, M. Effect of kernels on spatio-temporal patterns of a non-local prey-predator model. *Math. Biosci.* **2019**, *310*, 96–107. [CrossRef]

41. Merchant, S.M.; Nagata, W. Instabilities and spatiotemporal patterns behind predator invasions with nonlocal prey competition. *Theor. Popul. Biol.* **2011**, *80*, 289–297. [CrossRef] [PubMed]

42. Segal, B.L.; Volpert, V.A.; Bayliss, A. Pattern formation in a model of competing populations with nonlocal interactions. *Phys. D* **2013**, *253*, 12–22. [CrossRef]

43. Gourley, S.A. Travelling front solutions of a nonlocal Fisher equation. *J. Math. Biol.* **2000**, *41*, 272–284. [CrossRef]

44. Merchant, S.M.; Nagata, W. Selection and stability of wave trains behind predator invasions in a model with non-local prey competition. *IMA J. Appl. Math.* **2015**, *80*, 1155–1177. [CrossRef]

45. Deng, K.; Wu, Y. Global stability for a nonlocal reaction-diffusion population model. *Nonlinear Anal. Real World Appl.* **2015**, *25*, 127–136. [CrossRef]

46. Banerjee, M.; Volpert, V. Prey-predator model with a nonlocal consumption of prey. *Chaos* **2016**, *26*, 083120. [CrossRef] [PubMed]

47. Banerjee, M.; Volpert, V. Spatio-temporal pattern formation in Rosenzweig-MacArthur model: Effect of nonlocal interactions. *Ecol. Complex.* **2017**, *30*, 2–10. [CrossRef]

48. Pal, S.; Ghorai, S.; Banerjee, M. Analysis of a prey-predator model with non-local interaction in the prey population. *Bull. Math. Biol.* **2018**, *80*, 906–925. [CrossRef]

49. Bayliss, A.; Volpert, V.A. Complex predator invasion waves in a Holling-Tanner model with nonlocal prey interaction. *Phys. D* **2017**, *346*, 37–58. [CrossRef]

50. Han, B.-S.; Yang, Y.-H. On a predator-prey reaction-diffusion model with nonlocal effects. *Commun. Nonlinear Sci. Numer. Simul.* **2017**, *46*, 49–61. [CrossRef]

51. Tian, C.; Ling, Z.; Zhang, L. Nonlocal interaction driven pattern formation in a prey-predator model. *Appl. Math. Comput.* **2017**, *308*, 73–83. [CrossRef]

52. Banerjee, M.; Mukherjee, N.; Volpert, V. Prey-predator model with a nonlocal bistable dynamics of prey. *Mathematics* **2018**, *6*, 41. [CrossRef]

53. Pal, S.; Banerjee, M.; Ghorai, S. Spatio-temporal pattern formation in Holling-Tanner type model with nonlocal consumption of resources. *Int. J. Bifurc. Chaos* **2019**, *29*, 1930002. [CrossRef]

54. Ni, W.; Shi, J.; Wang, M. Global stability and pattern formation in a nonlocal diffusive Lotka-Volterra competition model. *J. Differ. Equ.* **2018**, *264*, 6891–6932. [CrossRef]

55. Han, R.; Dai, B.; Chen, Y. Pattern formation in a diffusive intraguild predation model with nonlocal interaction effects. *AIP Adv.* **2019**, *9*, 035046. [CrossRef]

56. Takeuchi, Y. *Global Dynamical Properties of Lotka-Volterra Systems*; World Scientific: Singapore, 1996.

57. Othmer, H.G.; Scriven, L.E. Interactions of reaction and diffusion in open systems. *Ind. Eng. Chem. Fundamen.* **1969**, *8*, 302–313. [CrossRef]

58. Murray, J.D. *Mathematical Biology*; Springer: Heidelberg/Berlin, Germany, 1989.

59. Banerjee, M.; Abbas, S. Existence and non-existence of spatial patterns in a ratio-dependent predator-prey model. *Ecol. Complex.* **2015**, *21*, 199–214. [CrossRef]

60. Hastings, A.; Abbott, K.C.; Cuddington, K.; Francis, T.; Gellner, G.; Lai, Y.-C.; Morozov, A.; Petrovskii, S.; Scranton, K.; Zeeman, M.L. Transient phenomena in ecology. *Science* **2018**, *361*, eaat6412. [CrossRef]

61. Barrio, R.A.; Varea, C.; Aragón, J.L.; Maini, P.K. A two-dimensional numerical study of spatial pattern formation in interacting Turing systems. *Bull. Math. Biol.* **1999**, *61*, 483–505. [CrossRef]

62. Neville, A.A.; Matthews, P.C.; Byrne, H.M. Interactions between pattern formation and domain growth. *Bull. Math. Biol.* **2006**, *68*, 1975–2003. [CrossRef]

<figure><img src="image_25.png" /></figure>