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

#### **Kalyan Manna 1, Vitaly Volpert 2,3,4 and Malay Banerjee 1,\***


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

**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). 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 *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 *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} = \begin{array}{c c} u\_1 \left( b\_1 - u\_1 - \alpha u\_2 - \epsilon v \right), \end{array} \tag{1}$$

$$\frac{du\_2}{dt} = \left(u\_2\left(b\_2 - \beta u\_1 - u\_2 - \mu v\right)\right) \tag{2}$$

$$\frac{dv}{dt} \quad = \ \ v \left( -b\_3 + d\epsilon u\_1 + d\mu u\_2 \right) \,, \tag{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 *bi* (*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 *α* and *β*, while the parameters  and *μ* 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*<sup>000</sup> = (0, 0, 0), *E*+<sup>00</sup> = (*b*1, 0, 0), *E*0+<sup>0</sup> = (0, *b*2, 0), *E*++<sup>0</sup> = *b*1−*αb*<sup>2</sup> <sup>1</sup>−*αβ* , *<sup>b</sup>*2−*βb*<sup>1</sup> <sup>1</sup>−*αβ* , 0 , *E*+0<sup>+</sup> = *<sup>b</sup>*<sup>3</sup> *<sup>d</sup>*  , 0, *db*<sup>1</sup><sup>−</sup>*b*<sup>3</sup> *d* <sup>2</sup> and *E*0++ = 0, *<sup>b</sup>*<sup>3</sup> *<sup>d</sup><sup>μ</sup>* , *db*2*μ*−*b*<sup>3</sup> *dμ*<sup>2</sup> . Also, this system admits a unique coexistence equilibrium point *E*+++ = (*u*∗ <sup>1</sup>, *u*<sup>∗</sup> <sup>2</sup>, *v*∗), where

$$\begin{array}{rcl} u\_1^+ &=& \frac{b\_3 \epsilon - db\_2 \epsilon \mu - ab\_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(a\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{array}$$

The feasibility conditions for the coexistence equilibrium point *E*+++ are given by (i) *b*<sup>3</sup> + *db*1*μ*<sup>2</sup> > *μ*(*db*<sup>2</sup> + *αb*3), *b*3*μ* + *db*<sup>2</sup><sup>2</sup> > (*db*1*μ* + *βb*3), *b*3*αβ* + *d*(*b*<sup>1</sup> + *b*2*μ*) > *b*<sup>3</sup> + *d*(*b*2*α*  + *b*1*βμ*) and <sup>2</sup> + *μ*<sup>2</sup> > (*α* + *β*) *μ*, or (ii) *b*<sup>3</sup> + *db*1*μ*<sup>2</sup> < *μ*(*db*<sup>2</sup> + *αb*3), *b*3*μ* + *db*<sup>2</sup><sup>2</sup> < (*db*1*μ* + *βb*3), *b*3*αβ* + *d*(*b*<sup>1</sup> + *b*2*μ*) < *b*<sup>3</sup> + *d*(*b*2*α*  + *b*1*βμ*) and <sup>2</sup> + *μ*<sup>2</sup> < (*α* + *β*) *μ*.

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 \quad = \quad 0,\tag{4}
$$

where

$$
\mu\_2 \quad = \quad \mu\_1^\* + \mu\_2^\*.\tag{5}
$$

$$a\_1 \quad = \quad (1 - a\beta)u\_1^\* u\_2^\* + d(\epsilon^2 u\_1^\* + \mu^2 u\_2^\*)v^\*,\tag{6}$$

*<sup>a</sup>*<sup>0</sup> <sup>=</sup> *<sup>d</sup>*{(<sup>2</sup> <sup>+</sup> *<sup>μ</sup>*2) <sup>−</sup> (*<sup>α</sup>* <sup>+</sup> *<sup>β</sup>*) *μ*}*u*<sup>∗</sup> 1*u*<sup>∗</sup> <sup>2</sup> *v*∗. (7)

Therefore, the equilibrium point *E*+++ is locally asymptotically stable provided the following Routh–Hurwitz criteria are satisfied:

$$a\_2 > 0, \ a\_0 > 0, \ a\_2 a\_1 - a\_0 > 0. \tag{8}$$

The condition *a*<sup>2</sup> > 0 is automatically satisfied for feasible *E*+++. However, we need to have the parametric restriction <sup>2</sup> + *μ*<sup>2</sup> > (*α* + *β*) *μ* in order to satisfy the condition *a*<sup>0</sup> > 0. Therefore, the parametric restriction <sup>2</sup> + *μ*<sup>2</sup> > (*α* + *β*) *μ* acts as a necessary condition for local stability. Furthermore, the model (1)–(3) undergoes Hopf bifurcation at  = <sup>∗</sup> when the following conditions are satisfied:

$$
\langle (a\_2 a\_1 - a\_0) |\_{\mathfrak{c} = \mathfrak{c}\_\*} = 0, \ \ \frac{d}{d\epsilon} (a\_2 a\_1 - a\_0)|\_{\mathfrak{c} = \mathfrak{c}\_\*} \neq 0,\tag{9}
$$

with the parametric restriction <sup>2</sup> + *μ*<sup>2</sup> > (*α* + *β*) *μ*.

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, \; a = 1.0, \; b\_2 = 1.0, \; \beta = 1.5, \; \mu = 1.0, \; b\_3 = 1.0, \; d = 0.5,\tag{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 1.** Single parameter bifurcation diagram for the density of first prey population (*u*1) with respect to the parameter  for the temporal model (1)–(3). Other parameter values are mentioned in (10).

#### **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} = \begin{aligned} d\_1 \frac{\partial^2 u\_1}{\partial x^2} + u\_1 \left( b\_1 - u\_1 - \alpha u\_2 - \epsilon v \right), \end{aligned} \tag{11}$$

$$\frac{\partial \mu\_2}{\partial t} = -d\_2 \frac{\partial^2 \mu\_2}{\partial x^2} + \mu\_2 \left( b\_2 - \beta \mu\_1 - \mu\_2 - \mu v \right), \tag{12}$$

$$\frac{\partial v}{\partial t} \quad = \quad d\_3 \frac{\partial^2 v}{\partial x^2} + v \left( -b\_3 + d \varepsilon u\_1 + d \mu u\_2 \right) \,, \tag{13}$$

subjected to non-negative initial conditions and periodic boundary conditions. Here, *u*1(*x*, *t*), *u*2(*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*, *<sup>L</sup>*], where *<sup>L</sup>*(<sup>&</sup>gt; <sup>0</sup>) <sup>∈</sup> <sup>R</sup>. The positive parameters *di* (*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*∗ <sup>1</sup>, *u*<sup>∗</sup> <sup>2</sup>, *v*∗) and then linearizing we find the following Jacobian matrix:

$$\mathcal{J}(k^2) \quad = \begin{bmatrix} -u\_1^\* - d\_1 k^2 & -\alpha u\_1^\* & -\varepsilon u\_1^\* \\ -\beta u\_2^\* & -u\_2^\* - d\_2 k^2 & -\mu u\_2^\* \\ d\varepsilon v^\* & d\mu v^\* & -d\_3 k^2 \end{bmatrix}, \tag{14}$$

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

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

where

$$\begin{array}{rcl} A(k^2) &=& (d\_1 + d\_2 + d\_3)k^2 + u\_1^\* + u\_2^\*, \\ B(k^2) &=& (d\_1 d\_2 + d\_2 d\_3 + d\_3 d\_1)k^4 + \{d\_1 u\_2^\* + d\_2 u\_1^\* + d\_3 (u\_1^\* + u\_2^\*)\} k^2 + (1 - a\beta) u\_1^\* u\_2^\* + d(\varepsilon^2 u\_1^\* + \mu^2 u\_2^\*) v^\*, \\ C(k^2) &=& d\_1 d\_2 d\_3 \bar{k}^6 + d\_3 (d\_1 u\_2^\* + d\_2 u\_1^\*) k^4 + \{d\_1 \mu^2 u\_2^\* v^\* + d\_2 \mu \varepsilon^2 u\_1^\* v^\* + d\_3 (1 - a\beta) u\_1^\* u\_2^\*\} k^2 + \{d(\varepsilon^2 + \mu^2) u\_1^\* u\_2^\* + d(\varepsilon^2 + \mu^2) u\_1^\* u\_2^\*\} k^4 + \{d\_1 \mu^2 u\_1^\* + d\_2 \mu \varepsilon^2 u\_1^\*\} k^6 + \{d\_1 \mu^2 u\_2^\* + d\_2 \mu \varepsilon^2 u\_1^\*\} k^2, \\ & & - d(a + \beta) \varepsilon \mu \} \, u\_1^\* u\_2^\* v^\*. \end{array}$$

Now, Re(*λ*) 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) - \mathbb{C}(k^2) > 0. \tag{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 *λ*1, *λ*<sup>2</sup> and *λ*<sup>3</sup> 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 = \quad -A(k^2),
\tag{17}
$$

$$
\lambda\_1 \lambda\_2 + \lambda\_2 \lambda\_3 + \lambda\_3 \lambda\_1 \quad = \quad B(k^2), \tag{18}
$$

$$
\lambda\_1 \lambda\_2 \lambda\_3 \quad = \quad -\mathbb{C}(k^2), \tag{19}
$$

$$-(\lambda\_1 + \lambda\_2)(\lambda\_2 + \lambda\_3)(\lambda\_3 + \lambda\_1) \quad = \quad A(k^2)B(k^2) - \mathbb{C}(k^2). \tag{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 *kT*, whereas other two eigenvalues still have negative real parts [58]. Since at the critical wavenumber *k* = *kT* we have one of the roots is equal to zero, without any loss of generality we assume

$$
\lambda\_1 \mid\_{k^2 = k\_T^2} = 0, \text{ } \text{Re}(\lambda\_2) \mid\_{k^2 = k\_T^2} < 0 \text{ and } \text{ } \text{Re}(\lambda\_3) \mid\_{k^2 = k\_T^2} < 0. \tag{21}
$$

Hence, at the critical wavenumber *k* = *kT* we obtain *C*(*k*<sup>2</sup> *<sup>T</sup>*) = 0. Further due to the above conditions (21) for Turing instability, we get *A*(*k*<sup>2</sup> *<sup>T</sup>*) > 0, *<sup>B</sup>*(*k*<sup>2</sup> *<sup>T</sup>*) > 0 and *<sup>A</sup>*(*k*<sup>2</sup> *T*)*B*(*k*<sup>2</sup> *<sup>T</sup>*) − *<sup>C</sup>*(*k*<sup>2</sup> *<sup>T</sup>*) = *A*(*k*<sup>2</sup> *T*)*B*(*k*<sup>2</sup> *<sup>T</sup>*) > 0. Thus, the system remains stable if *<sup>C</sup>*(*k*2) > 0 holds for all *<sup>k</sup>* 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 *kT* beyond the Turing bifurcation threshold. The expression of *C*(*k*2) can be rewritten as

$$\mathbb{C}(k^2) \quad \equiv \quad \mathbb{C}\_3(k^2)^3 + \mathbb{C}\_2(k^2)^2 + \mathbb{C}\_1(k^2) + \mathbb{C}\_{0\prime} \tag{22}$$

where *C*<sup>3</sup> > 0 since the diffusion coefficients are positive constants and *C*<sup>0</sup> > 0 from the local asymptotic stability condition of *E*+++. Also, the coefficient *C*<sup>2</sup> > 0 as the expression of it is given by *C*<sup>2</sup> = *d*3(*d*1*u*<sup>∗</sup> <sup>2</sup> + *d*2*u*<sup>∗</sup> <sup>1</sup> ). The minimum for *<sup>C</sup>*(*k*2) occurs at *<sup>k</sup>* = *kT* where

$$k\_T^2 \quad = \frac{-\mathcal{C}\_2 + \sqrt{\mathcal{C}\_2^2 - 3\mathcal{C}\_1\mathcal{C}\_3}}{3\mathcal{C}\_3},\tag{23}$$

as we have *dC*(*k*2) *<sup>d</sup>*(*k*2) <sup>=</sup> 0 and *<sup>d</sup>*2*C*(*k*2) *<sup>d</sup>*(*k*2)<sup>2</sup> <sup>&</sup>gt; 0 at *<sup>k</sup>* <sup>=</sup> *kT*. Now from the expression (23), we can notice that *<sup>k</sup>*<sup>2</sup> *T* is positive if *C*<sup>1</sup> < 0 or *C*<sup>2</sup> < 0 and *C*<sup>2</sup> <sup>2</sup> > 3*C*1*C*3. Since *C*<sup>2</sup> > 0 holds for our considered system, then the condition for the positivity of *k*<sup>2</sup> *<sup>T</sup>* reduces to *C*<sup>1</sup> < 0. Therefore, the Turing bifurcation boundary is given by

$$2\mathbf{C}\_2^3 - 9\mathbf{C}\_1\mathbf{C}\_2\mathbf{C}\_3 - 2(\mathbf{C}\_2^2 - 3\mathbf{C}\_1\mathbf{C}\_3)^{\frac{3}{2}} + 2\mathbf{7}\mathbf{C}\_0\mathbf{C}\_3^2 \quad = \quad 0. \tag{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, \; a = 1.0, \; b\_2 = 1.0, \; b = 1.5, \; \mu = 1.0, \; d = 0.5, \; b\_3 = 1.0, \; d\_1 = 0.1, \; d\_2 = 0.2,\tag{25}$$

and varied the parameters  and *d*3. For this set of parameter values, the temporal-Hopf bifurcation threshold is approximately given by  = 5.486. First, we present the temporal-Hopf and Turing bifurcation curves for the local spatiotemporal model (11)–(13) in (, *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 2.** Plots of temporal-Hopf and Turing bifurcation curves in (, *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).

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{array}{rclcrcl}\mu\_{1}(\mathbf{x},0) &=& \left\{ \begin{array}{rcl} u\_{1}^{\*} + \eta\_{1'} & |\mathbf{x}| < 10, \\ u\_{1'}^{\*} & & 10 \leq |\mathbf{x}| \leq L, \end{array} \right. \\\ \mu\_{2}(\mathbf{x},0) &=& \left\{ \begin{array}{rcl} u\_{2}^{\*} + \eta\_{2'} & |\mathbf{x}| < 10, \\ u\_{2'}^{\*} & & 10 \leq |\mathbf{x}| \leq L, \end{array} \right. \\\ v(\mathbf{x},0) &=& \left\{ \begin{array}{rcl} v^{\*} + \eta\_{3'} & |\mathbf{x}| < 10, \\ v^{\*} & & 10 \leq |\mathbf{x}| \leq L, \end{array} \right. \end{array} \tag{26}$$

where |*ηr*| 1, for *r* = 1, 2, 3.

Now, we consider a point (, *d*3)=(5.4, 9.0), which lies in the Turing region. For this choice of , 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 3.** (**a**) Plot of stationary Turing pattern exhibited by *u*1, (**b**) time evolution of the average densities of three species.  = 5.4 and *d*<sup>3</sup> = 9.0, other parameter values are mentioned in (25).

Now if we choose  = 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*<sup>3</sup> = 3.0 in Figure 5a, periodic in both space and time pattern for *d*<sup>3</sup> = 3.75 in Figure 5b, and periodic in space but stationary in time pattern for *d*<sup>3</sup> = 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*<sup>3</sup> for  = 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*<sup>3</sup> = 3.67. The periodic in both space and time pattern exists for *d*<sup>3</sup> ∈ [3.68, 3.91] and then stationary pattern emerges for *d*<sup>3</sup> > 3.91. Interestingly, the number of patches with high first prey density gradually decreases with the increment in the value of *d*3.

**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].

**Figure 5.** Plots of different types of patterns exhibited by *u*<sup>1</sup> for  = 6.0. Three different patterns are obtained for (**a**) *d*<sup>3</sup> = 3.0; (**b**) *d*<sup>3</sup> = 3.75; (**c**) *d*<sup>3</sup> = 9.0. Other parameter values are mentioned in (25).

**Figure 6.** Phase diagrams of spatial average densities of the associated three species corresponding to (**a**) Figure 5a, (**b**) Figure 5b.

**Figure 7.** Single parameter bifurcation diagram for the spatially averaged density of first prey population (< *u*<sup>1</sup> >) with respect to the diffusion coefficient *d*3. Here,  = 6.0 and other parameter values are mentioned in (25).

If we further increase the magnitude of  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*<sup>3</sup> = 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*<sup>3</sup> (for example, *d*<sup>3</sup> = 10.0). This transition from spatiotemporal chaotic to stationary pattern takes place via a pattern where patches move periodically for intermediate values of *d*<sup>3</sup> (for example, *d*<sup>3</sup> = 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.

(**a**) (**b**) (**c**)

**Figure 8.** Plots of different types of patterns exhibited by *u*<sup>1</sup> for  = 8.0. Three different patterns are obtained for (**a**) *d*<sup>3</sup> = 1.0; (**b**) *d*<sup>3</sup> = 9.5; (**c**) *d*<sup>3</sup> = 10.0. Other parameter values are mentioned in (25).

**Figure 9.** Phase diagrams of spatial average densities of the associated three species corresponding to (**a**) Figure 8a, (**b**) Figure 8b.

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  = 8.0 and *d*<sup>3</sup> = 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* ∈ [18, 000, 20, 000] and *t* ∈ [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 10.** Phase diagrams of spatial average densities of the associated three species for  = 8.0 and *d*<sup>3</sup> = 9.0 with (**a**) *t* ∈ [18, 000, 20, 000], (**b**) *t* ∈ [58, 000, 60, 000]. Other parameter values are mentioned in (25).

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*<sup>1</sup> = 0.1, *d*<sup>2</sup> = 0.2 and *d*<sup>3</sup> = 1.0, and varied the value of the parameter . 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 -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  = 8.5 approximately. For higher values of  (that is,  ∈ (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:


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 -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 11.** Single parameter bifurcation diagram for the spatially averaged density of first prey population (< *u*<sup>1</sup> >) with respect to the parameter  for the spatiotemporal model (11)–(13). Here, we have used pulse type initial population distributions for each value of . Other parameter values are mentioned in (25) with *d*<sup>3</sup> = 1.0.

**Figure 12.** Single parameter bifurcation diagram for the spatially averaged density of first prey population (< *u*<sup>1</sup> >) with respect to the parameter  for the spatiotemporal model (11)–(13). Here, backward continuation technique has been employed. Other parameter values are mentioned in (25) with *d*<sup>3</sup> = 1.0.

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 13.** (**a**) Plot of spatial pattern exhibited by *u*1, (**b**) corresponding phase portrait of the average densities of three species for  = 8.0, *d*<sup>3</sup> = 1.0. (**c**) Plot of spatial pattern exhibited by *u*1, (**d**) corresponding phase portrait of the average densities of three species for  = 8.0, *d*<sup>3</sup> = 10.0. Spatial domain size is [−50, 50] and other parameter values are mentioned in (25).

#### **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} \quad = \quad d\_1 \frac{\partial^2 u\_1}{\partial x^2} + u\_1 \left( b\_1 - f(u\_1) - au\_2 - \epsilon v \right), \tag{27}$$

$$\frac{\partial \mathfrak{u}\_2}{\partial t} = -d\_2 \frac{\partial^2 \mathfrak{u}\_2}{\partial \mathbf{x}^2} + \mathfrak{u}\_2 \left( b\_2 - \beta \mathfrak{u}\_1 - J(\mathfrak{u}\_2) - \mu \upsilon \right), \tag{28}$$

$$\frac{\partial v}{\partial t} \quad = \quad d\mathfrak{z} \frac{\partial^2 v}{\partial \mathbf{x}^2} + v \left( -b\mathfrak{z} + d\varepsilon u\_1 + d\mu u\_2 \right), \tag{29}$$

where

$$J(u\_r) \quad = \int\_{-\infty}^{\infty} u\_r(y, t) \phi\_r(x - y) dy, \ r = 1, 2, \ \lambda$$

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 *φ<sup>r</sup>* (*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 *φ<sup>r</sup>* (*r* = 1, 2) as even functions. Throughout this paper, we restrict ourselves only to the top-hat kernel functions defined as follows:

$$\phi\_r(\mathbf{x}) \;=\; \begin{cases} \frac{1}{2\delta\_r} & |\mathbf{x}| \le \delta\_r \\ 0, & \text{otherwise}, \end{cases} \tag{30}$$

for *r* = 1, 2. The newly introduced parameters *δ<sup>r</sup>* (*r* = 1, 2) represent the extent of the nonlocality by controlling the width of the kernel functions *φ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*∗ <sup>1</sup>, *u*<sup>∗</sup> <sup>2</sup>, *v*∗), we introduce the perturbation *u*<sup>1</sup> = *u*<sup>∗</sup> <sup>1</sup> <sup>+</sup> *<sup>ε</sup><sup>u</sup>*<sup>1</sup>*eλt*+*ikx*, *u*<sup>2</sup> = *u*<sup>∗</sup> <sup>2</sup> <sup>+</sup> *<sup>ε</sup><sup>u</sup>*<sup>2</sup>*eλt*+*ikx* and *<sup>v</sup>* <sup>=</sup> *<sup>v</sup>*<sup>∗</sup> <sup>+</sup> *<sup>ε</sup>ve <sup>λ</sup>t*+*ikx* with <sup>|</sup> *<sup>ε</sup>* | 1. Then the Jacobian matrix corresponding to the linearized system is given by

$$
\begin{split}
\widetilde{\mathcal{J}}(k,\delta\_{1},\delta\_{2}) &= \begin{bmatrix}
d\varepsilon v^{\*} & d\mu v^{\*} & -d\_{3}k^{2}
\end{bmatrix} \\ &= \begin{bmatrix}
d\varepsilon v^{\*} & d\mu v^{\*} & -d\_{3}k^{2}
\end{bmatrix}.
\end{split}
\tag{31}
$$

Therefore, the characteristic equation is given by

$$\left(\lambda^3 + \widetilde{A}(k, \delta\_1, \delta\_2)\lambda^2 + \widetilde{B}(k, \delta\_1, \delta\_2)\lambda + \widetilde{C}(k, \delta\_1, \delta\_2)\right) \\
= \begin{array}{c} 0, \\ \end{array} \tag{32}$$

where

*<sup>A</sup>*(*k*, *<sup>δ</sup>*1, *<sup>δ</sup>*2)=(*d*<sup>1</sup> <sup>+</sup> *<sup>d</sup>*<sup>2</sup> <sup>+</sup> *<sup>d</sup>*3)*k*<sup>2</sup> <sup>+</sup> *<sup>u</sup>*<sup>∗</sup> 1 sin(*kδ*1) *kδ*<sup>1</sup> + *u*∗ 2 sin(*kδ*2) *kδ*<sup>2</sup> , (33) *<sup>B</sup>*(*k*, *<sup>δ</sup>*1, *<sup>δ</sup>*2) = *<sup>d</sup>*3*k*<sup>2</sup> 0 (*d*<sup>1</sup> + *d*2)*k*<sup>2</sup> + *u*<sup>∗</sup> 1 sin(*kδ*1) *kδ*<sup>1</sup> + *u*∗ 2 sin(*kδ*2) *kδ*<sup>2</sup> 1 + *d*1*k*<sup>2</sup> + *u*<sup>∗</sup> 1 sin(*kδ*1) *kδ*<sup>1</sup> *d*2*k*<sup>2</sup> <sup>+</sup> *<sup>u</sup>*<sup>∗</sup> 2 sin(*kδ*2) *kδ*<sup>2</sup> +*d*(<sup>2</sup>*u*<sup>∗</sup> <sup>1</sup> + *<sup>μ</sup>*2*u*<sup>∗</sup> <sup>2</sup> )*v*<sup>∗</sup> − *αβu*<sup>∗</sup> 1*u*<sup>∗</sup> 2 = (*d*1*d*<sup>2</sup> + *d*2*d*<sup>3</sup> + *d*3*d*1)*k*<sup>4</sup> + 0 *d*1*u*<sup>∗</sup> 2 sin(*kδ*2) *kδ*<sup>2</sup> + *d*2*u*<sup>∗</sup> 1 sin(*kδ*1) *kδ*<sup>1</sup> + *d*<sup>3</sup> *u*∗ 1 sin(*kδ*1) *kδ*<sup>1</sup> + *u*∗ 2 sin(*kδ*2) *kδ*<sup>2</sup> 1 *<sup>k</sup>*<sup>2</sup> + sin(*kδ*1) sin(*kδ*2) *k*2*δ*1*δ*<sup>2</sup> <sup>−</sup> *αβ u*∗ 1*u*<sup>∗</sup> <sup>2</sup> + *d*(<sup>2</sup>*u*<sup>∗</sup> <sup>1</sup> + *<sup>μ</sup>*2*u*<sup>∗</sup> <sup>2</sup> )*v*∗, (34) *<sup>C</sup>*(*k*, *<sup>δ</sup>*1, *<sup>δ</sup>*2) = *<sup>d</sup>*3*k*<sup>2</sup> *d*1*k*<sup>2</sup> + *u*<sup>∗</sup> 1 sin(*kδ*1) *kδ*<sup>1</sup> *d*2*k*<sup>2</sup> <sup>+</sup> *<sup>u</sup>*<sup>∗</sup> 2 sin(*kδ*2) *kδ*<sup>2</sup> + *d* <sup>2</sup>*u*<sup>∗</sup> <sup>1</sup> *v*<sup>∗</sup> *d*2*k*<sup>2</sup> + *u*<sup>∗</sup> 2 sin(*kδ*2) *kδ*<sup>2</sup> + *dμ*2*u*<sup>∗</sup> <sup>2</sup> *v*<sup>∗</sup> *d*1*k*<sup>2</sup> + *u*<sup>∗</sup> 1 sin(*kδ*1) *kδ*<sup>1</sup> − *d*3*αβu*<sup>∗</sup> 1*u*<sup>∗</sup> <sup>2</sup> *<sup>k</sup>*<sup>2</sup> <sup>−</sup> *<sup>d</sup>*(*<sup>α</sup>* <sup>+</sup> *<sup>β</sup>*) *μu*<sup>∗</sup> 1*u*<sup>∗</sup> <sup>2</sup> *v*<sup>∗</sup> = *d*1*d*2*d*3*k*<sup>6</sup> + *d*<sup>3</sup> 0 *d*1*u*<sup>∗</sup> 2 sin(*kδ*2) *kδ*<sup>2</sup> + *d*2*u*<sup>∗</sup> 1 sin(*kδ*1) *kδ*<sup>1</sup> 1 *k*<sup>4</sup> + 0 *d*1*dμ*2*u*<sup>∗</sup> <sup>2</sup> *v*<sup>∗</sup> + *d*2*d* <sup>2</sup>*u*<sup>∗</sup> <sup>1</sup> *v*<sup>∗</sup> + *d*<sup>3</sup> sin(*kδ*1) sin(*kδ*2) *k*2*δ*1*δ*<sup>2</sup> <sup>−</sup> *αβ u*∗ 1*u*<sup>∗</sup> 2 1 *k*2 +*d* 0<sup>2</sup> sin(*kδ*2) *kδ*<sup>2</sup> <sup>+</sup> *<sup>μ</sup>*<sup>2</sup> sin(*kδ*1) *kδ*<sup>1</sup> <sup>−</sup> (*<sup>α</sup>* <sup>+</sup> *<sup>β</sup>*) *μ*<sup>1</sup> *u*∗ 1*u*<sup>∗</sup> <sup>2</sup> *v*∗. (35)

Using the Routh–Hurwitz criteria, we can say that the spatially homogeneous coexistence steady state *E*+++ = (*u*∗ <sup>1</sup>, *u*<sup>∗</sup> <sup>2</sup>, *v*∗) will be asymptotically stable if and only if

$$
\tilde{A}(k,\delta\_1,\delta\_2) > 0,\\
\tilde{C}(k,\delta\_1,\delta\_2) > 0,\\
\tilde{A}(k,\delta\_1,\delta\_2)\tilde{B}(k,\delta\_1,\delta\_2) - \tilde{C}(k,\delta\_1,\delta\_2) > 0. \tag{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:

$$
\widetilde{A}(k,\delta\_1,\delta\_2) > 0,\\
\widetilde{B}(k,\delta\_1,\delta\_2) > 0,\\
\widetilde{C}(k,\delta\_1,\delta\_2) = 0. \tag{37}
$$

Note that the expressions of *A*, *B* and *C* are transcendental in nature of their arguments *k*, *δ*<sup>1</sup> and *δ*2. The associated Turing bifurcation threshold can be obtained by solving the following two equations simultaneously

$$
\widetilde{\mathcal{C}}(k,\delta\_1,\delta\_2) = 0,\ \frac{\partial}{\partial k} \widetilde{\mathcal{C}}(k,\delta\_1,\delta\_2) = 0.\tag{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}\varepsilon^{2}u\_{1}^{\*})k^{2} + \left(\varepsilon^{2}\frac{\sin(k\delta\_{2})}{k\delta\_{2}} + \mu^{2}\frac{\sin(k\delta\_{1})}{k\delta\_{1}} - (\mathfrak{a} + \beta)\varepsilon\mu\right)u\_{1}^{\*}u\_{2}^{\*}\right\}v^{\*}}{d\_{1}d\_{2}k^{6} + \left\{d\_{1}u\_{2}^{\*}\frac{\sin(k\delta\_{2})}{\delta\_{2}} + d\_{2}u\_{1}^{\*}\frac{\sin(k\delta\_{1})}{\delta\_{1}}\right\}k^{3} + \left\{\frac{\sin(k\delta\_{1})\sin(k\delta\_{2})}{k^{2}\delta\_{1}\delta\_{2}} - a\beta\right\}u\_{1}^{\*}u\_{2}^{\*}k^{2}}} . \tag{39}$$

Now, differentiating *C*(*k*, *δ*1, *δ*2) with respect to *k*, we get

$$\begin{array}{rcl} \frac{\partial \hat{\phi}^{\circ}}{\partial \mathbf{k}} &=& d\_{3} \left[ \left\{ d\_{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})}{\delta\_{2}} + d\_{2} u\_{1}^{\*} \frac{\sin(k \delta\_{1})}{\delta\_{1}} \right\} k^{2} - \right] \\ &2 d\beta u\_{1}^{\*} u\_{2}^{\*} k + \left\{ \frac{\cos(k \delta\_{1}) \sin(k \delta\_{2})}{\delta\_{2}} + \frac{\sin(k \delta\_{1}) \cos(k \delta\_{2})}{\delta\_{1}} \right\} u\_{1}^{\*} u\_{2}^{\*} \right] + d \left[ 2 (d\_{1} \mu^{2} u\_{2}^{\*} + d\_{2} \varepsilon^{2} u\_{1}^{\*}) k + \\ &\left\{ \varepsilon^{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] \mathbf{v}^{\*}. \end{array} \tag{40}$$

Substituting the expression for *d*<sup>3</sup> in the second equation of (38), we obtain

$$\mathcal{F}(k, \delta\_1, \delta\_2) \quad = \quad 0,\tag{41}$$

where the expression for F, which is transcendental in nature, is given in Appendix A. Solving the Equation (41) for *k*, we obtain the critical wavenumber *kT*. But the transcendental nature of the Equation (41) prevents us from finding an analytical solution for *kT* and hence, we find the critical value numerically. Further, substituting the value of *kT* into the Equation (39) we obtain the critical value of the diffusion coefficient of predator *d<sup>T</sup>* <sup>3</sup> 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*<sup>2</sup> 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  − *d*<sup>3</sup> parameter space for *d*<sup>2</sup> = 0.2 (Figure 14a) and *d*<sup>2</sup> = 0.1 (Figure 14b). The other parameter values used in preparation of these diagrams were *b*<sup>1</sup> = 1.0, *α* = 1.0, *b*<sup>2</sup> = 1.0, *β* = 1.5, *μ* = 1.0, *d* = 0.5, *b*<sup>3</sup> = 1.0 and *d*<sup>1</sup> = 0.1. In these two figures, the black vertical straight line (at  = 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., *δ*<sup>1</sup> = *δ*<sup>2</sup> = 0), magenta curve represents the Turing curve for the nonlocal model where the nonlocality is incorporated only in the first prey species (i.e., *δ*<sup>2</sup> = 0) with *δ*<sup>1</sup> = 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., *δ*<sup>1</sup> = 0) with *δ*<sup>2</sup> = 1.0, and red curve denotes the Turing curve for the nonlocal model (27)–(29) with *δ*<sup>1</sup> = *δ*<sup>2</sup> = 1.0. From these two diagrams (i.e., Figure 14a,b), we can observe that the lower value of the parameter *d*<sup>2</sup> 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 14.** Plots of the temporal-Hopf and Turing bifurcation curves in (, *d*3)-parameter space for (**a**) *d*<sup>2</sup> = 0.2, and (**b**) *d*<sup>2</sup> = 0.1. Here, black vertical straight line corresponds to the temporal-Hopf threshold, blue curve represents the Turing curve for *δ*<sup>1</sup> = *δ*<sup>2</sup> = 0, magenta curve represents the Turing curve for *δ*<sup>1</sup> = 1.0 and *δ*<sup>2</sup> = 0, green curve represents the Turing curve for *δ*<sup>1</sup> = 0 and *δ*<sup>2</sup> = 1.0, and red curve represents the Turing curve for *δ*<sup>1</sup> = *δ*<sup>2</sup> = 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.

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*<sup>2</sup> = 0.2, *d*<sup>3</sup> = 9.0,  = 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, *δ*<sup>1</sup> = *δ*<sup>2</sup> = *δ*). 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 *δ*. Interestingly, the number of patches initially increased with the increment in *δ* and then decreased.

**Figure 15.** Color plot of density of the first prey species (*u*1) using the forward continuation technique. Here,  = 5.0, *d*<sup>3</sup> = 9.0 and other parameter values are mentioned in (25).

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 (*δ*). The resulting patterns for the first prey species are presented in Figure 16. The resulting patterns suggest that small values of *δ* retained the spatiotemporal distribution of the populations (e.g., Figure 16a for *δ* = 0.5), a slight increment in *δ* altered the spatiotemporal distribution by giving rise to periodic in both the space and time population distribution (e.g., Figure 16b for *δ* = 0.7), and further increments in *δ* resulted in the emergence of the stationary periodic in space population distribution (e.g., Figure 16c for *δ* = 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 *δ* to observe the effect of it on the chaotic dynamics of the local model. The corresponding numerical simulations suggest that the smaller values of *δ* 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 *δ* = 1.0). A gradual increment in *δ* resulted in a chaotic pattern where long patches broke into smaller patches (e.g., Figure 18b for *δ* = 2.0), and finally settled into a stationary periodic in space pattern for higher values of *δ* (e.g., Figure 18c for *δ* = 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 16.** Plots of different types of patterns exhibited by *u*<sup>1</sup> representing the effect of the extent of nonlocality for  = 6.0 and *d*<sup>3</sup> = 3.0. Three different patterns are obtained for (**a**) *δ* = 0.5; (**b**) *δ* = 0.7; (**c**) *δ* = 0.8. Other parameter values are the same as that of Figure 5a.

**Figure 17.** Phase diagrams of spatial average densities of the associated three species corresponding to (**a**) Figure 16a, (**b**) Figure 16b.

**Figure 18.** Plots of different types of patterns exhibited by *u*<sup>1</sup> representing the effect of the extent of nonlocality for  = 8.0 and *d*<sup>3</sup> = 1.0. Three different patterns are obtained for (**a**) *δ* = 1.0; (**b**) *δ* = 2.0; (**c**) *δ* = 2.5. Other parameter values are the same as that of Figure 8a.

**Figure 19.** Phase diagrams of spatial average densities of the associated three species corresponding to (**a**) Figure 18a, (**b**) Figure 18b.

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 *δ*<sup>1</sup> = *δ*<sup>2</sup> = 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 20.** Single parameter bifurcation diagram for the spatially averaged density of first prey population (< *u*<sup>1</sup> >) with respect to the parameter  for the nonlocal model (27)–(29). Here, *δ*<sup>1</sup> = *δ*<sup>2</sup> = 2.0, *d*<sup>3</sup> = 1.0 and other parameter values are mentioned in (25).

#### **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  are much large than the parameter *μ* 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 loose 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*<sup>1</sup> = <sup>2</sup>*πn*<sup>1</sup> *<sup>k</sup>* , where *<sup>n</sup>*<sup>1</sup> <sup>∈</sup> <sup>N</sup>, then this solution can be naturally extended by periodicity to the interval *<sup>L</sup>*<sup>2</sup> <sup>=</sup> <sup>2</sup>*πn*<sup>2</sup> *<sup>k</sup>* , where *n*<sup>2</sup> > *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 F(*k*, *δ*1, *δ*2) introduced in Section 4 is given by

$$\begin{cases} \begin{aligned} &f(\delta,\delta\_{1})=2\\ &4d\_{1}(d\_{1}d\_{1}u^{2}\_{2}+d\_{2}^{2}u\_{1}^{4})k^{8}+\left(d\_{1}^{2}l^{2}u\_{2}^{2}\cos(k\delta\_{2})+d\_{2}^{2}cu\_{1}^{2}\cos(k\delta\_{1})-6d\_{1}d\_{2}(a+\beta)\epsilon\mu u\_{1}^{u}u\_{2}^{\*}\right)k^{7}+\left(\beta d\_{1}d\_{2}u\_{1}^{2}\cos(k\delta\_{2})+d\_{1}^{2}\mu\epsilon\mu u\_{1}^{u}u\_{2}^{\*}\right)k^{6}+\left(\beta d\_{1}d\_{2}\epsilon\mu u\_{1}^{2}\right)k^{3}+\left(\left(d\_{1}+\beta\right)\epsilon\mu u\_{1}^{u}u\_{2}^{2}\right)k^{2}+\left(\left(d\_{1}+\beta\right)\epsilon\mu u\_{1}^{u}u\_{2}^{2}\right)k^{3}+\left(\left(d\_{1}+\beta\right)\epsilon\mu u\_{1}^{u}u\_{2}^{2}\right)k^{2}+\left(\left(d\_{1}+\beta\right)\epsilon\mu u\_{1}^{u}u\_{2}\right)k^{2}+\left(\beta d\_{1}\right)\epsilon\mu u\_{1}^{u}u\_{2}^{2}\left(d\_{1}+\beta\right)\epsilon\mu u\_{1}^{u}u\_{2}^{2}+\left(d\_{1}+\beta\right)\epsilon\mu u\_{1}^{u}u\_{2}^{2}+\left(d\_{1}+\beta\right)\epsilon\mu u\_{1}^{u}u\_{2}^{2}+\left(d\_{1}+\beta\right)\epsilon\mu u\_{1}^{u}u\_{2}^{2}+\left(d\_{1}+\beta\right)\epsilon\mu u\_{1}^{u}u\_{2}^{2}+\left(d\_{1}+\beta\right)\epsilon\mu u\_{1}^{u}u\_{2}^{2}+\left(d\_$$

#### **Appendix B**

Here, we briefly describe about the calculation of critical value *d<sup>T</sup>* <sup>3</sup> for inducing Turing instability with nonlocal intra-specific competition being present in only one prey species.

*Case* 1*: δ*<sup>1</sup> = 0*, δ*<sup>2</sup> = 0

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

$$\begin{array}{rcl}\check{C}(k,\delta\_{1})&=&d\_{1}d\_{2}d\_{3}k^{6}+d\_{3}\left\{d\_{1}u\_{2}^{\*}+d\_{2}u\_{1}^{\*}\frac{\sin(k\delta\_{1})}{k\delta\_{1}}\right\}k^{4}+\left\{d(d\_{1}\mu^{2}u\_{2}^{\*}+d\_{2}\epsilon^{2}u\_{1}^{\*})v^{\*}+d\_{3}\left(\frac{\sin(k\delta\_{1})}{k\delta\_{1}}-a\beta\right)u\_{1}^{\*}u\_{2}^{\*}\right\}k^{2}\\&+d\left\{\left(\epsilon^{2}+\mu^{2}\frac{\sin(k\delta\_{1})}{k\delta\_{1}}\right)-(a+\beta)\epsilon\mu\right\}u\_{1}^{\*}u\_{2}^{\*}v^{\*}.\end{array} \tag{A2}$$

Then, the Turing bifurcation threshold can be obtained by solving the following two equations simultaneously

$$
\tilde{\mathcal{C}}(k,\delta\_1) = 0,\\
\frac{\partial}{\partial k} \tilde{\mathcal{C}}(k,\delta\_1) = 0. \tag{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}\varepsilon^{2}u\_{1}^{\*})k^{2} + \left(\varepsilon^{2} + \mu^{2}\frac{\sin(k\delta\_{1})}{k\delta\_{1}} - (\mathfrak{a} + \beta)\varepsilon\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}} - a\beta\right\}u\_{1}^{\*}u\_{2}^{\*}k^{2}}.\tag{A4}$$

Now, substituting the above expression for *d*<sup>3</sup> in the second equation of (A3) and solving for *k*, we obtain the critical wavenumber *kT*. Further, substitution of *k* = *kT* into the Equation (A4) gives us the desired critical value *d<sup>T</sup>* 3 .

*Case* 2*: δ*<sup>1</sup> = 0*, δ*<sup>2</sup> = 0

Proceeding in a similar manner as described above, we obtain the following critical value *d<sup>T</sup>* <sup>3</sup> to induce Turing instability in this case:

$$d\_3^\Gamma = -\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 - (a + \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} - a\beta \right\} u\_1^\* u\_2^\* k\_T^2}. \tag{A5}$$

#### **References**


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

### *Article* **Persistence for a Two-Stage Reaction-Diffusion System**

**Robert Stephen Cantrell 1,\*, Chris Cosner <sup>1</sup> and Salomé Martínez <sup>2</sup>**


Received: 3 February 2020; Accepted: 5 March 2020; Published: 11 March 2020

**Abstract:** In this article, we study how the rates of diffusion in a reaction-diffusion model for a stage structured population in a heterogeneous environment affect the model's predictions of persistence or extinction for the population. In the case of a population without stage structure, faster diffusion is typically detrimental. In contrast to that, we find that, in a stage structured population, it can be either detrimental or helpful. If the regions where adults can reproduce are the same as those where juveniles can mature, typically slower diffusion will be favored, but if those regions are separated, then faster diffusion may be favored. Our analysis consists primarily of estimates of principal eigenvalues of the linearized system around (0, 0) and results on their asymptotic behavior for large or small diffusion rates. The model we study is not in general a cooperative system, but if adults only compete with other adults and juveniles with other juveniles, then it is. In that case, the general theory of cooperative systems implies that, when the model predicts persistence, it has a unique positive equilibrium. We derive some results on the asymptotic behavior of the positive equilibrium for small diffusion and for large adult reproductive rates in that case.

**Keywords:** reaction-diffusion; spatial ecology; population dynamics; stage structure; dispersal

**AMS Classification:** 92D40; 92D50; 35P15; 35K57

#### **1. Introduction**

The question of how dispersal interacts with spatial heterogeneity to influence population dynamics and species interactions has been studied extensively in recent years, specifically from the viewpoint of reaction-diffusion systems and related models—see, for example, [1–3] and the references cited therein. Most work on that topic assumes that each population is structured only by space and has only one mode of dispersal. However, populations are often structured by age, stage, or other attributes, and there may be variation among individuals in their dispersal rates or patterns. Here we will examine how the presence of a stage structure influences how diffusion rates influence population dynamics in a class of reaction-diffusion models for a population with two stages. In the case of a population with logistic growth, without age or stage structure, diffusing in a closed bounded spatially heterogeneous environment that is constant in time, it is well known that reaction-diffusion models predict that slower diffusion rates are advantageous relative to faster diffusion—see [4,5]. The results in [5] also hold for patch models. More broadly, a wide class of models arising in population genetics, population dynamics, and related areas display some version of the reduction principle, which says that dispersal, which causes faster mixing, typically reduces the rate of population growth—see [6]. However, the situation seems to be quite different in the case of stage structured populations. In [7], the authors considered a discrete-time patch model for a structured population and found that, in some cases, there was no selection against faster dispersal. The goal of the present paper is to use a spatially explicit reaction-diffusion model to understand how the spatial

distributions of habitats that are favorable for reproduction by adults and those that are favorable for survival and growth by juveniles affect whether faster diffusion is advantageous or harmful for a stage structured population. We will see that the answer depends on the details of the spatial distribution of favorable and unfavorable habitats.

The type of reaction-diffusion model we will study is

$$\begin{cases} \frac{\partial u}{\partial t} = d\_1 \Delta u + r(\mathbf{x})v - s(\mathbf{x})u - a(\mathbf{x})u - b(\mathbf{x})u^2 - c(\mathbf{x})uv \text{ in } \Omega, \ t > 0, \\\\ \frac{\partial v}{\partial t} = d\_2 \Delta v + s(\mathbf{x})u - e(\mathbf{x})v - f(\mathbf{x})v^2 - g(\mathbf{x})uv \text{ in } \Omega, \ t > 0, \\\\ \nabla u \cdot \nu = \nabla v \cdot \nu = 0 \text{ on } \partial\Omega, \ t > 0. \end{cases} \tag{1}$$

Ω is a bounded domain in R*N*, and *ν* is the outward unit normal to *∂*Ω, such that the system has Neumann boundary conditions, which are the no-flux boundary conditions for simple diffusion. In this system, *u* and *v* represent the population densities of juveniles and individuals that have reached reproductive age, i.e. adults, respectively, of the same species. Thus, the term *s*(*x*) represents the rate at which juveniles mature into adults, which is determined by the fraction of individuals that reach reproductive age and the rate at which they mature, while *r*(*x*) accounts for the local fecundity of adults such that *r*(*x*)*v*(*x*) describes that rate at which new juveniles are produced by an adult population with density *v* at location *x*. The terms *a*(*x*), *b*(*x*), *c*(*x*), *e*(*x*), *f*(*x*), and *g*(*x*) account for per-capita death rates and saturation factors due to logistic self-limitation. The diffusion coefficients *d*<sup>1</sup> and *d*<sup>2</sup> account for the the dispersal rates of juveniles and adults, respectively. The coefficients are all assumed to be nonnegative and continuous in Ω. This is the type of model for a stage structured population introduced in [8]. Related models with a different interpretation are discussed in [9,10] and in the references in those papers. The model expressed in Equation (1) is not an explicitly age-structured model. It assumes that individuals in the juvenile stage mature at some spatially dependent rate but does not track the age of individuals within each stage. Explicitly age structured models are considered in [11–13]. A different way of modeling an age-structured population, based on delayed reaction diffusion equations, is developed in [14]. Our focus here is on how spatial heterogeneity, dispersal, and stage structure interact, so we chose to use the simplest possible formulation of stage structure. In the case where *c* = *g* = 0, the system is cooperative, and the methods and results of [15,16] would apply to it. The linearization of Equation (1) around (0, 0) is cooperative, so the results of [15] apply to it; in particular, with a few technical assumptions, they imply that it has a principal eigenvalue.

The main questions we will address in this work are related to understanding the roles of the different functions and coefficients in Equation (1) in the persistence of the species. For the remainder of the paper we will focus primarily on understanding how the principal eigenvalue of the linearization of Equation (1) around (0, 0) depends on the coefficients and what that dependence means biologically. We will see that whether faster diffusion is harmful or helpful for the persistence of the population depends on the details of the distribution of habitats that are favorable for adult reproduction and those that are favorable for juvenile survival and maturation. In some cases slower diffusion is still an advantage, but sometimes faster diffusion turns out to be helpful, and sufficiently fast diffusion may even be necessary for persistence. The spatial distribution of habitats favorable to adult reproduction (*r*(*x*) large) relative to those favorable to juvenile development (*s*(*x*) large) turns out to be important in some cases. Our analysis here is similar in spirit to the sorts of results obtained for diffusive Lotka-Volterra competition models in [3,17–19]. In particular, we will examine the behavior of the system for small, large, and general diffusion rates. Related results for some epidemiological models are derived in [20,21].

The linearization of Equation (1) around (0, 0) has a principal eigenvalue whose sign determines whether the model predicts persistence or extinction. Since the sign of the principal eigenvalue of the linearization of Equation (1) around (0, 0) determines the fate that Equation (1) predicts for the population it describes, we will study in detail the following problem:

$$\begin{cases} \begin{array}{rcl} d\_1 \Delta \boldsymbol{\uprho} + r(\boldsymbol{x}) \boldsymbol{\uppsi} - (s(\boldsymbol{x}) + a(\boldsymbol{x})) \boldsymbol{\uprho} &= \boldsymbol{\uplambda} \boldsymbol{\uprho} \text{ in } \Omega, \\\ d\_2 \Delta \boldsymbol{\uppsi} + s(\boldsymbol{x}) \boldsymbol{\uprho} - \boldsymbol{\upepsilon}(\boldsymbol{x}) \boldsymbol{\uppsi} &= \boldsymbol{\uplambda} \boldsymbol{\uppsi} \text{ in } \Omega, \\\ \nabla \boldsymbol{\uprho} \cdot \boldsymbol{\upnu} = \nabla \boldsymbol{\uppsi} \cdot \boldsymbol{\upnu} &= \boldsymbol{0} \text{ on } \partial \Omega. \end{array} \end{cases} \tag{2}$$

#### **2. Basic Properties**

In this section, we discuss some basic properties of Equation (1). From now on, we assume that *<sup>r</sup>*,*s*, *<sup>a</sup>*, *<sup>b</sup>*, *<sup>c</sup>*,*e*, *<sup>f</sup>* , *<sup>g</sup>* <sup>∈</sup> *<sup>C</sup>α*(Ω¯ ), *<sup>∂</sup>*<sup>Ω</sup> is of class *<sup>C</sup>*2,*α*, and that the following hypotheses hold:

**Hypothesis 1** (**H1**)**.** *r*(*x*),*s*(*x*) ≥ 0 *in* Ω*, with r*(*xr*) = 0*, s*(*xs*) = 0 *for some xr*, *xs* ∈ Ω*.*

**Hypothesis 2** (**H2**)**.** *a*(*x*), *c*(*x*),*e*(*x*), *g*(*x*) ≥ 0 *in* Ω*.*

**Hypothesis 3** (**H3**)**.** *b*(*x*) > 0*, f*(*x*) > 0 *in* Ω*.*

The model expressed in Equation (1) has many mathematical features in common with the models discussed in [9,10] for populations where individuals can switch between two different movement modes. A key feature is that the linear part of Equation (1) is cooperative, so it will have a principal eigenvalue which determines the stability of the equilibrium (0, 0) and hence the persistence or extinction of the population. Another key feature of Equation (1) is that the nonlinearity is subhomogeneous. The maximum principle and existence of principal eigenvalues for cooperative linear systems such as the linear part of the right side of Equation (1) are derived in [15]. The general theory for systems such as Equation (1) is developed in [16] for the fully cooperative case (where *c* = *g* = 0, such that adults only compete with other adults and juveniles with other juveniles) and in the general case in [8–10,22]. As expected, the sign of the principal eigenvalue of the linearization of Equation (1) around (0, 0) gives us the relevant information to study the persistence of the species. If it is positive, the population will persist. If it is nonpositive, the population will go extinct. In the case where the coefficients of *c* and *g* in Equation (1) are zero such that the system is cooperative, the results and methods of [16] imply that, if the principal eigenvalue of the linear part is positive, then the system has a unique globally attractive equilibrium. If those coefficients are small, the methods of [9,10] can be applied to show that Equation (1) is asymptotically cooperative, and still has a unique globally attractive positive equilibrium. Combining results that are given in [8–10,15,16] or that follow directly by the same arguments used in those papers, we have the following:

**Lemma 1.** *The eigenvalue problem expressed in Equation (2) has a unique principal eigenvalue λ*<sup>1</sup> *that is characterized by having a positive eigenvector* (*ϕ*, *ψ*)*.*

**Lemma 2.** *If λ*<sup>1</sup> > 0*, then the system expressed in Equation (1) is persistent and has at least one positive equilibrium. If λ*<sup>1</sup> ≤ 0*, then* (0, 0) *is globally asymptotically stable in Equation (1).*

**Lemma 3.** *If λ*<sup>1</sup> > 0 *and c and g are sufficiently small, then the system expressed in Equation (1) has a unique globally attracting positive equilibrium.*

**Remark 1.** *In the case that c* = *g* = 0*, the system expressed in Equation (1) is cooperative and hence generates a monotone semi-flow on appropriate spaces.*

#### **3. The Case of** *d***1,** *d***<sup>2</sup> Small**

Following the approach in [23], we will establish the asymptotic behavior of the principal eigenvalue of Equation (2) when *d*1, *d*<sup>2</sup> are small and, in the fully cooperative case (where *c* ≡ *g* ≡ 0), the profile of the nonnegative solutions of the corresponding steady state system for Equation (1).

$$\begin{cases} \, \, d\_1 \Delta u + r(\mathbf{x}) \boldsymbol{\nu} - s(\mathbf{x}) \boldsymbol{u} - a(\mathbf{x}) \boldsymbol{u} - b(\mathbf{x}) \boldsymbol{u}^2 = \boldsymbol{0} \text{ in } \Omega, \\\, \, \, d\_2 \Delta \boldsymbol{\nu} + s(\mathbf{x}) \boldsymbol{u} - e(\mathbf{x}) \boldsymbol{\nu} - f(\mathbf{x}) \boldsymbol{\nu}^2 = \boldsymbol{0} \text{ in } \Omega, \\\, \, \nabla \boldsymbol{u} \cdot \boldsymbol{\nu} = \nabla \boldsymbol{\nu} \cdot \boldsymbol{\nu} = \boldsymbol{0} \text{ on } \partial \Omega, \end{cases} \tag{3}$$

as well. Related results are derived in [16]. We observe that the associated kinetic system, which corresponds to Equation (1), is given by

$$\begin{cases} u\_t = r(\mathbf{x})V(\mathbf{x}) - s(\mathbf{x})\mathcal{U}(\mathbf{x}) - a(\mathbf{x})\mathcal{U}(\mathbf{x}) - b(\mathbf{x})\mathcal{U}^2(\mathbf{x}) - c(\mathbf{x})\mathcal{U}(\mathbf{x})V(\mathbf{x}) = 0, \\\ v\_t = s(\mathbf{x})\mathcal{U}(\mathbf{x}) - e(\mathbf{x})V(\mathbf{x}) - f(\mathbf{x})V^2(\mathbf{x}) - g(\mathbf{x})\mathcal{U}(\mathbf{x})V(\mathbf{x}) = 0, \end{cases} \tag{4}$$

for each *x* ∈ Ω.

For each *x*, this system shares the same properties as Equation (1) given in Lemmas 1–3, which we state for convenience.

**Lemma 4.** *Set x* ∈ Ω*. The linearization around (0,0) of Equation* (4) *has a principal eigenvalue λ*1(*x*)*. Moreover,*


Observe that, when *d*<sup>1</sup> = *d*<sup>2</sup> = 0, the eigenvalues of the linearization around (0,0) of Equation (4) are the roots of det(*A*(*x*) − *λI*), with

$$A(\mathbf{x}) = \begin{bmatrix} -(s(\mathbf{x}) + a(\mathbf{x})) & r(\mathbf{x}) \\ s(\mathbf{x}) & -\varepsilon(\mathbf{x}) \end{bmatrix} \tag{5}$$

By a simple computation, we obtain that the maximum eigenvalue is given by

$$\Lambda(\mathbf{x}) = \frac{1}{2} \left[ -(\mathbf{s}(\mathbf{x}) + a(\mathbf{x}) + \varepsilon(\mathbf{x})) + \sqrt{(\mathbf{s}(\mathbf{x}) + a(\mathbf{x}) - \varepsilon(\mathbf{x}))^2 + 4r(\mathbf{x})\mathbf{s}(\mathbf{x})} \right],\tag{6}$$

which is positive provided that (*s*(*x*) + *a*(*x*))*e*(*x*) − *r*(*x*)*s*(*x*) < 0. Our first result, which is a direct application of Theorem 1.4 of [23] (Theorem A1 in the Appendix A), states that this is indeed the necessary and sufficient condition to have a positive principal eigenvalue when *d*<sup>1</sup> and *d*<sup>2</sup> are small.

**Proposition 1.** *The principal eigenvalue λ*<sup>1</sup> *of Equation* (2) *satisfies*

$$
\lambda\_1 \to \max\_{\mathbf{x} \in \overline{\Pi}} \Lambda(\mathbf{x}) \text{ as } d\_1, d\_2 \to 0. \tag{7}
$$

*Thus, there exists a δ* > 0 *such that if*

$$\min\_{\mathbf{x}\in\Pi} ((s(\mathbf{x}) + a(\mathbf{x}))e(\mathbf{x}) - r(\mathbf{x})s(\mathbf{x}))<0,\tag{8}$$

*the principal eigenvalue of Equation* (2) *is positive for all* 0 < *d*1, *d*<sup>2</sup> < *δ, while if*

$$\min\_{\mathbf{x}\in\widetilde{\mathbf{T}}} ((s(\mathbf{x}) + a(\mathbf{x}))e(\mathbf{x}) - r(\mathbf{x})s(\mathbf{x})) > 0,\tag{9}$$

*the principal eigenvalue is negative.*

As a consequence of this result, if Equation (9) holds, the unique nonnegative equilibrium of the system expressed in Equation (3) is (0, 0), and that equilibrium is globally attracting, whenever *d*1, *d*<sup>2</sup> are small; if Equation (8) holds, the system expressed in Equation (3) is persistent, and has a positive equilibrium for small *d*1, *d*2. If Equation (8) holds and *c* ≡ *g* ≡ 0, then by Lemma 3 the system expressed in Equation (3) has a unique globally attracting positive equilibrium, which we denote by (*ud*, *vd*) with *d* = (*d*1, *d*2).

Throughout the remainder of this section, we will assume that *c* ≡ *g* ≡ 0 in Ω, in which case, the system expressed in Equation (1) is cooperative.

The next result establishes the convergence of (*ud*, *vd*) to the unique nonnegative steady state (*U*(*x*), *V*(*x*)) of the kinetic system, which satisfies

$$\begin{cases} (\mathcal{U}(\mathbf{x}), V(\mathbf{x})) \text{ is positive where } (s(\mathbf{x}) + a(\mathbf{x}))e(\mathbf{x}) - r(\mathbf{x})s(\mathbf{x}) < 0, \\ (\mathcal{U}(\mathbf{x}), V(\mathbf{x})) = 0 \text{ where } (s(\mathbf{x}) + a(\mathbf{x}))e(\mathbf{x}) - r(\mathbf{x})s(\mathbf{x}) \ge 0. \end{cases} \tag{10}$$

**Theorem 1.** *Suppose that Equation* (8) *holds. Then* (*ud*, *vd*) → (*U*, *V*) *as d* → 0 *locally uniformly in* Ω*.*

To prove this theorem, we follow the proof of Theorem 1.5 of [23], specifically their Proposition 5.2 in [23] and its hypotheses, which are listed in [23] as (A1)–(A4) and are given in the Appendix A as (L1)-(L4) to avoid confusion with the equation labels there. We should point out that Assumptions (A2) and (A3) of [23] do not hold in our case, so we cannot apply that result directly. The difference is that we allow situations where the kinetic system expressed in Equation (4) has a positive equilibrium for some values of *<sup>x</sup>* <sup>∈</sup> <sup>Ω</sup>¯ but not for others, whereas Condition (A2) requires a positive equilibrium for the kinetic system for all *x*. For that reason, we need to construct a version of the arguments in [23] that is local in *x*. Condition (A3) in [23] is used only to prove the existence of a nontrivial subsolution for a system corresponding to Equation (3), which is independent of *d*1, *d*2. We show the existence of the analogous local subsolutions we need in our case in the next lemma.

**Lemma 5.** *Suppose that x*˜ ∈ Ω0*, where*

$$\Omega\_0 = \{ \mathbf{x} \in \Omega \; / \; (s(\mathbf{x}) + a(\mathbf{x})) \mathbf{e}(\mathbf{x}) - r(\mathbf{x})s(\mathbf{x}) < 0 \}. \tag{11}$$

*Then there exists <sup>d</sup>*<sup>0</sup> > <sup>0</sup>*, <sup>ρ</sup>*<sup>0</sup> > <sup>0</sup>*, and a function <sup>w</sup>*<sup>0</sup> > <sup>0</sup> *in <sup>B</sup>*(*x*˜, *<sup>ρ</sup>*) ⊂ <sup>Ω</sup>0*, which is a subsolution of Equation* (3) *for all* 0 < *d*1, *d*<sup>2</sup> < *d*0*.*

**Proof.** Let *p* = (*p*1, *p*2), a positive eigenvector of *A*(*x*˜) with *p*<sup>1</sup> + *p*<sup>2</sup> = 1, associated with its principal eigenvalue *σ*˜ > 0. Set *ε* > 0 and small. We can choose *ρ* > 0 such that *B*(*x*˜, *ρ*) ⊂ Ω<sup>0</sup> and

$$\begin{aligned} |a(\mathbf{x}) - \vec{a}| &< \varepsilon, \ |r(\mathbf{x}) - \vec{r}| &< \varepsilon, \\ |s(\mathbf{x}) - \vec{s}| &< \varepsilon \text{ and } |\varepsilon(\mathbf{x}) - \vec{e}| &< \varepsilon \text{ for all } \mathbf{x} \text{ in } \overline{B(\vec{x}, \rho)}, \end{aligned} \tag{12}$$

where *a*˜ = *a*(*x*˜), *r*˜ = *r*(*x*˜), *s*˜ = *s*(*x*˜), and *e*˜ = *e*(*x*˜). Set *η* > 0 as the principal eigenfunction associated with *λ* > 0, the principal eigenvalue of

$$
\Delta \eta + \lambda \eta = 0 \text{ in } B(\mathfrak{x}, \mathfrak{\rho}), \ \eta = 0 \text{ on } \partial B(\mathfrak{x}, \mathfrak{\rho}), \tag{13}
$$

with max*B*(*x*˜,*ρ*) *η* = 1. We claim that we can choose *δ*,*ε*, *ρ*, *d*<sup>0</sup> > 0 such that *δηp* is a subsolution of Equation (3) for all *d*1, *d*<sup>2</sup> < *d*0. For simplicity and to keep the notation consistent with that in [23], we define *F*(*x*, *u*, *v*)=(*F*1(*x*, *u*, *v*), *F*2(*x*, *u*, *v*)), with

$$F\_1(\mathbf{x}, \boldsymbol{u}, \boldsymbol{v}) = r\boldsymbol{v} - \boldsymbol{su} - a\boldsymbol{u} - b\boldsymbol{u}^2 \text{ and } F\_2(\mathbf{x}, \boldsymbol{u}, \boldsymbol{v}) = \boldsymbol{su} - \boldsymbol{ev} - f\boldsymbol{v}^2 = 0,\tag{14}$$

where we have omitted the variable *x* in *a*, *b*, *e*, *r*, *s*, *f* to shorten the expressions. Observe that

$$\begin{aligned} F\_1(\mathbf{x}, \delta \eta \, p) &= \delta \eta (\tilde{\sigma} p\_1 + (r - \tilde{r}) p\_2 - (s - \tilde{s}) p\_1 - (a - \tilde{a}) p\_1 - b \delta p\_1^2 \eta), \\ F\_2(\mathbf{x}, \delta \eta \, p) &= \delta \eta (\tilde{\sigma} p\_2 + (s - \tilde{s}) p\_1 - (e - \tilde{e}) p\_2 - f \delta p\_2^2 \eta), \end{aligned}$$

and using Equation (12), we obtain that, if we choose *ε* > 0 and a small *δ*, we have that

$$\begin{aligned} F\_1(\mathbf{x}, \delta \eta p) &\geq \delta \eta (\tilde{\sigma} p\_1 - \varepsilon p\_2 - 2\varepsilon p\_1 - b\delta p\_1^2 \eta) > \delta \eta \frac{\delta}{2} p\_1 \\ F\_2(\mathbf{x}, \delta \eta p) &\geq \delta \eta (\tilde{\sigma} p\_2 - \varepsilon p\_1 - \varepsilon p\_2 - f \delta p\_2^2) > \delta \eta \frac{\delta}{2} p\_2. \end{aligned} \tag{15}$$

Therefore, replacing these inequalities in Equation (3), we obtain

$$d\_1 \delta p\_1 \Delta \eta + F\_1(\mathbf{x}, \delta \eta p) \ge \delta \eta \left( -d\_1 p\_1 \lambda + \frac{\tilde{\sigma}}{2} p\_1 \right)$$

$$d\_2 \delta p\_2 \Delta \eta + F\_2(\mathbf{x}, \delta \eta p) \ge \delta \eta \left( -d\_2 p\_2 \lambda + \frac{\tilde{\sigma}}{2} p\_2 \right);$$

hence, if we set *d*<sup>0</sup> = *<sup>σ</sup>*˜ <sup>2</sup>*<sup>λ</sup>* . we obtain the desired result.

Using this lemma we can follow the proof of Proposition 5.2 in [23]. To facilitate our exposition, we will use the same notation. Set the operators *D* = diag (*d*1, *d*2), L = diag (Δ, Δ). To prove Theorem 1, we will state the needed lemmas, discussing their relationships with the lemmas in [23] leading to the proof of Proposition 5.2.

Suppose that Equation (8) holds, setting *w*<sup>0</sup> = *ηδp* as in Lemma 5, and *w*<sup>0</sup> = *M* where *M* > 0 is given in Assumption (A4) such that *F*1(*x*, *u*, *v*) ≤ −*cu* and *F*2(*x*, *u*, *v*) ≤ −*cv* for all *u*, *v* ≥ *M* and *x* ∈ Ω, with *c* > 0 fixed. Set *K* > 0 such that *K* + *∂uF*1(*x*, *u*, *v*) > 0 and *K* + *∂vF*2(*x*, *u*, *v*) > 0 for all <sup>0</sup> <sup>≤</sup> *<sup>u</sup>*, *<sup>v</sup>* <sup>≤</sup> *<sup>M</sup>*, and we define *<sup>z</sup>* <sup>=</sup> *<sup>w</sup><sup>k</sup>* as the unique solution of

$$\begin{cases} -D\mathcal{L}z + \mathcal{K}z = \mathcal{K}\mu + F(x,\mu) \text{ in } \Omega\_{\prime} \\\ \nabla z \cdot \nu = 0 \text{ on } \partial\Omega\_{\prime} \end{cases}$$

for *u* = *wk*−1.

**Lemma 6.** *Suppose that Equation* (8) *holds. For every k, we have <sup>w</sup>*<sup>0</sup> <sup>&</sup>lt; *<sup>w</sup>k*+<sup>1</sup> <sup>&</sup>lt; *<sup>w</sup><sup>k</sup>* <sup>&</sup>lt; *<sup>w</sup>*0*, and as <sup>k</sup>* <sup>→</sup> <sup>∞</sup>*, <sup>w</sup><sup>k</sup> converges uniformly to the unique positive solution w of Equation* (3)*, which satisfies w*<sup>0</sup> < *w* < *w<sup>k</sup> in* Ω *for all k* ≥ 0*.*

**Proof.** We will prove that *w*<sup>0</sup> < *w<sup>k</sup>* by induction. Suppose this is true for *k*. Observe that *w*<sup>0</sup> < *w*<sup>0</sup> by construction. In the set, *<sup>B</sup>*(*x*˜, *<sup>ρ</sup>*) <sup>⊂</sup> <sup>Ω</sup><sup>0</sup> as in Lemma <sup>5</sup> *<sup>w</sup>k*+<sup>1</sup> satisfies

$$-D\mathcal{L}(\overline{w}^{k+1} - \underline{w}^0) + K(\overline{w}^{k+1} - \underline{w}^0) = K(\overline{w}^k - \underline{w}^0) + F(\mathbf{x}, \overline{w}^k) - F(\underline{w}^0),$$

in *B*(*x*˜, *ρ*). By the induction hypothesis *w*<sup>0</sup> < *wk*, whence *Kw<sup>k</sup>* + *F*(*x*, *wk*) > 0; hence, by the strong maximum principle applied to each component, we have that *wk*+<sup>1</sup> > 0 in Ω¯ . Thus, we have

$$\begin{cases} -D\mathcal{L}(\overline{w}^{k+1} - \underline{w}^0) + K(\overline{w}^{k+1} - \underline{w}^0) > 0 \text{ in } B(\tilde{x}, \rho), \\\ \overline{w}^{k+1} - \underline{w}^0 > 0 \text{ in } \partial B(\tilde{x}, \rho), \end{cases}$$

such that we have that *<sup>w</sup>k*+<sup>1</sup> <sup>−</sup> *<sup>w</sup>*<sup>0</sup> <sup>&</sup>gt; 0 in *<sup>B</sup>*(*x*˜, *<sup>ρ</sup>*). The remainder of the proof is a standard monotone iteration argument, just as in the proof of Lemma 5.3 of [23]. We observe that

$$\begin{cases} -D\mathcal{L}(\overline{w}^1 - \overline{w}^0) + K(\overline{w}^1 - \overline{w}^0) = K\overline{w}^0 + F(\mathbf{x}, \overline{w}^0) - K\overline{w}^0 < 0 \text{ in } \Omega, \\\ \nabla[\overline{w}^1 - \overline{w}^0] \cdot \nu = 0 \text{ on } \partial\Omega, \end{cases}$$

Thus, by the strong maximum principle, we have *w*<sup>1</sup> < *w*0.

Similarly, if *w<sup>k</sup>* < *wk*−1, then

$$\begin{cases} -D\mathcal{L}(\overline{\boldsymbol{w}}^{k+1} - \overline{\boldsymbol{w}}^{k}) + K(\overline{\boldsymbol{w}}^{k+1} - \overline{\boldsymbol{w}}^{k}) = \\\ K\overline{\boldsymbol{w}}^{k} + F(\boldsymbol{x}, \overline{\boldsymbol{w}}^{k}) - K\overline{\boldsymbol{w}}^{k-1} - F(\boldsymbol{x}, \overline{\boldsymbol{w}}^{k-1}) < 0 \text{ in } \Omega, \\\ \nabla[\overline{\boldsymbol{w}}^{k+1} - \overline{\boldsymbol{w}}^{k}] \cdot \nu = 0 \text{ on } \partial\Omega. \end{cases}$$

By induction, the sequence {*wk*} is decreasing, and it is bounded below by *max*{0, *<sup>w</sup>*0(*x*)}, so by standard elliptic theory, it converges to a nonnegative nontrivial solution of Equation (3) as *k* → ∞. Since by Lemma 3 the nontrivial nonnegative solution of Equation (3) is unique, it coincides with the one constructed as the limit of the sequence {*wk*}.

Define *<sup>W</sup>*<sup>0</sup> <sup>=</sup> *<sup>w</sup>*<sup>0</sup> and *<sup>W</sup>k*+<sup>1</sup> <sup>=</sup> *<sup>W</sup><sup>k</sup>* <sup>+</sup> *<sup>F</sup>*(*x*, *<sup>W</sup><sup>k</sup>* ) in Ω. Following the proof of Lemmas 5.6 and 5.7 in [23], we can prove the following result.

**Lemma 7.** *Suppose that Equation* (8) *holds. For every k, we have*

$$
\underline{w}^0 < \overline{W}^{k+1} < \overline{W}^k.
$$

*Then, <sup>W</sup><sup>k</sup> converges locally uniformly to W*<sup>∞</sup> *as k* <sup>→</sup> <sup>∞</sup>*, with*

$$\mathcal{W}^{\infty} = (\mathcal{U}(\mathfrak{x}), V(\mathfrak{x})) \text{ in } \Omega\_{0\prime} \text{ and } \mathcal{W}^{\infty}(\mathfrak{x}) = 0 \text{ in } \Omega \text{ } \backslash \Omega\_{0\prime}.$$

*where* Ω<sup>0</sup> *is given by Equation* (11)*.*

**Proof.** Observe that, by Equation (15), the function *w*<sup>0</sup> is a subsolution of the kinetic system. Repeating the proof of Lemma 5.6 in [23], or following the arguments leading to the monotonicity of the proof of Lemma <sup>6</sup> above, we have that the sequence {*W<sup>k</sup>* } is monotone decreasing and bounded below by *<sup>w</sup>*0. Therefore, *<sup>W</sup><sup>k</sup>* <sup>→</sup> *<sup>W</sup>*<sup>∞</sup> pointwise, which satisfies *<sup>F</sup>*(*x*, *<sup>W</sup>*∞) = 0, i.e. a nonnegative equilibrium of the kinetic system. Therefore, if at some *x* ∈ Ω we have that (*s*(*x*) + *a*(*x*))*e*(*x*) − *r*(*x*)*s*(*x*) ≥ 0, then *<sup>W</sup>*∞(*x*) = 0. On the other hand, *<sup>W</sup>*∞(*x*) <sup>≥</sup> *<sup>w</sup>*0(*x*) <sup>&</sup>gt; 0 in *<sup>B</sup>*(*x*˜, *<sup>ρ</sup>*) <sup>⊂</sup> <sup>Ω</sup><sup>0</sup> for *<sup>x</sup>* <sup>∈</sup> <sup>Ω</sup>0. Since *<sup>x</sup>* is arbitrary and the sequence does not depend on *w*0, we obtain that *W*∞(*x*)=(*U*(*x*), *V*(*x*)), the unique positive kinetic equilibrium, whenever *<sup>x</sup>* <sup>∈</sup> <sup>Ω</sup>0. Particularly, *<sup>W</sup>*<sup>∞</sup> is continuous. Using Theorem 5.8 of [23], we obtain that the convergence is uniform in any compact set of Ω.

**Lemma 8.** *For each k, as d*1, *<sup>d</sup>*<sup>2</sup> <sup>→</sup> <sup>0</sup>*, we have that <sup>w</sup><sup>k</sup> converges to <sup>W</sup><sup>k</sup> uniformly in* <sup>Ω</sup>¯ *.*

The proof of this result is the same as the one of Lemma 5.5 in [23].

**Proof of Theorem 7.** Using a diagonal argument, and Lemma (8), the unique positive solution *w* of Equation (3) converges to *W*<sup>∞</sup> as *d*1, *d*<sup>2</sup> → ∞.

#### **4. The Case of** *d***<sup>1</sup> and** *d***<sup>2</sup> Large**

We will start by giving a proof of a result that is well known as a "folk theorem." It is stated in slightly more generality than is needed for the specific application. For *i* = 1, ... , *N*, let *L<sup>i</sup>* denote the operator

$$L^i u = \nabla \cdot \mu\_i(\mathbf{x}) [\nabla u - \mu \nabla u\_i(\mathbf{x})] \quad \text{for} \quad \mathbf{x} \in \Omega \tag{16}$$

with no-flux boundary conditions

$$\left[\nabla\mu - \mu\nabla a\_i\right] \cdot \nu = 0 \quad \text{for} \quad \mathfrak{x} \in \partial\Omega. \tag{17}$$

Assume that *<sup>μ</sup>i*(*x*) <sup>≥</sup> *<sup>μ</sup>*<sup>0</sup> <sup>&</sup>gt; 0 on <sup>Ω</sup>¯ for all *<sup>i</sup>*. Let *<sup>A</sup>* = (*aij*(*x*)) be an *<sup>N</sup>* <sup>×</sup> *<sup>N</sup>* irreducible matrix with *aij* ≥ 0 if *i* = *j*. Consider the eigenvalue problem

$$d\_i L^i \varrho\_i + \sum\_{j=1}^N a\_{ij} \varrho\_j = \lambda \varrho\_j, \quad i = 1 \ldots N \tag{18}$$

where *di* > 0 for all *i* and *ϕ<sup>i</sup>* satisfies the boundary condition expressed in Equation (17) for each *i*. Note that, if we let Φ*<sup>i</sup>* = *exp*(−*αi*(*x*))*ϕi*, then Φ*<sup>i</sup>* satisfies Neumann boundary conditions such that the system expressed in Equation (18) rewritten in terms of the variables Φ*<sup>i</sup>* is still cooperative. Because of the classical boundary conditions, the usual results on elliptic regularity and on maximum principles for cooperative systems from [15,23] can be applied to the system for the Φ*<sup>i</sup>* values, such that the system and hence Equation (18) will have a principal eigenvalue under suitable conditions on the domain Ω and the coefficients. This idea has been used in models for single populations without an age structure or competing pairs of such populations—see, for example, [2,24].

Furthermore, we have *L<sup>i</sup>* (*exp*(*αi*(*x*)) = 0 such that the principal eigenvalue of *L<sup>i</sup>* is zero, and the eigenfunction is a multiple of *exp*(*αi*). Let *A* be the matrix defined by

$$\overline{A}\_{ij} := \frac{\int\_{\Omega} a\_{ij} \exp(\alpha\_i) dx}{\int\_{\Omega} \exp(\alpha\_i) dx}. \tag{19}$$

Denote the principal eigenvalue of Equation (18) as *<sup>λ</sup>*1(*d*) where *<sup>d</sup>* = (*d*1, ... , *dN*). Denote the principal eigenvalue of *A* as Λ.

**Lemma 9.** *Suppose that, for some <sup>γ</sup>* <sup>∈</sup> (0, 1)*, the coefficients of Equation (18) satisfy <sup>α</sup>* <sup>∈</sup> *<sup>C</sup>*2,*γ*(Ω)*, <sup>μ</sup>* <sup>∈</sup> *<sup>C</sup>*1,*γ*(Ω)*, and aij* <sup>∈</sup> *<sup>C</sup>γ*(Ω) *for <sup>i</sup>*, *<sup>j</sup>* <sup>=</sup> <sup>1</sup> ... *N, and that <sup>∂</sup>*<sup>Ω</sup> *is of class <sup>C</sup>*2,*γ. Suppose further that <sup>A</sup> is irreducible. If min*{*di* : *<sup>i</sup>* <sup>=</sup> 1, . . . *<sup>N</sup>*} → <sup>∞</sup>*, then <sup>λ</sup>*1(*d*) <sup>→</sup> <sup>Λ</sup>*.*

**Proof.** Choose any sequence *<sup>d</sup><sup>n</sup>* = (*d*1*n*, ... *dNn*) such that *min*{*din* : *<sup>i</sup>* <sup>=</sup> 1, ... *<sup>N</sup>*} → <sup>∞</sup>. Choose any subsequence, then renumber it as *dn*. Let *λ<sup>n</sup>* be the principal eigenvalue of Equation (18) corresponding to *<sup>d</sup><sup>n</sup>* and let *<sup>ϕ</sup>in*(*x*) > 0 be the *<sup>i</sup>*th component of the eigenvector, where the eigenvector is normalized by *max*{*ϕin*(*x*) : *x* ∈ Ω, *i* = 1, ... , *N*} = 1. Integrating the *i*th equation of Equation (18) over Ω and summing over *i* yields

$$\lambda\_{\Pi} \int\_{\Omega} \sum\_{i=1}^{N} \varrho\_{in}(\mathbf{x}) d\mathbf{x} = \int\_{\Omega} \sum\_{i,j=1}^{N} a\_{ij}(\mathbf{x}) \varrho\_{jn}(\mathbf{x}) d\mathbf{x} \le A\_1 \int\_{\Omega} \sum\_{i=1}^{N} \varrho\_{in}(\mathbf{x}) d\mathbf{x} \le$$

where *A*<sup>1</sup> is a constant, depending only on *A*. It follows that *λ<sup>n</sup>* is uniformly bounded from above. Similarly, *λ<sup>n</sup>* is uniformly bounded from below. Thus, any subsequence of *λ<sup>n</sup>* itself has a convergent subsequence. It then follows from dividing the *i*th equation of Equation (18) by *din* that *L<sup>i</sup> ϕin* is uniformly bounded, and *Liϕin* → 0 as *n* → ∞. By elliptic regularity, the sequence *ϕin* is uniformly bounded in *W*2,*p*(Ω) for any *p* < ∞, then by Sobolev embedding, it has a subsequence that is convergent in *C*1(Ω) and weakly convergent in *W*2,*p*(Ω). This will be true for any *i*. Taking a further subsequence if necessary and renumbering again, we obtain a sequence where *λ<sup>n</sup>* → *λ*<sup>∗</sup> for some *λ*<sup>∗</sup> and *ϕin* → *ϕ*<sup>∗</sup> *<sup>i</sup>* for all *i*, with *Liϕ*<sup>∗</sup> *<sup>i</sup>* = 0. We then must have *ϕ*<sup>∗</sup> *<sup>i</sup>* = *ciexp*(*αi*) for some nonnegative constant *ci*, and with *max*{*ϕ*<sup>∗</sup> *<sup>i</sup>* (*x*) : *x* ∈ Ω, *i* = 1 ... *N*} = 1. Integrating Equation (18) over Ω and using the no-flux boundary conditions gives

$$\sum\_{j=1}^{N} \left[ \int\_{\Omega} a\_{\vec{i}\vec{j}}(\mathbf{x}) \, \boldsymbol{\varrho}\_{\vec{j}}^{\*}(\mathbf{x}) d\mathbf{x} \right] = \lambda^{\*} \int\_{\Omega} \boldsymbol{\varrho}\_{\vec{i}}^{\*} , \quad \vec{i} = 1 \dots N , \tag{20}$$

such that *<sup>N</sup>*

$$\sum\_{j=1}^{N} \left[ \frac{\int\_{\Omega} a\_{ij}(\mathbf{x}) \exp(a\_j(\mathbf{x})) d\mathbf{x}}{\int\_{\Omega} \exp(a\_i(\mathbf{x})) d\mathbf{x}} \right] \mathbf{c}\_j = \lambda^\* \mathbf{c}\_i, \quad i = 1 \ldots N. \tag{21}$$

It follows that (*c*1, ... , *cN*) must be a nontrivial nonnegative eigenvector of *A* with the normalization prescribed by *max*{*ciexp*(*αi*(*x*)) : *x* ∈ Ω, *i* = 1, ... *N*} = 1. These last conditions uniquely determine the limits of the subsequence of the original subsequence {*λ*(*dn*), *ϕn*}. Since every subsequence of the original sequence {*λ*(*dn*), *ϕn*} has a subsequence converging to the values determined by Equation (21), the same must be true for the original sequence. Since the original sequence of values {*dn*} could be any increasing sequence that approaches infinity as *n* → ∞, the conclusion of the lemma follows.

In the specific system expressed in Equation (1) that we consider, *L<sup>i</sup>* = Δ, such that *α<sup>i</sup>* and *μ<sup>i</sup>* are constants. In that case, we have *Aij* = *aij*, where *aij* is the average of *aij* over Ω. Denote the averages of the coefficients in Equation (1) by *r*¯, *s*¯, etc. Calculations analogous to those in Equation (5), Equation (6), and the related discussion then yield the following:

**Corollary 1.** *Suppose that the hypotheses of Lemma 9 are satisfied. There exists a D* > 0 *such that if*

$$\vec{\varepsilon}(\vec{s} + \vec{a}) - \vec{r}\vec{s} < 0\_{\prime\prime}$$

*the principal eigenvalue λ*<sup>1</sup> *of Equation (2) is positive for all d*1, *d*<sup>2</sup> > *D, while if*

$$
\vec{e}(\vec{s} + \vec{a}) - \vec{r}\vec{s} > 0,
$$

*the principal eigenvalue is nonpositive.*

**Remark 2.** *In the ODE system corresponding to Equation (1) with coefficients averaged over* Ω*, one can compute R*<sup>0</sup> *as r*¯*s*¯/[*e*¯(*s*¯ + *a*¯)] *via the methods of [25]. The first inequality in Corollary 1 is equivalent to R*<sup>0</sup> > 1*, while the second is equivalent to R*<sup>0</sup> < 1*. By writing R*<sup>0</sup> = [*r*¯/(*s*¯ + *a*¯)][*s*¯/*e*¯]*, we can interpret the condition for persistence as saying that the products of the ratios of the growth terms over the loss terms for adults and juveniles should be greater than* 1 *for persistence.*

#### **5. General Diffusion Rates**

**Case1**: Persistence or extinction for all diffusion rates

*Mathematics* **2020**, *8*, 396

#### **Proposition 2.** *If*

$$\int\_{\Omega} \sqrt{rs} \, dx - \frac{1}{2} \int\_{\Omega} (s + a + e) dx > 0 \tag{22}$$

*then λ*<sup>1</sup> > 0 *for all positive diffusion rates.*

*If*

$$\min\_{x \in \overline{\Omega}} \left[ 4(s(x) + a(x))e(x) - (r(x) + s(x))^2 \right] > 0 \tag{23}$$

*then λ*<sup>1</sup> < 0 *for all positive diffusion rates.*

**Proof.** If we divide the first equation in Equation (2) by *ϕ* and integrate over Ω, using Green's formula to integrate the term Δ*ϕ*/*ϕ*, we obtain the inequality

$$|\Omega|\lambda\_1 \ge \int\_{\Omega} r \left(\frac{\Psi}{\Psi}\right) d\mathbf{x} - \int\_{\Omega} (\mathbf{s} + a) d\mathbf{x}.\tag{24}$$

Similarly, if we divide the second equation by *ψ* and integrate we obtain

$$|\Omega|\lambda\_1 \ge \int\_{\Omega} \mathbf{s}\left(\frac{\boldsymbol{\varrho}}{\boldsymbol{\Psi}}\right) d\mathbf{x} - \int\_{\Omega} \boldsymbol{\varrho} \, d\mathbf{x}.\tag{25}$$

If we add Equations (24) and (25) and divide by 2, we obtain

$$
\lambda\_1 \ge \frac{1}{2|\Omega|} \left( \int\_{\Omega} \left[ r \left( \frac{\Psi}{\varrho} \right) + s \left( \frac{\varrho}{\Psi} \right) \right] d\mathbf{x} - \int\_{\Omega} (s + a + \mathfrak{e}) d\mathbf{x} \right). \tag{26}
$$

By Cauchy's inequality, *rz* + *sz*−<sup>1</sup> ≥ <sup>2</sup> <sup>√</sup>*rs* for all *<sup>z</sup>* <sup>&</sup>gt; 0, so from Equation (26) we obtain

$$\lambda\_1 \ge \frac{1}{|\Omega|} \left[ \int\_{\Omega} \sqrt{r\mathbf{s}} \, d\mathbf{x} - \frac{1}{2} \int\_{\Omega} (\mathbf{s} + \mathbf{a} + \mathbf{e}) d\mathbf{x} \right] \tag{27}$$

Thus, *λ*<sup>1</sup> > 0 if Equation (22) holds, so the first part of Proposition 2 holds. Going in the other direction, if we multiply the first equation of Equation (2) by *ϕ* and integrate, using integration by parts on the *ϕ*Δ*ϕ* term, and similarly multiply the second equation by *ψ* and integrate, and then add the results, we get

$$
\lambda\_1 \int\_{\Omega} (\varrho^2 + \Psi^2) d\mathbf{x} \le \int\_{\Omega} [-(s+a)\varrho^2 + (r+s)\varrho\psi - \varepsilon\psi^2] d\mathbf{x} \,\tag{28}
$$

The integrand on the right side of Equation (28) is a quadratic form in *ϕ* and *ψ*, which will be negative definite if

$$4(s+a)e \succ (r+s)^2,\tag{29}$$

so *λ*<sup>1</sup> < 0 if Equation (23) holds, which proves the second part of Proposition 2.

Remarks: Note that the first integral in Equation (22) is what appears in the formula for the Bhattacharyya coefficient [26,27], which is used to compare how well probability distributions match each other. Specifically, if two probability distributions *P* and *Q* have probability density functions *<sup>p</sup>*(*x*) and *<sup>q</sup>*(*x*) for *<sup>x</sup>* <sup>∈</sup> *<sup>U</sup>* <sup>⊂</sup> <sup>R</sup>*n*, the Bhattacharyya coefficient is

$$BC(P,Q) = \int\_{\mathcal{U}} \sqrt{p(x)q(x)}dx.$$

For any *P* and *Q*, 0 ≤ *BC*(*P*, *Q*) ≤ 1. If *BC*(*P*, *Q*) = 1, then *P* and *Q* are the same, that is, *p* = *q* a.e. If *BC*(*P*, *Q*) = 0, then the supports of *p* and *q* are disjoint. If we write *r*(*x*) = *r*0*ρ*(*x*) and *s*(*x*) = *s*0*σ*(*x*)

such that <sup>Ω</sup> *<sup>ρ</sup>*(*x*) = <sup>Ω</sup> *σ*(*x*) = 1, then we can compute *r*<sup>0</sup> = *r*¯|Ω| and *s*<sup>0</sup> = *s*¯|Ω|. We can treat *ρ* and *σ* as if they were probability density functions for distributions *R* and *S*. We then have

$$\int\_{\Omega} \sqrt{r} \text{d}x = |\Omega| \sqrt{r} \text{SBC}(\mathbb{R}, \mathbb{S}). \tag{30}$$

The maximum of *BC*(*R*, *S*) is 1, corresponding to the case where *r* and *s* are multiples of each other, and the minimum is 0, corresponding to the case where the supports of *r* and *s* are disjoint. Thus, the degree to which *ρ* and *σ* match each other has a strong impact on the estimate for *λ* in Equation (27).

Using Equation (30) and the fact that *BC*(*R*, *S*) ≤ 1 in Equation (22) shows that Equation (22) implies 2√*r*¯*s*¯ <sup>&</sup>gt; [(*s*¯ <sup>+</sup> *<sup>a</sup>*¯) + *<sup>e</sup>*¯]. Squaring both sides and using Cauchy's inequality implies *<sup>e</sup>*¯(*s*¯ <sup>+</sup> *<sup>a</sup>*¯) <sup>−</sup> *r*¯*s*¯ < 0 as in the first case of Corollary 1. Similarly, if Equation (22) holds, then 2*r*(*x*)*s*(*x*) > (*s*(*x*) + *a*(*x*)) + *e*(*x*) for some *x* ∈ Ω, and it then follows in the same way that the inequality in the first case of Proposition 1 holds. If Equation (23) holds, then Equation (29) holds, and then by Cauchy's inequality the second case of Proposition 1 holds. Thus, the conditions expressed in Equation (22) and Equation (23) in Proposition 2, which imply *λ*<sup>1</sup> > 0 or *λ*<sup>1</sup> < 0 for all diffusion rates, also imply some of the corresponding conditions we have obtained for either large or small diffusion rates.

In the situation where the spatial distributions of habitat quality *r* for reproduction by adults and *s* for survival and maturation of juveniles into adults are perfectly correlated, such that *r*(*x*) = *r*1*s*(*x*) for some constant *r*1, the eigenvalue problem expressed in Equation (2) can be rewritten as a weighted symmetric eigenvalue problem by multiplying the second equation in Equation (2) by *r*1, which yields

$$\begin{aligned} d\_1 \Delta \varrho(\mathbf{x}) - (s(\mathbf{x}) + a(\mathbf{x})) \varrho + r(\mathbf{x}) \psi &= \lambda \varrho \\ d\_2 r\_1 \Delta \psi + r(\mathbf{x}) \varrho - r\_1 \varepsilon(\mathbf{x}) \psi &= \lambda r\_1 \psi. \end{aligned} \tag{31}$$

The principal eigenvalue for Equation (31) has a variational characterization of *λ*<sup>1</sup> as

$$\lambda\_1 = \max\_{\boldsymbol{\varrho}, \boldsymbol{\varrho} \in \mathcal{W}^{1,2}(\Omega)} \frac{\int\_{\Omega} (-d\_1 |\nabla \boldsymbol{\varrho}|^2 - d\_2 r\_1 |\nabla \boldsymbol{\psi}|^2 - (s + a) \boldsymbol{\varrho}^2 + 2r \boldsymbol{\varrho} \boldsymbol{\psi} - r\_1 c \boldsymbol{\psi}^2) d\mathbf{x}}{\int\_{\Omega} (\boldsymbol{\varrho}^2 + r\_1 \boldsymbol{\psi}^2) d\mathbf{x}}. \tag{32}$$

It follows in that case that *λ*<sup>1</sup> is decreasing in both *d*<sup>1</sup> and *d*2, such that slower diffusion is advantageous.

**Case 2**: Asymptotic behavior for large reproductive rates

Suppose that *r*(*x*) = *nr*0(*x*) and that *s*(*x*)*r*0(*x*) > 0 for *x* ∈ Ω<sup>0</sup> with Ω<sup>0</sup> = ∅; therefore, there is a region where both the adult reproduction rate and the juvenile maturation rate are positive. The factor *n* scales the reproductive rate of adults in regions where *r*0(*x*) > 0. For any fixed diffusion rates, it turns out that, for sufficiently large values of the scaling coefficient *n*, the principal eigenvalue of Equation (2) is positive, so the system expressed in Equation (1) is persistent. We will characterize the asymptotic behavior of the principal eigenvalue as *n* → ∞. If we make the further assumption that *g* ≡ *c* ≡ 0, then the system expressed in Equation (1) is cooperative, so for an *n* large enough that the principal eigenvalue of Equation (2) is positive, Equation (1) has a unique positive equilibrium, and we will characterize the behavior of that equilibrium as *n* → ∞ as well in that case.

Let *λ<sup>n</sup>* <sup>1</sup> denote the principal eigenvalue for Equation (2) with *r*(*x*) = *nr*0(*x*). Observe that, since *<sup>s</sup>*(*x*)*r*0(*x*) > 0 for *<sup>x</sup>* ∈ <sup>Ω</sup><sup>0</sup> and <sup>Ω</sup><sup>0</sup> = <sup>∅</sup>, Proposition <sup>2</sup> implies that *<sup>λ</sup><sup>n</sup>* <sup>1</sup> > 0 for *n* sufficiently large, and in fact by Equation (27), *λ<sup>n</sup>* <sup>1</sup> → ∞ as *n* → ∞. The following proposition states the asymptotic behavior of *λn* <sup>1</sup> as *n* → ∞.

**Proposition 3.** *If r*(*x*) = *nr*0(*x*) *and s*(*x*)*r*0(*x*) > 0 *for x* ∈ Ω<sup>0</sup> *with* Ω<sup>0</sup> = ∅*, then*

$$\lim\_{n \to \infty} \frac{\lambda\_1^n}{\sqrt{n}} \to \max\_{x \in \overline{\Omega}} (\sqrt{r\_0(x)s(x)}).$$

**Proof.** We start by noting that, if *λ<sup>n</sup>* <sup>1</sup> , then (*ϕn*, *ψn*) are the principal eigenvalue and corresponding eigenfunction of Equation (2) for *r*(*x*) = *nr*0(*x*); therefore, *λ<sup>n</sup>* 1/ <sup>√</sup>*n*, *<sup>ϕ</sup>*5*<sup>n</sup>* <sup>=</sup> *<sup>ϕ</sup>n*, and *<sup>ψ</sup>*5*<sup>n</sup>* <sup>=</sup> <sup>√</sup>*nψ<sup>n</sup>* are the principal eigenvalue and eigenfunction of the following problem:

$$\begin{cases} \frac{d\_1}{\sqrt{n}} \Delta \hat{\boldsymbol{\mu}} - \frac{(\boldsymbol{s}(\boldsymbol{x}) + \boldsymbol{a}(\boldsymbol{x}))}{\sqrt{n}} \hat{\boldsymbol{\eta}} + r\_0(\boldsymbol{x}) \hat{\boldsymbol{\psi}} &= \hat{\lambda} \hat{\boldsymbol{\eta}} \text{ in } \Omega, \\\ \frac{d\_2}{\sqrt{n}} \Delta \hat{\boldsymbol{\psi}} - \frac{\boldsymbol{e}(\boldsymbol{x})}{\sqrt{n}} \hat{\boldsymbol{\psi}} + \boldsymbol{s}(\boldsymbol{x}) \hat{\boldsymbol{\varphi}} &= \hat{\lambda} \hat{\boldsymbol{\psi}} \text{ in } \Omega, \\\ \nabla \hat{\boldsymbol{\varphi}} \cdot \boldsymbol{\nu} = \nabla \hat{\boldsymbol{\psi}} \cdot \boldsymbol{\nu} &= 0 \text{ on } \partial \Omega. \end{cases} \tag{33}$$

Considering the elliptic operators *L*1*u* = *d*1Δ*u* − (*s*(*x*) − *a*(*x*))*u* and *L*2*v* = *d*2Δ*v* − *e*(*x*)*v*, *<sup>D</sup>* <sup>=</sup> diag √1 *<sup>n</sup>* , <sup>√</sup><sup>1</sup> *n* and L = diag(*L*1, *L*2), the system expressed in Equation (33) satisfies the hypothesis of Theorem 1.4 of [23]. Thus, *n* → ∞

$$\frac{\lambda\_1^n}{\sqrt{n}} \to \max\_{x \in \overline{\Omega}} \lambda(A(x)),$$

where

$$A(x) = \begin{pmatrix} 0 & r\_0(x) \\ s(x) & 0 \end{pmatrix},$$

which has eigenvalues ± *r*0(*x*)*s*(*x*), from whence the result follows.

In the case where *g* ≡ *c* ≡ 0 such that Equation (1) is cooperative, Lemma 3 implies that Equation (1) has a unique positive equilibrium if the principal eigenvalue of Equation (2) is positive. The next result states the asymptotic behavior of the unique positive equilibrium of Equation (1) for *n* large in that case.

**Proposition 4.** *Suppose the hypotheses of Proposition 3 are satisfied and that g* ≡ *c* ≡ 0*. Let* (*un*, *vn*) *be the unique positive solution of Equation* (3)*. Then n*<sup>−</sup> <sup>2</sup> <sup>3</sup> (*un*, *vn*) <sup>→</sup> (*U*∞, *<sup>V</sup>*∞) *uniformly in* <sup>Ω</sup>¯ *where*

$$\begin{split} \, ^\circ L^\infty(\mathbf{x}) &= \frac{r\_0(\mathbf{x})^{\frac{2}{3}}}{b(\mathbf{x})^{\frac{2}{3}}} \frac{\mathbf{s}(\mathbf{x})^{\frac{1}{3}}}{f(\mathbf{x})^{\frac{1}{3}}}, \quad \, ^\circ V^\infty(\mathbf{x}) = \frac{\mathbf{s}(\mathbf{x})^{\frac{2}{3}}}{f(\mathbf{x})^{\frac{2}{3}}} \frac{r\_0(\mathbf{x})^{\frac{1}{3}}}{b(\mathbf{x})^{\frac{1}{3}}} \, \, \text{when } \mathbf{s}(\mathbf{x}) r\_0(\mathbf{x}) > 0 \\\\ \, ^\circ L^\infty(\mathbf{x}) &= 0, \, V^\infty(\mathbf{x}) = 0 \, \text{when } \mathbf{s}(\mathbf{x}) r\_0(\mathbf{x}) = 0. \end{split} \tag{34}$$

**Proof.** To prove this result, we use a different scaling. After some simple computations, we find that

$$(w\_n, z\_n) = \left(\mu\_n n^{-\frac{2}{3}}, v\_n n^{-\frac{1}{3}}\right) \dots$$

is the unique positive solution of the scaled system

$$\begin{cases} \begin{aligned} &n^{-\frac{2}{3}}[d\_1\Delta w - (s(\mathbf{x}) + a(\mathbf{x}))w] + r\_0(\mathbf{x})z - b(\mathbf{x})w^2 = 0 \text{ in } \Omega, \\ &n^{-\frac{1}{3}}[d\_2\Delta z - \varepsilon(\mathbf{x})z] + s(\mathbf{x})w - f(\mathbf{x})z^2 = 0 \text{ in } \Omega, \\ &\nabla w \cdot \nu = \nabla z \cdot \nu = 0 \text{ on } \partial\Omega. \end{aligned} \end{cases} \tag{35}$$

We set the operators

$$L\_1 w = d\_1 \Delta w - (s(x) + a(x))w \text{ and } L\_2 z = d\_2 \Delta z - e(x)z,$$

$$D = \text{diag}\left(n^{-\frac{2}{3}}, n^{-\frac{1}{3}}\right) \text{ and } \mathcal{L} = (L\_1, L\_2), \text{ and}$$

$$F(x, w, z) = (r\_0(x)z - b(x)w^2, s(x)w - f(x)z^2).$$

To prove the proposition, we can follow the same steps as in the proof of Theorem 1. Indeed, the principal eigenvalue of the linearization around (0, 0) of the associated kinetic system of Equation (35) is *s*(*x*)*r*0(*x*), and when it is positive, the kinetic equilibrium is given by the right hand side of Equation (34). This concludes the proof.

#### **6. Conclusions**

The most fundamental conclusion from our analysis is that reaction-diffusion models for populations with a stage structure in spatially heterogeneous environments do not necessarily predict that slower diffusion is advantageous for persistence. This is in contrast to the case where populations are structured only by spatial location, where a version of the reduction principle [6] applies; in a competition between otherwise identical populations with different diffusion rates, the prediction is that "the slower diffuser wins" [4,5]. The mechanism underlying this observation is that, in our structured model, the regions where it is possible for adults to produce offspring may be separated from those where juveniles can survive and mature into adults. The conditions we find that imply persistence generally require that the product *r*(*x*)*s*(*x*) of the reproductive rate of adults and the maturation rate of juveniles be sufficiently large relative to their death rates. For slow diffusion, the condition for persistence is that *r*(*x*)*s*(*x*) > *e*(*x*)(*s*(*x*)+(*a*(*x*)) at some point *x* ∈ Ω. For fast diffusion, it is *r*¯*s*¯ > *e*¯(*s*¯ + *a*¯) where *r*¯,*s*¯,*e*¯, and *a*¯ are the spatial averages of those quantities. If the spatial distributions of *r* and *s* are closely correlated and are large in a few places but small in most, such that the maximum of *rs* is large but the averages *r*¯ and *s*¯ are small, the condition for persistence with slow diffusion may be satisfied, while the condition with fast diffusion may fail. In that type of environment, slow diffusion is clearly favored. Furthermore, if *r* and *s* are perfectly correlated in the sense that they are multiples of each other, the principal eigenvalue determining the growth rate of the population at low density is decreasing with respect to the diffusion rates, as in the case of unstructured populations in heterogeneous environments. On the other hand, if both *r* and *s* are large on some regions but very small outside of them, and the regions where they are large are disjoint (that is, separated from each other), then the product *rs* could be small everywhere, but the averages *r*¯ and *s*¯ could be large. In that case, the condition for persistence with small diffusion may fail, but the condition with fast diffusion may be satisfied, such that fast diffusion is favored.

We found that a sufficient condition for persistence for all diffusion rates is

$$\int\_{\Omega} \sqrt{rs} \, d\mathbf{x} - \frac{1}{2} \int\_{\Omega} (s + a + e) d\mathbf{x} > 0.$$

The first term can be written as <sup>√</sup>*r*¯*s*¯|Ω|*BC*(*r*(*x*)/*r*¯,*s*(*x*)/*s*¯), where *BC* denotes the Bhattacharrya coefficient (see [26,27]), which measures how closely probability densities match each other. For distributions that are equal to each other, *BC* = 1, but for distributions that are mutually exclusive in the sense that the regions where they are positive do not intersect, *BC* = 0. This observation again shows that the degree to which the spatial distributions of *r* and *s* match each other is significant in determining the predictions of the model expressed in Equation (1).

Finally, we found that, if we scale the adult reproductive rate as *r*(*x*) = *nr*0(*x*) and there is some overlap between the distributions of *r* and *s* such that *s*(*x*)*r*0(*x*) > 0 on some subset of Ω with positive measure, then for any fixed diffusion rates the system expressed in Equation (1) will be persistent if *n* is sufficiently large. This means a population with any diffusion rates can persist if there is even a modest overlap between the regions where adults can reproduce and where juveniles can mature, provided that the reproductive rate of adults is sufficiently large. We characterized the asymptotic behavior as *n* → ∞ of the principal eigenvalue of Equation (2). In the cooperative case where *g* ≡ *c* ≡ 0, the system expressed in Equation (1) will have a unique positive equilibrium if it is persistent, and in that case we also characterized the asymptotic behavior as *n* → ∞ of the equilibrium.

There are several directions for further research on the general topic of this paper. It would be of interest to take the approach of [4] and consider the competition between two stage structured populations described by systems such as Equation (1) that differ only in their diffusion rates. That would be somewhat challenging because it would involve systems of four equations, but at least in the cooperative case where *c* = *g* = 0 the general theory of monotone dynamical systems as in [28] and some of the ideas and methods of [10] would apply. It would also be interesting to consider models with an explicit age structure, as introduced in [11] and studied in [12,13]. Finally, it would be interesting albeit challenging to consider the case of time-periodic environments with spatial heterogeneity. Temporal variation alone is sufficient to cause faster diffusion to be favored in such environments in some cases (see [29]), but even without a stage structure, the time-dependent case is challenging, and there are many open questions.

**Author Contributions:** Conceptualization, R.S.C. C.C., and S.M.;mathematical analysis, R.S.C., C.C. and S.M.; writing–original draft preparation, C.C. and S.M.; writing–review and editing, R.S.C., C.C., and S.M.; funding acquisition, R.S.C., C.C. and S.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received external funding as noted in the subsequent Acknowledgments .

**Acknowledgments:** R.S.C. and C.C. are supported in part by NSF Awards DMS 15-14792 and 18-53478. S.M. is supported by CONICYT + PIA/Concurso de Apoyo a Centros Científicos y Tecnológicos de Excelencia con Financiamiento Basal AFB170001.

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

#### **Appendix A**

In [23], the authors considered the equilibria and dynamics of the system

$$\begin{cases} \frac{\partial u}{\partial t} = D \mathcal{L}u + F(x, u) \text{ in } \Omega \times (0, \infty), \\\\ \mathcal{B}u = 0 \text{ on } \partial \Omega \times (0, \infty) \end{cases} \tag{A1}$$

where *u* = (*u*1, ... , *un*)*<sup>T</sup>* is a vector of smooth functions, Ω is a bounded domain in R*<sup>N</sup>* with smooth boundary, *u* = (*u*1(*x*), ... *un*)*<sup>T</sup>* is a vector of smooth functions, *D* = *diag*(*d*1, ... *dn*) is a diagonal matrix of positive constants, L = *diag*(*L*1, ... *<sup>L</sup>n*) is a diagonal matrix of second order uniformly strongly elliptic operators of the form

$$L^i = \sum\_{j,k=1}^N \mathfrak{a}^i\_{jk} \frac{\partial^2}{\partial \mathfrak{x}\_j \partial \mathfrak{x}\_k} + \sum\_{j=1}^N \beta^i\_j \frac{\partial}{\partial \mathfrak{x}\_j} + \gamma^i$$

with smooth coefficients, and B = (*B*1, ... *Bn*), where for each *i*, *Bi* defines a Dirichlet, Neumann, or Robin boundary condition. (They include Neumann as a case of Robin.) They also considered the associated linearized problem, which they wrote as

$$\begin{cases} \,^\*D\mathcal{L}\phi + Au\phi = -\lambda\phi \text{ in } \Omega, \\\\ \,^\*B\phi = 0 \text{ on } \partial\Omega, \end{cases} \tag{A2}$$

where *A* = (*aij*) is an *n* × *n* matrix of smooth functions with *aij* ≥ 0 for *i* = *j*, and *φ* = (*φ*1(*x*),... *φn*(*x*))*<sup>T</sup>* is a vector of smooth functions.

Note that, in our notation, we use the opposite sign convention to the one used in [23], such that what they denote as −*λ*, we denote as *λ*.

For details on what the specific smoothness assumptions require, see [23]. Under those assumptions, by the Perron-Frobenius theorem, for each *x* ∈ Ω, the matrix *A* has a principal eigenvalue, which in our notation we denote as Λ(*x*). The first major result of [23] is their Theorem 1.4, which can be stated as follows:

**Theorem A1.** *(Theorem 1.4 of [23]) The principal eigenvalue λ*<sup>1</sup> *of the system expressed in Equation* (A2) *with Dirichlet, Neumann, or Robin boundary conditions satisfies*

$$\lim\_{\max\{d\_1,\ldots,d\_n\}\to 0} \lambda\_1 = -\max\_{x\in\Omega} \Lambda(x).$$

Except for the adjustments needed for our different notation, that theorem applies directly to our system in all cases.

The second major result of [23] gives conditions under which the system expressed in Equation (A1) with small diffusion rates has the same dynamics as the kinetic system

$$\frac{d\mathcal{U}\_i}{dt} = F\_i(\mathbf{x}, \mathcal{U}\_1, \dots, \mathcal{U}\_n) \text{ for } i = 1, \dots, n. \tag{A3}$$

The conditions can be stated as follows, noting that we have replaced the A used in [23] with L to avoid any confusion with equation numbers in this Appendix:

**(L1)** *∂Fi*/*∂Uj* ≥ 0 (i.e. the systems expressed in Equations (A1) and (A3) are cooperative).

**(L2)** For each *<sup>x</sup>*<sup>0</sup> <sup>∈</sup> <sup>Ω</sup>¯ , the system expressed in Equation (A3) has a unique positive equilibrium *α*(*x*0), which is globally asymptotically stable among positive solutions and is locally linearly stable, and *α*(*x*) depends continuously on *x*.

**(L3)** There is a *<sup>δ</sup>*<sup>0</sup> <sup>&</sup>gt; 0 such that, for *<sup>j</sup>* <sup>=</sup> 1, ... , *<sup>n</sup>*, *Fj*(*x*, *<sup>U</sup>*)/*Uj* <sup>&</sup>gt; *<sup>δ</sup>*<sup>0</sup> for all *<sup>x</sup>* <sup>∈</sup> <sup>Ω</sup>¯ provided 0 < *Ui* ≤ *δ*<sup>0</sup> for *i* = 1, . . . , *n*.

**(L4)** There is a *δ* <sup>0</sup>, *M* > 0 such that for *j* = 1, ... , *n*, *Fj*(*x*, *U*)/*Uj* < −*δ* <sup>0</sup> for all *<sup>x</sup>* <sup>∈</sup> <sup>Ω</sup>¯ provided *Ui* ≥ *M* for *i* = 1, . . . , *n*.

The second result of [23] that we use is Proposition 5.2, which can be stated as follows:

**Theorem A2.** *(Proposition 5.2 of [23]) For any positive steady state wd of Equation* (A1)*, we have that w* → *α uniformly in* Ω *as max*{*d*1,..., *dn*} → 0*, with the α given in Assumption (L2).*

#### **References**


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

### *Article* **Using** *G***-Functions to Investigate the Evolutionary Stability of Bacterial Quorum Sensing**

#### **Anne Mund, Christina Kuttler and Judith Pérez-Velázquez \***

Zentrum Mathematik, Technische Universität München Boltzmannstr. 3, 85748 Garching, Germany; mund@ma.tum.de (A.M.); kuttler@ma.tum.de (C.K.)

**\*** Correspondence: cerit@ma.tum.de

Received: 30 October 2019; Accepted: 13 November 2019; Published: 15 November 2019

**Abstract:** In ecology, *G*-functions can be employed to define a growth function *G* for a population *b*, which can then be universally applied to all individuals or groups *bi* within this population. We can further define a strategy *vi* for every group *bi*. Examples for strategies include diverse behaviour such as number of offspring, habitat choice, and time of nesting for birds. In this work, we employ *G*-functions to investigate the evolutionary stability of the bacterial cooperation process known as quorum sensing. We employ the *G*-function ansatz to model both the population dynamics and the resulting evolutionary pressure in order to find evolutionary stable states. This results in a semi-linear parabolic system of equations, where cost and benefit are taken into account separately. Depending on different biological assumptions, we analyse a variety of typical model functions. These translate into different long-term scenarios for different functional responses, ranging from single-strategy states to coexistence. As a special feature, we distinguish between the production of public goods, available for all subpopulations, and private goods, from which only the producers can benefit.

**Keywords:** Evolutionary dynamics; G-function; Quorum Sensing; Public Goods; semi-linear parabolic system of equations

#### **1. Introduction**

Since its discovery in 1977 [1], quorum sensing (QS) has received increasing attention as a mechanism of bacterial communication. Nowadays, we know that QS coordinates a range of behaviours at the population level [2]. One of those is the production of iron-scavenging molecules, so-called siderophores, in *Pseudomonas aeruginosa*. The bacteria secrete the siderophores and other virulence factors into the surrounding medium, where their benefit can be shared by all cells in the local population, leading to the term *public goods* [3,4], e.g., with the focus on kin selection in [5] or the conflict between individual and group interests in [6]. However, such a cooperative system is vulnerable to exploitation by non-cooperative cheaters. Gaining all the benefits of QS without incurring the metabolic costs, these cheater mutants have a growth advantage compared to the wild-type. Studies have shown that such mutants arise in laboratory conditions as well as in infection models [7–13] and act as social cheats. A few examples of these social cheats can be found in [3,14,15] (all considering the Gram-negative *P. aeruginosa*), while Pollitt et al. [16] provided an example from the Gram-positive *S. aereus*.

This gives rise to questions about the evolutionary stability of QS. Once cheaters arise (e.g., through loss-of-function mutation), they should theoretically outgrow the producers as they have more resources available to invest in cell division. Cheaters have been shown to outcompete producers both in vitro and in vivo, but QS seems to be evolutionary stable in natural systems nevertheless. In contrast to the

aforementioned public goods, *private goods* are only accessible to the producing cell itself. They are hence innately protected from cheaters and provide their benefit exclusively for cells with a functioning QS system. Apart from extracellular molecules, QS also controls the production of proteins which act within the cell. In *P. aeruginosa*, one such protein is nucleoside hydrolase (Nuh). As Nuh is involved in metabolising adenosine, only bacteria with intact signal receptors can digest this carbon source. In this way, cooperation via QS provides a private fitness benefit to cooperating cells if adenosine is available as carbon source [17]. Schuster et al. [18] discussed the benefits attained from private goods. Specifically, on the question why public goods are ultimately stable, they suggested that, in fact, the private good is, directly or indirectly, involved in cooperative behaviour.

Extracellular molecules do not provide benefit indiscriminately, but are limited by diffusion and habitat structure. Kümmerli et al. [19] found a negative correlation between habitat structure and water solubility of siderophores, a class of secreted enzymes under QS control in a wide range of bacteria. For highly structured environments such as animal tissues, water solubility of siderophores is high, while microstructures in the environment naturally limit the resulting diffusion. Conversely, water solubility of siderophores is low in unstructured environments such as water habitats. This leads to siderophores clinging to each other as well as to lipid membranes. In this way, a fraction of the siderophores stay with their producer and provide some private benefit.

Several mechanisms, such as kin selection [20] and policing [21], have been described that could explain the evolutionary stability of cooperation and QS despite the advantages cheaters have in such a system.

In addition to experimental approaches, there have been many different mathematical modelling approaches to study QS, from stochastic cellular automata [22] to classic individual-based models [23], including differential equations [24]. In this paper, we develop a mathematical model to investigate the evolutionary stability of QS, which results in a semi-linear parabolic system of equation, following Frank [25] and Mund et al. [26], who used an age-dependent model as basis to check evolutionary stability. We explore the population dynamics between secretors and cheaters and are also interested in the evolutionary pressure that arises through those population dynamics. As such, we employ the *G*-function ansatz that was introduced by Cohen et al. [27] and further developed since then [28–31]. With the *G*-function ansatz, one finds a function *G* which describes the growth rate of a subpopulation. This function is dependent on the *strategy* of this subpopulation and both population numbers and strategies of other existing subpopulations (see Figure 1). A *strategy* can be taken to be any behaviour of interest which can be assigned a real number.

In our case, we take the QS-intensity of the subpopulation as the strategy. A strategy value of zero therefore denotes a QS cheater, while we assume a value of one to be the normal output of QS active secretors. If we take *vi* as the strategy of the subpopulation, *v* is the vector of strategy values for all subpopulations. Accordingly, *bi* and *b* represent the population numbers of subpopulation *i* and the complete vector, respectively. The population dynamics for subpopulation *i* thus follows the equation

$$
\dot{b}\_i'(t) = G(\upsilon\_i(t), \upsilon(t), b(t)) \cdot b\_i(t). \tag{1}
$$

At the same time, the per-capita growth function *G* also governs the evolutionary dynamics of the population. If we calculate the derivative of *G* with respect to the subpopulation strategy, we find the direction of higher per-capita growth. This is the direction evolutionary forces will drive the population, leading to

$$
\psi\_i(t) = \partial\_{\boldsymbol{\mu}} G(\boldsymbol{\mu}, \boldsymbol{\upsilon}(t), \boldsymbol{\nu}(t))|\_{\boldsymbol{\mu} = \boldsymbol{\upsilon}\_i(t)}.\tag{2}
$$

To simplify notation from here on, we define

$$
\partial\_1 G(\upsilon\_i, \upsilon, b) := \partial\_{\mathbb{H}} G(\mu, \upsilon, b)|\_{\mathbb{H} \to \upsilon\_i}.\tag{3}
$$

**Figure 1.** Schematic representation of influences on the growth rate of a population: the growth rate *G* of a subpopulation will be mainly influenced by three factors: its own strategy *vi*, the environment defined by the complete strategy set *v* and the whole population *b*.

This framework allows us to find strategy values which are resistant against others in such a way that, under the influence of natural selection, no other strategy can invade the population. These strategies are called evolutionary stable strategies (ESS) [32].

Many experiments have shown the importance of spatial arrangements when considering bacterial interaction [33,34]. As such, we aim to include spatial variables in our model. This can be done in various ways, but we concentrate on the one that was introduced by Kronik and Cohen [31]. This leads to a system of partial differential equations.

In [35], we showed the existence and uniqueness of weak and classical solutions to the associated semi-linear parabolic system in Equations (1) and (2) with Robin boundary conditions using the coupled upper-lower solution approach. In this paper, we analyse a variety of typical model functions. These translate into different long-term scenarios for different functional responses, ranging from single-strategy states to coexistence. In Section 2, we introduce the basic model and some assumptions we work with, before considering private goods models in Section 3. These models are then analysed numerically in Section 4. In Section 5, we summarise and interpret our findings.

#### **2. The Basic Model: Using the G-Function Approach**

In this case *G* takes the form:

$$\mathcal{G}(\upsilon\_i, \upsilon, b) = B(\upsilon, b) \cdot \mathcal{C}(\upsilon\_i) - \mu \|b\|\_1. \tag{4}$$

For a formulation with carrying capacity instead of population-dependent death rate, we can transform this equation into

$$\mathcal{G}(\upsilon\_i, \upsilon, b) = \mathcal{B}(\upsilon, b)\mathcal{C}(\upsilon\_i) \left(1 - \frac{\mu ||b||\_1}{\mathcal{B}(\upsilon, b)\mathcal{C}(\upsilon\_i)}\right),\tag{5}$$

which gives us *B*(*v*, *b*)*C*(*vi*)/*μ* as the capacity for the population. We can now look at the derivative of *G* with respect to *vi*.

$$
\partial\_1 G(\upsilon\_i, \upsilon, b) = B(\upsilon, b) \cdot \mathbb{C}'(\upsilon\_i). \tag{6}
$$

We make some regularity assumptions on *G* and its derivative to simplify the mathematical analysis. These do not restrict our biological application, as our model growth functions do not include discontinuous behaviour and are smooth.


These assumptions ensure that Equation (21) does not exhibit divergent behaviour (*v* → ∞ and/or *b* → ∞), which is not biologically plausible. If, e.g., *∂*1*G*(*vi*, *v*, *b*) ≥ *γ* > 0 ∀*vi*, then *c* + *λt* with *λ* = *γ* · *ε* is a lower solution of Equation (21b), as

$$
\epsilon \partial\_1 G(\lambda t + c, v, b) - \lambda \ge 0.
$$

It follows that *vi* ≥ *c* + *λt* and thus *vi* → ∞. This would mean ever-increasing strategy values without some kind of trade-off, which we do not find in nature.

In the special case of QS, we take *u* to be 0, requiring *∂*1*G*(0, *v*, *b*) = 0. This keeps *v* from leaving the biologically meaningful parameter range R<sup>+</sup> <sup>0</sup> (production cannot be lower than 0).

#### *2.1. Cost and Benefit*

To model the *G*-function for QS, one of the avenues we can take is to divide the impact of *vi*, *v*, and *b* into two parts: a growth term influenced by *vi* and a benefit provided by *v*, *b*, and possibly also *vi*. This reflects the fact that production of QS molecules is costly to the individual, while the resulting factors are public goods and therefore provide benefit to all bacteria. The additional dependence on *vi* can be seen as a form of private benefit and is discussed in detail when it occurs.

One important thing to note is that the growth term is actually *reduced* with rising *vi*, as increased public good production incurs increased metabolic costs. In this way, less energy is retained for reproduction. We denote the cost term by *<sup>C</sup>* : <sup>R</sup> <sup>→</sup> <sup>R</sup>+, and make the following assumptions:

(VI) As public good production is costly, *C* is strictly monotonically decreasing for *vi* ∈ [0, ∞).

(VII) When producing public goods, the growth rate is reduced by a certain factor:

$$0 < \mathcal{C}(v\_i) < 1 \quad \text{for } v\_i \neq 0,\tag{7}$$

While Item (VI) is clear from our assumptions on QS, Item (VII) is not as immediately clear. We introduce Item (VII) because we use this factor multiplicatively for *G*. Thus, a value of 1 would signify

unimpeded growth, while a value between 0 and 1 reduces growth. In this way, we assume that QS costs alone do not lead to negative growth rates.

The benefit of QS is provided by secreting extracellular proteins. We denote it by *<sup>B</sup>* : <sup>R</sup>*<sup>m</sup>* <sup>×</sup> <sup>R</sup>*<sup>m</sup>* <sup>→</sup> <sup>R</sup><sup>+</sup> and make two main assumptions here:

(VIII) There is a limit to how much benefit can be obtained,

$$\lim\_{v \to \infty} B(v, b) = B\_{\text{max}}.\tag{8}$$

(IX) There is no benefit if no public goods are produced,

$$B(v, b) = 0, \quad \text{if } b = 0, \sum\_{i} b\_{i} v\_{i} = 0. \tag{9}$$

Equation (8) models a saturation behaviour for the benefit—even if the cells were producing an infinite amount of extracellular protein, the benefit that can be derived is still capped through saturation of enzymes or similar phenomena. Equation (9) ensures that there is no benefit from QS when there are no living bacteria, or all of them have stopped participating in QS.

#### *2.2. Analysis*

We know from Assumption (VI) that *C*(*vi*) is strictly monotonically decreasing for positive *vi*. Combined with Assumption (V), we have *C* (*vi*) < 0 ∀*vi* > 0, *C* (0) = 0 and there exists only *vi* = 0 as an ESS. As the only stable stationary point, *vi* = 0 will attract all other solutions, leading to the decline of QS.

#### **3. Models with Private Benefit**

Following up on the *private goods*-hypothesis explained in Section 1, we now assume that there is a private benefit associated with producing the public goods, e.g. a small percentage of the produced enzymes may cling to the producing bacteria. We model this by adding a term *B*(*vi*), which is a benefit term that is solely dependent on the strategy of the subpopulation itself. This *B*(*vi*) should fulfil similar assumptions to those we made of *B*(*v*, *b*), which is why we choose to denote it by the same letter. Overall, we can write for *G*:

$$G(\upsilon\_i, \upsilon, b) = (B(\upsilon, b) + B(\upsilon\_i)) \cdot \mathbb{C}(\upsilon\_i) - \mu \|b\|\_1 \,. \tag{10}$$

It follows that

$$
\partial\_1 \mathcal{G}(\upsilon\_i, \upsilon\_i b) = \left( B(\upsilon, b) + B(\upsilon\_i) \right) \cdot \mathcal{C}'(\upsilon\_i) + \mathcal{B}'(\upsilon\_i) \cdot \mathcal{C}(\upsilon\_i), \tag{11}
$$

hence for an ESS it must hold that

$$-\mathcal{C}'(\upsilon\_i)\left(B(\upsilon,b) + B(\upsilon\_i)\right) = B'(\upsilon\_i)\mathcal{C}(\upsilon\_i).\tag{12}$$

To proceed with a more in-depth analysis, we choose a specific cost function which satisfies all our assumptions,

$$\mathbb{C}(v\_i) = \exp\left(-Kv\_i^2\right),\tag{13}$$

and look at two specific types of model terms for *B*(*vi*).

#### *3.1. First Type of Model Terms with Monotonicity Property*

The first type of model terms consists of those functions *B*(*vi*), for which *B* (*vi*)/*vi* is monotonically decreasing in *vi*. This is the case for all possible concave functions (e.g., those with *B*(0) = 0, *B*(∞) = *B*¯ > 0), as well as for

$$B(\upsilon\_i) = \frac{\upsilon\_i^2}{\upsilon\_i^2 + a^2}, \tag{10.1}$$

Equation (13) allows us to write Equation (12) as

$$\mathcal{B}\left(\mathcal{B}(\upsilon, b) + \mathcal{B}(\upsilon\_i)\right)(2\mathcal{K}) = \frac{\mathcal{B}'(\upsilon\_i)}{\upsilon\_i} =: \mathcal{R}(\upsilon\_i) \tag{14}$$

for *v*¯*<sup>i</sup>* = 0. Since we know the left-hand side of this equation to be monotonically increasing, a monotonically decreasing right-hand side means that there can only be one solution to Equation (14). Whether or not a solution exists at all can be determined by looking at the limit behaviour.

We know that *B*(*vi*) should exhibit a saturation for *vi* → ∞, forcing

$$\lim\_{\upsilon\_i \to \infty} \frac{B'(\upsilon\_i)}{\upsilon\_i} = 0.\tag{15}$$

We also take a look at the limiting values at 0:

$$\lim\_{v\_l \to 0} \left( B(v, b) + B(v\_l) \right)(2K) = 2KB(v, b) =: B\_{0\prime} \tag{16a}$$

$$\lim\_{\upsilon\_i \to 0} \frac{\mathcal{B}'(\upsilon\_i)}{\upsilon\_i} = \begin{cases} \mathcal{B}'^{\prime}(0) & \text{if } \mathcal{B}'(0) = 0 \\ +\infty & \text{if } \mathcal{B}'(0) > 0 \end{cases} \tag{16b}$$

It follows that for *B* (0) = 0 there exists a positive stationary *v*¯*<sup>i</sup>* if and only if *B* (0) > 2*KB*(*v*, *b*), while it always exists if *B* (0) > 0. If such a positive *v*¯*<sup>i</sup>* exists, it is stable (see Appendix B.2.1). If, on the other hand, *B* (0) = 0 but *B* (0) < 2*KB*(*v*, *b*), then *v*¯*<sup>i</sup>* = 0 is an ESS.

#### *3.2. Second Type of Model Terms: Hill Equations*

For a second type of model terms, we consider Hill equations. These relatively simple functions are often used in mathematical modelling for cooperative binding of molecules. In the case of QS, Dockery and Keener [36] employed them in the modelling of QS signal regulation. We look at terms of the form

$$B(\upsilon\_i) = \frac{\upsilon\_i^h}{\upsilon\_i^h + a^h} \tag{17}$$

with *h* ≥ 3, as the cases *h* = 1 and *h* = 2 fall into the first type of model terms. Then, we have

$$R(\upsilon\_i) = \frac{B'(\upsilon\_i)}{\upsilon\_i} = \frac{ha^h \upsilon\_i^{h-2}}{(\upsilon\_i^h + a^h)^2}. \tag{18}$$

Instead of being monotonically decreasing as before, *R*(*vi*) has the following properties:


For a better understanding of the ideas behind, we assume for a moment that *B*(*v*, *b*) > 0. If we can now show that *R*(*v*∗ *<sup>i</sup>* ) > 2*K*(*B*(*v*, *b*) + *B*(*v*<sup>∗</sup> *<sup>i</sup>* )), two positive stationary solutions *vi* to Equation (12) exist.

$$R(v\_i^\*) = \frac{h \cdot \left(\frac{h-2}{h+2}\right)^{\frac{h-2}{h}}}{a^2 \cdot \left(\frac{2h}{h+2}\right)^2} \qquad\qquad\qquad\qquad B(v\_i^\*) = \frac{h-2}{2h}$$

In that way, *R*(*v*∗ *<sup>i</sup>* ) > 2*K*(*B*(*v*, *b*) + *B*(*v*<sup>∗</sup> *<sup>i</sup>* )) is equivalent to the condition

$$\begin{aligned} \left(B(\upsilon, b)4h + 2(h - 2)\right)2Ka^2 &< (h - 2)^{\frac{h - 2}{h}}(h + 2)^{\frac{2}{h}}\\ \Rightarrow \quad Ka^2 &< \frac{(h - 2)^{\frac{h - 2}{h}}(h + 2)^{\frac{2}{h}}}{4(B(\upsilon, b)2h + h - 2)} \end{aligned} \tag{19}$$

which unfortunately cannot be simplified much further, without assumptions on the parameters *K* and *a*. We remind ourselves that *K* is the cost of cooperation, while *a* denotes the amount of signalling necessary to gain half of the maximal (private) benefits. Equation (19) thus gives upper limits to these terms, dependent on the steepness of the activation curve *h* as well as the public benefit *B*(*v*, *b*). A larger public benefit is actually detrimental to the existence of positive stationary strategies *vi*.

At the same time, there cannot be more than two intersections in [0, +∞) by Rolle's theorem (see, e.g., [37], p. 134), since

$$\begin{split} \left(2KB(\boldsymbol{\upsilon},\boldsymbol{b}) + B(\boldsymbol{\upsilon}\_{l})\right)' - R'(\boldsymbol{\upsilon}\_{i}) &= 2KB'(\boldsymbol{\upsilon}\_{i}) - R'(\boldsymbol{\upsilon}\_{l}) \\ &= \frac{d^{h}l\boldsymbol{\upsilon}\_{i}^{h-3}}{(\boldsymbol{\upsilon}\_{i}^{h} + a^{h})^{3}} \left(2K\boldsymbol{\upsilon}\_{i}^{2}(\boldsymbol{\upsilon}\_{i}^{h} + a^{h}) + (h+2)\boldsymbol{\upsilon}\_{i}^{h} - (h-2)a^{h}\right) \\ &= \frac{d^{h}l\boldsymbol{\upsilon}\_{i}^{h-3}}{(\boldsymbol{\upsilon}\_{i}^{h} + a^{h})^{3}} \left(2K\boldsymbol{\upsilon}\_{i}^{h+2} + (h+2)\boldsymbol{\upsilon}\_{i}^{h} + 2Ka^{h}\boldsymbol{\upsilon}\_{i}^{2} - (h-2)a^{h}\right), \end{split} \tag{20}$$

where the term in parenthesis is strictly monotonically increasing and the equation therefore only has one *ξ* ∈ (0, +∞) for which (2*KB*(*v*, *b*) + *B*(*ξ*)) − *R* (*ξ*) = 0. Both *v*¯*<sup>i</sup>* = 0 and the larger of both positive stationary solutions are stable regardless of parameter values (see Appendix B.2.1).

If *B*(*v*, *b*) = 0, which occurs if for example the existing subpopulations do not engage in QS, the situation is different. As both *R*(0) = 0 and 2*KB*(*v*, *b*) + *B*(0) = 0, there can only be one positive intersection. By using Property 5, we know that *R*(*vi*) is increasing more quickly than 2*KB*(*v*, *b*) + *B*(*vi*) for small *vi*. It follows that a positive intersection exists. As we know that for *B*(*v*, *b*) = 0 there is an intersection at *vi* = 0 and from Equation (20) that there can only be two intersections in [0, +∞), there can be no other positive intersection. Additionally, one can prove this positive intersection to be stable, while the zero solution is unstable (see Appendix B.2.1).

In biological terms, our results indicate that the surrounding subpopulations have a profound impact on the evolution of cooperativity. If one of the subpopulations is cooperating (*B*(*v*, *b*) = 0), the others will experience bistable behaviour—they might cooperate themselves or not participate in QS, depending on starting strategy. If, on the other hand, all of the subpopulations do not cooperate at first (*B*(*v*, *b*) = 0), then the evolutionary pressure will drive them towards the positive ESS, meaning they pick up QS with time.

#### **4. Numerical Simulations**

The full set of equations we are working with therefore reads

$$\dot{b}\_{l}(\mathbf{x},t) = G(v\_{l}(\mathbf{x},t), v(\mathbf{x},t), b(\mathbf{x},t)) \cdot b\_{l}(\mathbf{x},t) + D \,\triangle b\_{l}(\mathbf{x},t) \tag{21a}$$

$$\psi\_{i}(\mathbf{x},t) = \varepsilon \partial\_{1} \mathbb{G}(v\_{i}(\mathbf{x},t), v(\mathbf{x},t), b(\mathbf{x},t)) + D \left( \frac{2 \nabla b\_{i}(\mathbf{x},t) \cdot \nabla v\_{i}(\mathbf{x},t)}{b\_{i}(\mathbf{x},t)} + \triangle v\_{i}(\mathbf{x},t) \right), \tag{21b}$$

noting that ˙ *bi*(*x*, *<sup>t</sup>*) = *<sup>∂</sup>bi ∂t* .

We used our own experimental data from [34] to derive values for the parameters in our model, which can be found in Table A1. Our aim here was to obtain rough estimates which replicate the general qualitative behaviour, as opposed to fit the experimental data perfectly. To solve the system numerically, we employed the method of lines [38] and used an explicit Runge–Kutta solver on the discretised system. The simulations were run in Matlab.

We assumed a closed system in which the existing subpopulations are mixed but not completely homogeneous in the beginning. This leads to Neumann boundary conditions and the initial distribution displayed in Figure 2.

**Figure 2.** Visual representation of initial conditions for numerical simulations. For initial conditions, we used patches. The two subpopulations (dark and light patches) were mixed but not homogeneously distributed.

We also added in a feature of QS that we have neglected before: the costs of QS depend on the signal level, as QS products are only produced when the signal level is high enough. We can achieve this by modifying the cost term from Equation (13) to

$$\mathcal{K}(\upsilon\_i, \upsilon, b) = \exp\left(-\mathbb{K} \cdot (\mathcal{B}(\upsilon, b) \cdot \upsilon\_i)^2\right),\tag{22}$$

as *B*(*v*, *b*) is a measure for the level of QS. We have not done so before, since it does not change the analysis qualitatively, but leads to unnecessarily long and obfuscating terms.

#### *4.1. Basic Model*

We have seen in our calculations that *vi* = 0 is the only stable strategy in this situation. We should therefore have both *vi* → 0 and *bi* → 0. However, when looking at the simulation results for 3000 h as displayed in Figure 3, we notice that this convergence is very slow indeed. A quick calculation of *∂*1*G* and *G* for our parameter values confirms that both are close to zero. Thus, even though QS is theoretically unstable, both cheaters and producers persist alongside each other for a very long time, albeit at different densities.

(**a**) Population dynamics (**b**) Strategy dynamics **Figure 3.** Long-term behaviour of two populations with different start strategies, using a *G*-function without additions. Both populations numbers are slowly converging in time course towards zero, with producers (dark line) having the lower numbers.

#### *4.2. Models with Private Benefit*

In Section 3 we distinguish between two types of model functions. Depending on the Hill coefficient, Hill-terms can fall into either one. We investigated the long-term behaviour for *h* = 1, 2 or 3.

As predicted, the behaviour in these cases is quite different. For *h* ≤ 2, *B*(*vi*) falls into the first type of model functions. For this type, we postulate that there can only be one positive stationary point for *v* and if it is stable, and then the zero solution is unstable. This is the case for the parameters we calculated. Additionally, we found that for *h* = 1 the zero solution is not a stationary point. Thus, both strategies converge towards the stable positive equilibrium, one from above and the other from below. Population 2, which started out as cheaters, gains QS functionality, while there is reduced production from Population 1. Since reducing production is slower to reach the stable point in this parameter constellation, Population 1 succumbs to the population pressure from Population 2 and dies out (see Figure 4). This happens on a rather short time frame of less than 200 h.

(**b**) Strategy dynamics

**Figure 4.** Evolution of two populations with different start strategies, using a *G*-function with private benefit and hill factors of *h* = 1 (line), 2 (dashed) and 3 (dash-dotted). For *h* = 1, Population 1 (dark colour) dies out rather quickly, while Population 2 (light colour) gains the QS functionality, albeit at a low level, and remains at a stable population level. If *h* = 2, Population 1 reduces its QS strategy to a lower, but stable value. Population 2 is unable to compete and dies out. For *h* = 3, Populations 1 and 2 coexistt in a stable way with similar numbers. One population is QS active while the other is not.

For *h* = 2, we can see in Figure 4 that once again the strategies converge to a positive value. However, this time the zero solution is a stationary point, although an unstable one. Hence, Population 2 cannot gain QS functionality by starting out with *v*<sup>2</sup> = 0. The end result is the extinction of Population 2, while Population 1 remains stable.

Overall, for *h* ≤ 2, we found that one population is driven to extinction, while the other stays at a stable level with QS intact at lower levels.

When *h* > 2, *B*(*vi*) falls into the second type of model functions. That means there might be a stable positive strategy in addition to a stable zero solution. We found this to be the case for *h* = 3 for our parameters, as one can see in Figure 4. This means that Populations 1 and 2 remain at stable population and strategy levels, with Population 1 being QS active while Population 2 consists of cheaters. In this scenario, producers and non-producers can live side-by-side indefinitely.

#### **5. Discussion**

We have applied the *G*-function ansatz to the QS dynamics of *P. aeruginosa* and analysed the resulting model both analytically and numerically. We have found that the type of *G*-function considered changes the expected outcome drastically.

If we assume that there is no private benefit to QS at all, then the long-term outcome is preassigned and independent of parameters. Whether QS is inactive in all bacteria after some time, or the subpopulation of secretors dies out, the effect remains the same: QS is evolutionary unstable.

As soon as we assume that there is some kind of private benefit, the situation changes. Depending on parameter values, QS might still be unstable, but there can also be positive evolutionary stable secretion values. For the first type of model functions, there is only one stable secretion value. It follows that in the long term the population will unitise (Figure 4). The situation remains much the same if, instead of a direct benefit, we consider a reduction in death rate (results not shown). In contrast, private benefit functions of the second type admit bistable behaviour. This allows secretors and cheaters to coexist indefinitely in addition to single subpopulation states. We have also found that bacteria are driven towards QS if there is no pre-existing cooperation in their environment.

Since the left-hand-side of Equation (12) has a value of zero for *vi* = 0, whereas the right-hand-side has a value greater than zero for *vi* = 0, there could be any number of stable strategies (or none) if one chooses different functional responses.

Our model functions have simplified the actual process of QS a great deal, but one could also include direct dependencies on the signal molecules and enzymes produced. Such models with abiotic components were considered by Cohen et al. [30]. While their short-term dynamics might differ from the scenarios we have considered, their long-term outcome does not. They are also difficult to analyse without resorting to numerical methods, which is why we have concentrated on the simplified functions.

The *G*-function has allowed us to determine the long-term evolution of QS by considering the short-term population dynamics. By analysing different functional responses, we have explored some possible outcomes. One could also further refine the population dynamics for different real world scenarios. The *G*-function serves to always model the resulting evolutionary pressure, which tells us something about the way mutation and selection are going to drive the population. As both are happening on a fast time-scale in bacteria, they have a big impact on the interaction between bacterial populations as well as on their interaction with humans (see, e.g., [3]). The *G*-function can thus be a valuable tool in modelling bacterial interactions.

For completeness, in Appendix C, we show the spatial evolution of two populations with different start strategies, using a G-function without additions (Figure A1) and with private benefit and a hill factor of *h* = 1 (Figure A2).

Lastly, we would like to refer the interested reader to the following associated publications which complement this study: Mund et al. [34] and Mund et al. [35].

**Author Contributions:** Conceptualization, C.K. and J.P.-V.; Data curation, A.M.; Formal analysis, A.M.; Funding acquisition, C.K.; Investigation, A.M., C.K. and J.P.-V.; Methodology, A.M.; Software, A.M.; Supervision, C.K. and J.P.-V.; Validation, A.M.; Visualization, A.M.; Writing—original draft, A.M., C.K. and J.P.-V.; Writing—review and editing, A.M., C.K. and J.P.-V.

**Funding:** This work was supported by the German Research Foundation (DFG) and the Technical University of Munich within the funding programme Open Access Publishing.

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

#### **Appendix A. Variables and Parameters**


**Table A1.** Table of all occurring parameters.

#### **Appendix B. Proofs**

*Appendix B.1. Properties of R*(*Vi*) *for the Second Type of Model Terms*

We first prove that *R*(*vi*) has exactly one maximum as well as calculate its maximal value. To that end, we take the derivative of *R*

$$\mathcal{R}'(v\_l) = \frac{a^h h v\_l^{h-3} \left(-(h+2)v\_l^h + (h-2)a^h\right)}{\left(v\_l^h + a^h\right)^3} \tag{A1}$$

and search for critical points

$$0 = d^h l \upsilon\_i^{h-3} \left( -(h+2)\upsilon\_i^h + (h-2)d^h \right). \tag{A2}$$

For *h* > 3, there exists a critical point at *vi* = 0. However, since *R*(0) = 0 and *R*(*vi*) > 0 for *vi* > 0, this critical point is clearly a minimum on the bounded space *vi* <sup>∈</sup> <sup>R</sup><sup>+</sup> <sup>0</sup> . The other critical point satisfies

$$v\_i^\* = \sqrt[k]{\frac{h-2}{h+2}} \cdot a\_i$$

and is thus unique in R<sup>+</sup> <sup>0</sup> . It remains to show that *v*<sup>∗</sup> *<sup>i</sup>* is indeed a maximum. However, knowing that *R*(0) = 0 and lim*vi*→<sup>∞</sup> *R*(*vi*) = 0, we can conclude that it can only be a maximum.

We prove Property 5 by induction.

**base case** We can derive from the definition of *R*(*vi*) that

$$\begin{aligned} R'(v\_i) &= \frac{a^h h v\_i^{h-3} \left( -(h+2) v\_i^h + (h-2) a^h \right)}{\left( v\_i^h + a^h \right)^3} \\ &= \frac{v\_i^{h-(1+2)} f(v\_i)}{\left( v\_i^h + a^h \right)^{1+2}} \end{aligned}$$

where *f*(*vi*) = *ahh* −(*<sup>h</sup>* + <sup>2</sup>)*v<sup>h</sup> <sup>i</sup>* + (*<sup>h</sup>* − <sup>2</sup>)*a<sup>h</sup>* and thus *<sup>f</sup>*(0) = *ha*2*h*(*<sup>h</sup>* − <sup>2</sup>) > 0.

**inductive step** Assuming that the statement holds for *<sup>R</sup>*(*n*)(*vi*) and that *<sup>n</sup>* <sup>+</sup> <sup>1</sup> <sup>≤</sup> *<sup>h</sup>* <sup>−</sup> 2, we have

$$\begin{split} R^{(n+1)}(\upsilon\_{i}) &= \left( (\upsilon\_{i}^{h} + a^{h})^{n+2} (h - (2+n)) \upsilon\_{i}^{h - (n+3)} f(\upsilon\_{i}) \right. \\ &\left. + (\upsilon\_{i}^{h} + a^{h})^{n+2} f'(\upsilon\_{i}) \upsilon\_{i}^{h - (n+2)} \right. \\ &\left. - (n+2) h (\upsilon\_{i}^{h} + a^{h})^{n+1} f(\upsilon\_{i}) \upsilon\_{i}^{2h - (n+3)} \right) / \left( \upsilon\_{i}^{h} + a^{h} \right)^{2n+4} \\ &= \frac{\upsilon\_{i}^{h - (n+1+2)} g(\upsilon\_{i})}{(\upsilon\_{i}^{h} + a^{h})^{n+1+2}} . \end{split}$$

where

$$g(v\_i) = (v\_i^h + a^h) \left( (h - (2 + n))f(v\_i) + f'(v\_i)v\_i \right) - (n + 2)hf(v\_i)v\_i^h.$$

and thus *<sup>g</sup>*(0) = *<sup>a</sup>h*(*<sup>h</sup>* − (<sup>2</sup> + *<sup>n</sup>*))*f*(0) > 0 since *<sup>n</sup>* < *<sup>h</sup>* − 2 and *<sup>f</sup>*(0) > 0.

It follows that *<sup>R</sup>*(*n*)(0) = 0 for *<sup>n</sup>* <sup>&</sup>lt; *<sup>h</sup>* <sup>−</sup> 2 and *<sup>R</sup>*(*h*−<sup>2</sup>)(0) <sup>&</sup>gt; 0. At the same time, a similar argument shows that

$$\left(2K(\mathcal{B}(v,b) + \mathcal{B}(v\_i))\right)^{(n)} = 2KB^{(n)}(v\_i) = \frac{\upsilon\_i^{h-n} \mathfrak{L}Kf(\upsilon\_i)}{(\upsilon\_i^h + a^h)^{1+n}}.$$

As such, <sup>2</sup>*K*(*B*(*v*, *<sup>b</sup>*) + *<sup>B</sup>*(*vi*)) (*n*) = 0 for *n* < *h*.

*Appendix B.2. Stability of v*¯*<sup>i</sup>*

Appendix B.2.1. Models with Private Benefit

#### First Type of Model Terms

For *B*(*vi*) in the first type of model terms, we assume that *B* (*vi*)/*vi* is monotonically decreasing. It follows that

$$0 > \left(\frac{B'(v\_i)}{v\_i}\right)' = \frac{v\_i B''(v\_i) - B'(v\_i)}{v\_i^2}$$

$$\implies \quad 0 > B''(v\_i) - \frac{B'(v\_i)}{v\_i} \quad \forall v\_i > 0. \tag{A3}$$

Looking at the second derivative of *G*, we have

$$\begin{split} \partial\_1^2 G(v\_l, v, b) &= \left( B(v, b) + B(v\_l) \right) \mathcal{C}''(v\_l) + 2B'(v\_l) \mathcal{C}'(v\_l) + B''(v\_l) \mathcal{C}(v\_l) \\ &= \left( \left( B(v, b) + B(v\_l) \right) \left( -2K + 4K^2 v\_l^2 \right) + 2B'(v\_l)(-2K)v\_l + B''(v\_l) \right) \mathcal{e}^{-Kv\_l^2} \\ \partial\_1^2 G(0, v, b) &= B(v, b)(-2K) + \mathcal{B}''(0), \end{split}$$

from which we can gather that *v*¯*<sup>i</sup>* = 0 is stable if *B* (0) < 2*KB*(*v*, *b*). We use Equation (14) and substitute for *B*(*v*, *b*) + *B*(*v*¯*i*)

$$\partial\_1^2 G(\vartheta\_l, \upsilon\_i b) = \left(\frac{B'(\vartheta\_i)}{2K\tilde{\upsilon}\_i} \left(-2K + 4K^2 \vartheta\_l^2\right) + 2B'(\vartheta\_l)(-2K)\upsilon\_i + B''(\upsilon\_i)\right) \mathbf{e}^{-K\vartheta\_i^2} $$

$$= \left(-\frac{B'(\vartheta\_i)}{\vartheta\_i} + B''(\vartheta\_i) - 2K\vartheta\_l B'(\vartheta\_i)\right) \mathbf{e}^{-K\tilde{\upsilon}\_l^2}.\tag{A4}$$

As *B* (*vi*) > 0, stability of positive *v*¯*<sup>i</sup>* immediately follows from Equation (A3).

Second Type of Model Terms

We start out by assuming *B*(*v*, *b*) > 0. We know that *vi* = 0 is stable if *B* (0) < 2*KB*(*v*, *b*) and that for *h* ≥ 3 it holds that *B* (0) = 0. Hence, *vi* = 0 is an ESS regardless of parameter values.

A positive stationary solution *v*¯*<sup>i</sup>* is stable, if

$$\begin{aligned} 0 &> -\left(\frac{1}{\vec{v}\_i} + 2K\boldsymbol{\sigma}\_i\right)B'(\boldsymbol{\sigma}\_i) + B''(\boldsymbol{\sigma}\_i) \\ \Rightarrow \quad 0 &> -\left(\frac{1}{\vec{v}\_i} + 2K\boldsymbol{\sigma}\_i\right)\frac{ha^h\boldsymbol{\sigma}\_i^{h-1}}{(\boldsymbol{\sigma}\_i^h + a^h)^2} + \frac{ha^h\boldsymbol{\sigma}\_i^{h-2}\left(-(h+1)\boldsymbol{\sigma}\_i^h + (h-1)a^h\right)}{(\boldsymbol{\sigma}\_i^h + a^h)^3} \\ \Rightarrow \quad 0 &> -\left(\frac{1}{\vec{v}\_i} + 2K\boldsymbol{\sigma}\_i\right)\boldsymbol{\sigma}\_i(\boldsymbol{\sigma}\_i^h + a^h) - (h+1)\boldsymbol{\sigma}\_i^h + (h-1)a^h. \end{aligned}$$

We can rearrange this inequality to

$$
\left(\frac{1}{\mathcal{O}\_i} + 2K\vartheta\_i\right)\vartheta\_i(\vartheta\_i^h + a^h) + (h+1)\vartheta\_i^h > (h-1)a^h.
$$

The larger of both positive stationary solutions fulfils *v*¯*<sup>i</sup>* > *v*<sup>∗</sup> *<sup>i</sup>* = *<sup>h</sup> h*−<sup>2</sup> *<sup>h</sup>*+<sup>2</sup> · *<sup>a</sup>*, if it exists. It follows that

$$\begin{aligned} \left(\frac{1}{\mathcal{O}\_l} + 2K\phi\_l\right)\vartheta\_i(\vartheta\_i^h + a^h) + (h+1)\vartheta\_i^h\\ &> \left(1 + 2K\vartheta\_i^2\right)\left(\frac{h-2}{h+2}\cdot a^h + a^h\right) + (h+1)\cdot \frac{h-2}{h+2}\cdot a^h\\ &> \frac{2h}{h+2}\cdot a^h + \frac{(h+1)(h-2)}{h+2}\cdot a^h\\ &= \frac{2h + h^2 - h - 2}{h+2}\cdot a^h\\ &= (h-1)\cdot a^h. \end{aligned}$$

Hence, the larger of both positive solutions is always stable, if it exists. It follows immediately from the stability of the zero solution that the smaller of both positive solutions must be unstable.

In the other case, we first note that both *B* (0) = 0 and *B*(*v*, *b*) = 0, which means we cannot gain insight into the stability of the zero solution directly through Equation (A4). However, we know that, for small enough *ε* > 0, *R*(*ε*) > 2*KB*(*ε*) and that there exists only one positive *ξ* for which 2*KB* (*ξ*) = *R* (*ξ*). Additionally, we know that *ξ* lies in the open interval between zero and the positive intersection *v*¯*<sup>i</sup>* and as such is unequal to *v*¯*i*. It can thus follow that

$$R'(\boldsymbol{\upsilon}\_{i}) < 2KB'(\boldsymbol{\upsilon}\_{i})$$

$$R'(\boldsymbol{\upsilon}\_{i}) = \left(\frac{B'(\boldsymbol{\upsilon}\_{i})}{\boldsymbol{\upsilon}\_{i}}\right)' = \frac{B''(\boldsymbol{\upsilon}\_{i})}{\boldsymbol{\upsilon}\_{i}} - \frac{B'(\boldsymbol{\upsilon}\_{i})}{\boldsymbol{\upsilon}\_{i}^{2}}$$

$$\Rightarrow \quad \frac{B''(\boldsymbol{\upsilon}\_{i})}{\boldsymbol{\upsilon}\_{i}} - \frac{B'(\bar{\boldsymbol{\upsilon}}\_{i})}{\bar{\boldsymbol{\upsilon}}\_{i}^{2}} < 2KB'(\boldsymbol{\upsilon}\_{i})$$

$$-2KB'(\boldsymbol{\upsilon}\_{i})\boldsymbol{\upsilon}\_{i} - \frac{B'(\boldsymbol{\upsilon}\_{i})}{\boldsymbol{\upsilon}\_{i}} + B''(\boldsymbol{\upsilon}\_{i}) < 0,$$

which is exactly the condition for stability from Equation (A4). For *B*(*v*, *b*) = 0, we thus have an unstable zero solution and a stable positive strategy.

#### **Appendix C. Spatial Evolution of Two Populations with Different Start Strategies**

In this section, we show the spatial evolution of two populations with different start strategies, using a G-function without additions (Figure A1) and with private benefit and a hill factor of *h* = 1 (Figure A2).

**Figure A1.** Long-term behaviour of two populations with different start strategies, using a *G*-function without additions.

**Figure A2.** Evolution of two populations with different start strategies, using a *G*-function with private benefit and a hill factor of *h* = 1. Note the shorter timescale in this plot.

#### **References**



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