*Article* **Cross Diffusion Induced Turing Patterns in a Tritrophic Food Chain Model with Crowley-Martin Functional Response**

#### **Nitu Kumari \* and Nishith Mohan**

School of Basic Sciences, Indian Institute of Technology Mandi, Mandi, Himachal Pradesh 175001, India; nishithmohan888@gmail.com

**\*** Correspondence: nitu@iitmandi.ac.in

Received: 20 January 2019; Accepted: 23 February 2019; Published: 1 March 2019

**Abstract:** Diffusion has long been known to induce pattern formation in predator prey systems. For certain prey-predator interaction systems, self diffusion conditions ceases to induce patterns, i.e., a non-constant positive solution does not exist, as seen from the literature. We investigate the effect of cross diffusion on the pattern formation in a tritrophic food chain model. In the formulated model, the prey interacts with the mid level predator in accordance with Holling Type II functional response and the mid and top level predator interact via Crowley-Martin functional response. We prove that the stationary uniform solution of the system is stable in the presence of diffusion when cross diffusion is absent. However, this solution is unstable in the presence of both self diffusion and cross diffusion. Using a priori analysis, we show the existence of a inhomogeneous steady state. We prove that no non-constant positive solution exists in the presence of diffusion under certain conditions, i.e., no pattern formation occurs. However, pattern formation is induced by cross diffusion because of the existence of non-constant positive solution, which is proven analytically as well as numerically. We performed extensive numerical simulations to understand Turing pattern formation for different values of self and cross diffusivity coefficients of the top level predator to validate our results. We obtained a wide range of Turing patterns induced by cross diffusion in the top population, including floral, labyrinth, hot spots, pentagonal and hexagonal Turing patterns.

**Keywords:** cross diffusion; Turing patterns; non-constant positive solution

#### **1. Introduction**

The two primary objectives of studying systems ecology are to get an understanding of the dynamics of the ecological systems and the nature of the forces, which determine the community structure. Earlier, the classical ecological models focused on two species interactions. Continuous time models of two interacting species have been studied extensively in the literature (see [1] and references therein). Mathematically, only two basic patterns are exhibited in these models—approach to a steady state or to a limit cycle [2]—but the ecological models of nature exhibit very complex dynamical behavior. Price et al. [3] proposed that community behavior must be based on three or more trophic levels. Three species continuous time models have observed more complex dynamical behavior [4–7]. All these studies depend upon the classical types of functional responses, which include Holling Types I, II and III and Holling-Tanner ratio dependent functional response.

The Crowley-Martin [8] functional response is one of the predator dependent functional responses, i.e., the functional response is function of both prey and predator abundance because of predator interference. The basic assumption behind the formulation of Crowley-Martin functional response is that there is a decrease in predator feeding rate because of high predator density even if there is a high prey density. Hence, the effect of predator interference on the feeding rate is important, even when an

individual predator is handling or searching for a prey at a given instant of time. Here, *v* and *r* are population density of prey and predator, respectively. The parameters *w*2, *d*, *b* are positive parameters that describe the effects of capture rate, handling time and magnitude of interference among predators on the feeding rate, respectively. The per capita rate in this formulation is given by:

$$f(v, r) = \frac{w\_2 v}{1 + dv + br + bdvr}.$$

The Crowley-Martin functional response is used for datasets that indicate an asymptotic feeding rate that is affected by predator density. A continuous time model to analyze the dynamics of a three species food chain with Crowley-Martin functional response was studied by Upadhyay et al. [9]. The model system in [9] is based on the assumption that the species under consideration are spatially homogeneously distributed; however, in nature, the species distribution is always inhomogeneous. Therefore, to model a realistic food chain scenario, reaction diffusion mechanism should be considered, as employed in [10–12]. A qualitative analysis of a two species model with Crowley-Martin functional response and diffusion is carried out in [13]. In the above model, the movement of a species is determined solely by their own characteristics, i.e., the movements of the species in these models are physically affected by the population pressures due to the mutual interference among the individuals of the same species. Therefore, such models take into account only the self diffusion of the species in concern. However, in the case of prey-predator interacting systems, the movement of one species may affect the movement or motility of the other. The above models fail to describe such situation, which is more realistic. In many systems, the predators develop migratory strategies to take advantage over the prey. In the case of a three species prey-predator model, such a migratory behavior depends on the concentration of both predators (i.e., the mid level predator and the top level predator). This leads to *cross diffusive system*, in addition to each species natural tendency to diffuse i.e., self diffusion; such models are studied in [14–18]. The cross diffusive coefficient may have positive or negative values. The positive cross diffusivity coefficient represents that one species tends to move in the direction of lower concentration of another species. The negative cross diffusion term represents that the population flux of one species is in the direction of higher concentration of the other species.

In the present paper, we consider a tritrophic food chain model in which the species interact via Crowley-Martin and Holling type II functional response at different levels. We introduce a nonlinear cross diffusion among top level and mid level predators. The objective is to understand the effect of both types of diffusion self and cross on the stability of the system. An analytic approach is adopted to understand the formation of inhomogeneous steady state solutions. Numerical simulation depicting spatial pattern formation was performed for a wide range of parameters to understand the predation behavior due to the effect of cross and self diffusion.

In Section 2, we formulate the temporal model using various assumptions and spatially extend it to formulate a nonlinear reaction diffusion system and with cross diffusivity (a strongly coupled parabolic system). In Section 3, we perform a linear stability analysis of the temporal model and obtain the conditions under which the stationary uniform solution is locally asymptotically stable. In Section 4, we provide analytical proof of stability of the stationary uniform solution in the presence of diffusion and in the absence of cross diffusion but unstable in the presence of cross diffusion. In Section 5, we provide an analytical proof of the existence of inhomogeneous steady states. Section 6 deals with the numerical simulation of the cross diffusive system and explains Turing pattern formation for different values of cross diffusive and self diffusive coefficients of the top level predator.

#### **2. Model System**

We first consider a temporal model governing the dynamics of tritrophic food chain, consisting of Holling type-II and Crowley Martin type functional responses, defined by a system of differential equations [9]. Here, *U*(*T*) is the population density of the lowest trophic level species (prey), *V*(*T*) is the population density of the middle trophic level species (intermediate predator) and *R*(*T*) is the

population density of the highest trophic level species (top predator), at any time *T*. We made the following basic assumptions for formulation of the proposed model system.

**Assumption 1.** *The prey U grows with an intrinsic growth rate a*1*, and has a carrying capacity of K in the absence of predator V. D is a measure of the extent to which the environment can provide protection to U and w is the maximum value which per capita reduction rate of U can attain. In addition, the mid level predator V predates on <sup>U</sup> in accordance with Holling Type II functional response, wUV <sup>U</sup>* <sup>+</sup> *<sup>D</sup>. Therefore, considering the above assumptions, we formulate the first equation of the model system as:*

$$\frac{d\mathcal{U}}{dT} = a\_1 \mathcal{U} \left( 1 - \frac{\mathcal{U}}{K} \right) - \frac{w \mathcal{U} \mathcal{V}}{\mathcal{U} + D} \tag{1}$$

**Assumption 2.** *The intermediate or mid level predator V has a natural death rate of a*2*. w*<sup>1</sup> *describes the maximum value which per capita reduction rate of V can attain, D*<sup>1</sup> *is similar to D of U and p represents the internal competition coefficient among the members of V. The predation of V over U is governed by Holling Type II functional response. The top level predator R predates on V in accordance with Crowley-Martin functional response. The constants w*2*, b and d are the saturating Crowley-Martin type functional response parameters, where w*<sup>2</sup> *describes capture rate effect, b represents handling time of prey and d is the magnitude of interference among predators. The second equation of the model system is given as:*

$$\frac{dV}{dT} = -a\_2 V - pV^2 + \frac{w\_1 LV}{U + D\_1} - \frac{w\_2 VR}{1 + dV + bR + bdVR} \tag{2}$$

**Assumption 3.** *The top level predator R dies at a natural death rate of c. As the predating behavior of R over V is described by the Crowley-Martin functional response. w*<sup>3</sup> *is the saturating Crowley-Martin type functional response parameter similar to w*2*. The third equation of the model system takes the following form:*

$$\frac{dR}{dT} = -cR + \frac{w\_3 V R}{1 + dV + bR + bdVR} \tag{3}$$

All parameters described above take only positive values and the model system is a three species food chain model involving a hybrid type of prey dependent and predator dependent functional response. Therefore, using the above assumptions, the model system takes the following form:

$$\begin{aligned} \frac{d\mathcal{U}}{dT} &= a\_1 \mathcal{U} \left( 1 - \frac{\mathcal{U}}{K} \right) - \frac{w \mathcal{U} \mathcal{V}}{l \mathcal{I} + D} \\ \frac{d\mathcal{V}}{dT} &= -a\_2 \mathcal{V} - p\mathcal{V}^2 + \frac{w\_1 \mathcal{U} \mathcal{V}}{\mathcal{U} + D\_1} - \frac{w\_2 \mathcal{V} \mathcal{R}}{1 + dV + b\mathcal{R} + bd \mathcal{V} \mathcal{R}} \\ \frac{dR}{dT} &= -c\mathcal{R} + \frac{w\_3 \mathcal{V} \mathcal{R}}{1 + dV + b\mathcal{R} + bd \mathcal{V} \mathcal{R}} \end{aligned} \tag{4}$$

The model system described by Equation (4) consists of thirteen parameters, which makes the analysis quite cumbersome, therefore the number of parameters are reduced by rescaling the model system. The model is rescaled, using the following new variables and parameters:

$$\begin{aligned} t &= a\_1 T, \quad u = \frac{U}{K}, \quad v = \frac{wV}{a\_1 K}, \quad r = \frac{ww\_2 R}{a\_1^2 dK}, \quad w\_4 = \frac{D}{K}, \quad w\_5 = \frac{a\_2}{a\_1}, \quad w\_6 = \frac{w\_1}{a\_1}, \quad w\_7 = \frac{D\_1}{K}, \\\ w\_8 &= \frac{a\_1 b}{w\_2}, \quad w\_9 = \frac{a\_1^2 b dK}{w w\_2}, \quad w\_{10} = \frac{w}{a\_1 dK}, \quad w\_{11} = \frac{c}{a\_1}, \quad w\_{12} = \frac{w\_3}{a\_1 d}, \quad c = \frac{pw}{a\_1^2 K} \end{aligned} \tag{5}$$

The rescaled system is as follows:

$$\begin{aligned} \frac{du}{dt} &= u\left(1 - u - \frac{v}{u + w\_4}\right) \\ \frac{dv}{dt} &= v\left(-w\_5 + \frac{w\_6 u}{u + w\_7} - cv - \frac{r}{v + \left(w\_8 + w\_9 v\right)r + w\_{10}}\right) \\ \frac{dr}{dt} &= r\left(-w\_{11} + \frac{w\_{12} v}{v + \left(ws + w\_9 v\right)r + w\_{10}}\right) \end{aligned} \tag{6}$$

The model system proposed in Equation (4) and the rescaled form of the model given in Equation (6) are based on the assumption that the species under consideration are homogeneously distributed but in nature the species distribution is always inhomogeneous. Therefore, to model a realistic food chain scenario, we consider the model system with diffusion. Consider the spatially explicit three species predator prey food chain model system. At any location (*x*, *y*) and time *t*, the interaction of three species populations *u*(*x*, *y*, *t*), *v*(*x*, *y*, *t*) and *r*(*x*, *y*, *t*) can be modeled with the reaction diffusion equations given in Equation (7). Here, Δ denotes two dimensional Laplacian operator given by <sup>Δ</sup> <sup>≡</sup> *<sup>∂</sup>*<sup>2</sup> *<sup>∂</sup>x*<sup>2</sup> <sup>+</sup> *∂*2 *∂y*<sup>2</sup> , where, *<sup>x</sup>*, *<sup>y</sup>* <sup>∈</sup> <sup>Ω</sup> and *<sup>t</sup>* <sup>&</sup>gt; 0. <sup>Ω</sup> is a bounded domain in <sup>R</sup><sup>2</sup> with a smooth boundary *∂*Ω.

$$\begin{aligned} \frac{\partial u}{\partial t} - d\_1 \Delta u &= u \left( 1 - u - \frac{v}{u + w\_4} \right) \\ \frac{\partial v}{\partial t} - d\_2 \Delta v &= v \left( -w\_5 + \frac{w\_6 u}{u + w\gamma} - \varepsilon v - \frac{r}{v + (w\_8 + w\gamma v)} r + w\_{10} \right) \\ \frac{\partial r}{\partial t} - d\_3 \Delta r &= r \left( -w\_{11} + \frac{w\_{12} v}{v + \left( w\_8 + w\_9 v \right) r + w\_{10}} \right) \\ \frac{\partial u}{\partial n} = \frac{\partial v}{\partial n} &= \frac{\partial r}{\partial n} = 0. \end{aligned} \tag{7}$$

The boundary conditions *<sup>∂</sup><sup>u</sup> <sup>∂</sup><sup>n</sup>* <sup>=</sup> *<sup>∂</sup><sup>v</sup> <sup>∂</sup><sup>n</sup>* <sup>=</sup> *<sup>∂</sup><sup>r</sup> <sup>∂</sup><sup>n</sup>* <sup>=</sup> 0 are directional derivative normal to the boundary and *n* is the outward normal to Ω.

In the above system, the movement of a species is determined solely by their own characteristics, i.e., the movements of the species can be physically affected by the population pressures due to the mutual interference between the individuals of the same species. Therefore, we refer to the constants *d*1, *d*<sup>2</sup> and *d*<sup>3</sup> as the *self diffusion rates* of species *u*, *v* and *r*, respectively. However, in the case of a prey-predator interacting system, the movement of one species may affect the movement or motility of the other. Hence, the above model cannot be used to describe such a situation. There are chances of development of migratory strategies among the predators to take advantage over the prey. This behavior can be described by taking into account the effect of *cross diffusion*. Along with the self tendency to move, i.e., self diffusion, the species also migrate or move under the influence of other. Such behavior takes into account concentration of both predators—the mid and the top level predators—constituting a cross diffusive system. The coefficient of cross diffusion may be positive or negative. The positive cross diffusion coefficient represents that one species tends to move in the direction of lower concentration of another species. The negative cross diffusion term represents that the population flux of one species is in the direction of higher concentration of the other. Therefore, we propose the interaction of the above three species population model with cross diffusion as follows:

$$\begin{aligned} \frac{\partial u}{\partial t} - d\_1 \Delta u &= u \left( 1 - u - \frac{v}{u + w\_4} \right) \\ \frac{\partial v}{\partial t} - d\_2 \Delta v &= v \left( -w\_5 + \frac{w\_6 u}{u + w\_7} - cv - \frac{r}{v + (w\_8 + w\_9 v)} r + w\_{10} \right) \\ \frac{\partial r}{\partial t} - d\_3 \Delta \left( r + d\_4 v r \right) &= r \left( -w\_{11} + \frac{w\_{12} v}{v + (w s + u y v)} r + w\_{10} \right) \\ \frac{\partial u}{\partial n} = \frac{\partial v}{\partial n} = \frac{\partial r}{\partial n} = 0. \end{aligned} \tag{8}$$

where <sup>Δ</sup> denotes two dimensional Laplacian operator given by <sup>Δ</sup> <sup>≡</sup> *<sup>∂</sup>*<sup>2</sup> *<sup>∂</sup>x*<sup>2</sup> <sup>+</sup> *∂*2 *∂y*<sup>2</sup> , (*x*, *y*) ∈ Ω.

The nonlinear diffusion terms in the equation governing the dynamics of top level predator *r* implies that the direction of dispersion of top predator contains a self diffusion term by which it moves from a region of its higher concentration to a region of its lower concentration and a cross diffusive term. The top predator *r* diffuses with the flux

$$\mathbf{J} = -\nabla \left( d\_3 r + d\_3 d\_4 \upsilon r \right) = -d\_3 d\_4 r \nabla \upsilon - \left( d\_3 + d\_3 d\_4 \upsilon \right) \nabla r.$$

where the −*d*3*d*4∇*r* < 0, the −*d*3*d*4*r*∇*v* part of the flux **J** is directed towards the decreasing population density of mid level predator *v*. This is because, in certain prey-predator systems, several prey form a huge group to protect themselves from the attacking predators. In addition, as in the model system the predation of top predator is impossible, to extensively study the role of cross diffusion on the model system, we induce it at *r*, i.e., the top predator.

#### **3. Linear Stability Analysis of Temporal Model**

Under the assumption that the model system in Equation in Equation (6) has the unique positive stationary uniform solution, which we denote by **u**∗ = (*u*∗, *v*∗,*r*∗) (where **u**∗ is a three coordinate vector), we derive the conditions under which it is locally asymptotically stable.

**Theorem 1.** *The positive equilibrium solution given by* **u**∗ = (*u*∗, *v*∗,*r*∗) *of the model system in Equation (6) is locally asymptotically stable if the parameters satisfy,*

$$\frac{r^\* v^\* \left(1 + w \wp r^\*\right)}{\left[v^\* + \left(w\_8 + w\_9 v^\*\right) r^\* + w\_{10}\right]^2} < e \quad \text{and} \quad \frac{v^\*}{\left(u^\* + w\_4\right)^2} < 1. \tag{9}$$

**Proof.** We consider the following notation given in Equation (10):

$$\mathbf{G}(\mathbf{u}) = \begin{pmatrix} G\_1(\mathbf{u}) \\ G\_2(\mathbf{u}) \\ G\_3(\mathbf{u}) \end{pmatrix} = \begin{pmatrix} u\mathbf{g}\_1(\mathbf{u}) & \cong \, u \left(1 - u - \frac{v}{u + w\_4}\right) \\ v\mathbf{g}\_2(\mathbf{u}) & \cong \, v \left(-w\_5 + \frac{w\_6 u}{u + w\_7} - \varepsilon v - \frac{r}{v + \left(w\_8 + w\_9 v\right)r + w\_{10}}\right) \\ r\mathbf{g}\mathbf{s}(\mathbf{u}) & \cong \, r \left(-w\_{11} + \frac{w\_{12} v}{v + \left(w\_8 + w\_9 v\right)r + w\_{10}}\right) \end{pmatrix} \tag{10}$$

Calculating **Gu**(**u**∗), and putting, **Gu**(**u**∗) = 0.

$$\mathbf{G}\_{\mathbf{u}}(\mathbf{u}^{\*}) = \begin{pmatrix} \frac{u^{\*}v^{\*}}{(u^{\*} + w\_{4})^{2}} - u^{\*} & \frac{-u^{\*}}{u^{\*} + w\_{4}} & 0\\ \frac{w\_{6}w\_{5}v^{\*}}{(u^{\*} + w\_{7})^{2}} & -\varepsilon v^{\*} + \frac{r^{\*}v^{\*}(1 + w\_{9}r^{\*})}{[v^{\*} + (w\_{8} + w\_{9}r^{\*})r^{\*} + w\_{10}]^{2}} & -\frac{v^{\*}(v^{\*} + w\_{10})}{[v^{\*} + (w\_{8} + w\_{9}r^{\*})r^{\*} + w\_{10}]^{2}}\\ 0 & \frac{r^{\*}w\_{12}(w\_{9}r^{\*} + w\_{10})}{[v^{\*} + (w\_{8} + w\_{9}r^{\*})r^{\*} + w\_{10}]^{2}} & \frac{-w\_{12}v^{\*}r^{\*}(w\_{8} + w\_{9}r^{\*})}{[v^{\*} + (w\_{8} + w\_{9}r^{\*})r^{\*} + w\_{10}]^{2}} \end{pmatrix} \tag{11}$$

Here, *ρ*(*λ*) denotes the characteristic polynomial of **Gu**(**u**∗) given by:

$$
\rho(\lambda) = \lambda^3 + H\_1 \lambda^2 + H\_2 \lambda + H\_3,
$$

where

$$H\_1 = \frac{w\_{12}v^\*r^\*\left(ws + w\wp v^\*\right)}{\left[v^\* + \left(w\_8 + w\wp v^\*\right)r^\* + w\_{10}\right]^2} + v^\*\left(\varepsilon - \frac{r^\*v^\*(1+w\wp r^\*)}{\left[v^\* + \left(w\_8 + w\wp v^\*\right)r^\* + w\_{10}\right]^2}\right) + u^\*\left(1 - \frac{v^\*}{\left(u^\* + w\_4\right)^2}\right)$$

*H*<sup>2</sup> =*u*<sup>∗</sup> <sup>1</sup> <sup>−</sup> *<sup>v</sup>*<sup>∗</sup> (*u* + *w*4)<sup>2</sup> *w*12*v*∗*r*∗(*w*<sup>8</sup> + *w*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] <sup>2</sup> + *v*<sup>∗</sup> - *<sup>e</sup>* <sup>−</sup> *<sup>r</sup>*∗*v*∗(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 & <sup>+</sup> *<sup>w</sup>*6*w*7*u*∗*v*<sup>∗</sup> (*u*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*4)(*u*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*7)<sup>2</sup> <sup>+</sup> *<sup>v</sup>*∗*r*∗*w*12(*v*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*10)(*w*8*r*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*10) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] <sup>2</sup> <sup>+</sup> *<sup>w</sup>*12*v*∗2*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 × - *<sup>e</sup>* <sup>−</sup> *<sup>r</sup>*∗*v*∗(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 *<sup>H</sup>*<sup>3</sup> <sup>=</sup> *<sup>w</sup>*6*w*7*u*∗*v*∗2*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) (*u* + *w*4)(*u* + *w*7)<sup>2</sup> [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] <sup>2</sup> + *v*∗*r*∗*w*12(*v*<sup>∗</sup> + *w*10)(*w*8*r*<sup>∗</sup> + *w*10) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 + - *w*12*v*∗2*r*∗(*w*<sup>8</sup> + *w*9) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 *<sup>e</sup>* <sup>−</sup> *<sup>r</sup>*∗*v*∗(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 & × *u*∗ <sup>1</sup> <sup>−</sup> *<sup>v</sup>*<sup>∗</sup> (*u*<sup>∗</sup> + *w*4)<sup>2</sup> &

By using the criteria stated in Theorem 1, it is easy to verify that *H*1, *H*2, *H*<sup>3</sup> > 0 and *H*1*H*<sup>2</sup> − *H*<sup>3</sup> > 0. Therefore, from the Routh-Hurwitz criteria that the three roots *ρ*(*λ*) = 0 have negative real parts, the positive equilibrium solution **u**∗ of the model system in Equation (6) is locally asymptotically stable under the condition in Equation (9).

#### **4. Analysis of Spatially Extended Model**

In this section, we focus on the spatiotemporal models in Equation (7) in presence of diffusion but absence of cross diffusion and Equation (8) in presence of both self and cross diffusion. Our objective is to derive the conditions under which the positive equilibrium solution is stable with diffusion but without cross diffusion but unstable with diffusion and cross diffusion both. This phenomenon is called *cross diffusion driven instability* [16,18].

#### *4.1. Without Cross Diffusion*

We now consider the system in Equation (7). To discuss the local stability analysis of the system of parabolic Equations (7) and (8), we use the notation used in [16].

**Notation 1.** *Let* 0 = *μ*<sup>1</sup> < *μ*<sup>2</sup> < ... *be the eigenvalues of* −Δ *on* Ω *under no-flux boundary conditions, and E* (*μi*) *be the space of eigenfunctions corresponding to μi. We define the following space decomposition*

$$\text{A } (i) \quad \mathbf{X}\_{\overline{\mathbf{i}}} \coloneqq \left\{ \mathbf{c} \, \, \middle| \, \, \phi\_{\overline{\mathbf{i}}} \colon \mathbf{c} \in \mathbf{R}^3 \right\} \text{ where } \left\{ \phi\_{\overline{\mathbf{i}}} \right\} \text{ are orthonormal basis of } \mathbb{E}\left(\mu\_{\overline{\mathbf{i}}}\right) \text{ for } \mathbf{j} = 1, \dots, dim \text{E}\left(\mu\_{\overline{\mathbf{i}}}\right).$$

$$\begin{array}{rcl} \text{(iii)} & \text{\textbf{X}} & := & \left\{ \mathbf{u} \in \left[ \mathbb{C}^{1}(\Omega) \right]^{\circ} : \partial\_{\boldsymbol{n}} \boldsymbol{n} = \stackrel{\scriptstyle \mathbf{\mathcal{I}}^{\circ} \boldsymbol{\mathcal{I}}^{\circ} \boldsymbol{\mathcal{I}} = \boldsymbol{0} \text{ or } \partial \Omega \right\} \text{, and so } \mathsf{X} & = & \bigoplus\_{l=1}^{\infty} \mathsf{X}\_{\mathbf{i}l} \text{ where } \mathsf{X}\_{\mathbf{i}} = \bigoplus\_{l=1}^{\infty} \mathsf{X}\_{\mathbf{i}l}.\\ & & \bigoplus\_{j=1}^{\dim \mathsf{E}(\mu\_{i})} \mathsf{X}\_{\mathbf{i}j}. \end{array}$$

**Theorem 2.** *The positive equilibrium solution* **u**∗ *of the model system in Equation (7) is locally asymptotically stable if*

$$\frac{r^\* v^\* \left(1 + w\_{\theta} r^\*\right)}{\left[v^\* + \left(w\_{\theta} + w\_{\theta} v^\*\right) r^\* + w\_{10}\right]^2} < e \quad \text{and} \quad \frac{v^\*}{\left(u^\* + w\_4\right)^2} < 1. \tag{12}$$

**Proof.** The spatially extended model system without cross diffusion in Equation (7) has been linearized at **u**∗ and expressed as [16]:

$$\mathbf{u}\_{\mathbf{t}} = \left(D\boldsymbol{\Delta} + \mathbf{G}\_{\mathbf{u}}(\mathbf{u}^\*)\right)\mathbf{u},\tag{13}$$

where,

$$D = \begin{bmatrix} d\_1 & 0 & 0 \\ 0 & d\_2 & 0 \\ 0 & 0 & d\_3 \end{bmatrix}$$

Notation 1 suggests that **Xi** is invariant under the operator *D*Δ + **Gu**(**u**∗), and *λ* is an eigenvalue of this operator on **Xi** if and only if it is an eigenvalue of the matrix in *D*Δ + **Gu**(**u**∗). The characteristic polynomial of the matrix *μiD* + **Gu**(**u**∗) is given by

$$
\psi\_{\lambda} = \lambda^3 + A\_1 \lambda^2 + A\_2 \lambda + A\_3,
$$

where

*<sup>A</sup>*<sup>1</sup> <sup>=</sup>*μ<sup>i</sup>* (*d*<sup>1</sup> <sup>+</sup> *<sup>d</sup>*<sup>2</sup> <sup>+</sup> *<sup>d</sup>*3) <sup>+</sup> *<sup>w</sup>*12*v*∗*r*<sup>∗</sup> (*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] <sup>2</sup> + *v*<sup>∗</sup> - *<sup>e</sup>* <sup>−</sup> *<sup>r</sup>*∗*v*∗(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 + *u*∗ <sup>1</sup> <sup>−</sup> *<sup>v</sup>*<sup>∗</sup> (*u*<sup>∗</sup> + *w*4)<sup>2</sup> *A*<sup>2</sup> =*μid*<sup>1</sup> *<sup>μ</sup>i*(*d*<sup>2</sup> <sup>+</sup> *<sup>d</sup>*3) + *<sup>w</sup>*12*v*∗*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] <sup>2</sup> + *v*<sup>∗</sup> - *<sup>e</sup>* <sup>−</sup> *<sup>r</sup>*∗*v*∗(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 & + *u*∗ <sup>1</sup> <sup>−</sup> *<sup>v</sup>*<sup>∗</sup> (*u*<sup>∗</sup> + *w*4)<sup>2</sup> *<sup>μ</sup>i*(*d*<sup>2</sup> <sup>+</sup> *<sup>d</sup>*3) + *<sup>w</sup>*12*v*∗*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 + *v*∗ - *<sup>e</sup>* <sup>−</sup> *<sup>r</sup>*∗*v*∗(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 & <sup>+</sup> *<sup>w</sup>*6*w*7*u*∗*v*<sup>∗</sup> (*u*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*4)(*u*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*7)<sup>2</sup> <sup>+</sup> *<sup>μ</sup>*<sup>2</sup> *<sup>i</sup> d*2*d*<sup>3</sup> <sup>+</sup> *<sup>μ</sup>id*2*w*12*v*∗*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] <sup>2</sup> <sup>+</sup> *<sup>v</sup>*∗*r*∗*w*12(*v*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*10)(*w*8*r*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*10) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 + - *<sup>μ</sup>id*3*v*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*12*v*∗2*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 *<sup>e</sup>* <sup>−</sup> *<sup>r</sup>*∗*v*∗(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 *A*<sup>3</sup> =*μid*<sup>3</sup> *w*6*w*7*u*∗*v*∗ (*u*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*4)(*u*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*7)<sup>2</sup> <sup>+</sup> *<sup>w</sup>*6*w*7*u*∗*v*∗2*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) (*u*<sup>∗</sup> + *w*4)(*u* + *w*7)<sup>2</sup> [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 + *μ*2 *<sup>i</sup> <sup>d</sup>*2*d*<sup>3</sup> <sup>+</sup> *<sup>μ</sup>id*2*w*12*v*∗*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] <sup>2</sup> <sup>+</sup> *<sup>v</sup>*∗*r*∗*w*12(*v*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*10)(*w*8*r*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*10) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 + - *<sup>μ</sup>id*3*v*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*12*v*∗2*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 *<sup>e</sup>* <sup>−</sup> *<sup>r</sup>*∗*v*∗(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 & × *μid*<sup>3</sup> + *u*<sup>∗</sup> <sup>1</sup> <sup>−</sup> *<sup>v</sup>*<sup>∗</sup> (*u*<sup>∗</sup> + *w*4)<sup>2</sup> &

By using the stated criteria in Theorem 2, it seems clear that *A*1, *A*2, *A*<sup>3</sup> > 0 and *A*1*A*<sup>2</sup> − *A*<sup>3</sup> > 0. Therefore, it follows from the Routh-Hurwitz criteria that, for each *i* ≥ 1, all three roots *λi*,1, *λi*,2, *λi*,3 of *ψi*(*λ*) = 0 have negative real parts. Hence, the positive equilibrium solution **u**<sup>∗</sup> of Equation (7) is locally asymptotically stable under the stated condition in Equation (12).

From Theorem 2, it is clear that addition of self diffusion terms to the temporal model system in Equation (6) results in stable positive equilibrium solution under the stated condition. Therefore, diffusion driven instability has not yet occurred. Now, we now consider the effect of cross diffusion on the model system.

#### *4.2. With Cross Diffusion*

In this section, we study in detail the system given in Equation (8).

**Theorem 3.** *Provided that Theorem 2 holds, assuming d*<sup>4</sup> > 0 *and considering the following inequality holds,*

$$\begin{cases} \varepsilon - \frac{r^\*(1+w\_9r^\*)}{\left[v^\* + (w\_8+w\_9r^\*)\,r^\* + w\_{10}\right]^2} \Bigg\{ \left\{1 - \frac{v^\*}{\left(u^\* + w\_4\right)^2} \right\} + \frac{w\_6w\_7}{\left(u^\* + w\_4\right)\left(u^\* + w\_7\right)^2} \\ \quad < \frac{d\_4r^\*(r^\* + w\_{10})}{\left(1 + d\_4v^\*\right)\left[v^\* + \left(w\_8 + w\_9r^\*\right)\,r^\* + w\_{10}\right]^2} \Bigg\} \\ \end{cases} \tag{14}$$

*if <sup>μ</sup>*<sup>2</sup> <sup>&</sup>lt; *<sup>μ</sup>*˜*, where <sup>μ</sup>*<sup>2</sup> *is as explained in Notation 1 and <sup>μ</sup> is given by Equation (24), then there exists <sup>d</sup>*<sup>∗</sup> <sup>3</sup> > 0*, such that, when d*<sup>3</sup> ≥ *d*<sup>∗</sup> <sup>3</sup>*, the positive equilibrium solution* **u**<sup>∗</sup> *of the cross diffusive system in Equation (8) becomes unstable.*

We denote Φ(**u**)=(*d*1*u*, *d*2*v*, *d*3*r*(1 + *d*4*v*))*T*. Linearizing the system in Equation (8) at **u**∗, we have

$$\mathbf{u}\_{\mathbf{t}} = \left(\Phi\_{\boldsymbol{\upmu}} \boldsymbol{\upLambda} + \mathbf{G}\_{\mathbf{u}}(\mathbf{u}^\*)\right) \mathbf{u},\tag{15}$$

where,

$$
\Phi\_{\mathbb{M}} = \begin{bmatrix} d\_1 & 0 & 0 \\ 0 & d\_2 & 0 \\ 0 & d\_3 d\_4 r^\* & d\_3 + d\_3 d\_4 v^\* \end{bmatrix}
$$

The characteristic polynomial of −*μi*Φ*<sup>u</sup>* + **Gu**(**u**∗) is given by,

$$
\psi\_{\lambda} = \lambda^3 + B\_1 \lambda^2 + B\_2 \lambda + B\_3 \tag{16}
$$

where

$$\begin{aligned} B\_1 &= \mu\_i (d\_1 + d\_2 + d\_3) + \mu\_i d\_3 d\_4 v^\* + \frac{w\_{12} v^\* r^\* (w\_8 + w\_9 v^\*)}{\left[v^\* + (w\_8 + w\_9 v^\*) \, r^\* + w\_{10}\right]^2} \\ &+ v^\* \left(\varepsilon - \frac{r^\* v^\* (1 + w\_9 r^\*)}{\left[v^\* + (w\_8 + w\_9 v^\*) \, r^\* + w\_{10}\right]^2}\right) + u^\* \left(1 - \frac{v^\*}{(u^\* + w\_4)^2}\right), \end{aligned}$$

*B*<sup>2</sup> = *<sup>μ</sup>i*(*d*<sup>2</sup> <sup>+</sup> *<sup>d</sup>*3) + *<sup>μ</sup>id*3*d*4*v*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*12*v*∗*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 + *v*∗ - *<sup>e</sup>* <sup>−</sup> *<sup>r</sup>*∗*v*∗(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 & × 0 *μid*<sup>1</sup> + *u*<sup>∗</sup> <sup>1</sup> <sup>−</sup> *<sup>v</sup>*<sup>∗</sup> (*u*<sup>∗</sup> + *w*4)<sup>2</sup> 1 + *μ*<sup>2</sup> *<sup>i</sup> <sup>d</sup>*2*d*3(<sup>1</sup> <sup>+</sup> *<sup>d</sup>*4*v*∗) + *<sup>μ</sup>id*3*ev*∗(<sup>1</sup> <sup>+</sup> *<sup>d</sup>*4*v*∗) + *<sup>w</sup>*12*ev*∗2*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 <sup>+</sup> *<sup>μ</sup>id*2*w*12*v*∗*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] <sup>2</sup> <sup>+</sup> *<sup>w</sup>*12*v*∗*r*∗(*v*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*10)(*w*8*r*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*10) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*r*∗) *v*<sup>∗</sup> + *w*10] 4 <sup>−</sup> *<sup>μ</sup>id*3*r*∗*v*∗(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] <sup>2</sup> (<sup>1</sup> <sup>+</sup> *<sup>d</sup>*4*v*∗) <sup>−</sup> *<sup>μ</sup>id*3*d*4*r*∗*v*∗(*v*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*10) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 <sup>−</sup> *<sup>w</sup>*12*v*∗2*r*∗2(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗)(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*r*∗) *v*<sup>∗</sup> + *w*10] 4 *<sup>B</sup>*<sup>3</sup> <sup>=</sup> *<sup>μ</sup>id*3*w*6*w*7*u*∗*v*<sup>∗</sup> (*u*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*4)(*u*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*7)<sup>2</sup> (<sup>1</sup> <sup>+</sup> *<sup>d</sup>*4*v*∗) + *<sup>w</sup>*6*w*7*u*∗*v*∗2*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*∗) (*u*<sup>∗</sup> + *w*4)(*u*<sup>∗</sup> + *w*7)<sup>2</sup> [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 0 *μid*<sup>1</sup> + *u*<sup>∗</sup> <sup>1</sup> <sup>−</sup> *<sup>v</sup>*<sup>∗</sup> (*u*<sup>∗</sup> + *w*4)<sup>2</sup> 1 × *μ*2 *<sup>i</sup> d*2*d*3(1 + *d*4*v*∗) + *μid*3*ev*∗(1 + *d*4*v*∗) <sup>+</sup> *<sup>w</sup>*12*ev*∗2*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] <sup>2</sup> <sup>+</sup> *<sup>μ</sup>id*2*w*12*v*∗*r*∗(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 <sup>+</sup> *<sup>w</sup>*12*v*∗*r*∗(*v*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*10)(*w*8*r*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*10) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*r*∗) *v*<sup>∗</sup> + *w*10] <sup>4</sup> <sup>−</sup> *<sup>μ</sup>id*3*d*4*r*∗*v*∗(*v*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*10) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 <sup>−</sup> *<sup>μ</sup>id*3*r*∗*v*∗(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] <sup>2</sup> (<sup>1</sup> <sup>+</sup> *<sup>d</sup>*4*v*∗) <sup>−</sup> *<sup>w</sup>*12*v*∗2*r*∗2(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗)(*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*r*∗) *v*<sup>∗</sup> + *w*10] 4 &

If *λ*1(*μi*), *λ*2(*μi*), *λ*3(*μi*) are the three roots of *ψ*(*λ*) = 0, then

$$
\lambda\_1(\mu\_i) . \lambda\_2(\mu\_i) . \lambda\_3(\mu\_i) = -B\_3. \tag{17}
$$

For at least one **Re**(*λi*(*μi*)) > 0, it suffices to show that *B*<sup>3</sup> < 0. In addition,

$$B\mathfrak{z} = \det(\mu\_i \phi\_i - \mathbf{G}\_\mathbf{u}(\mathbf{u}^\*)) \tag{18}$$

Hence, we have,

$$B\_3 = b\_3 \mu\_i^3 + b\_2 \mu\_i^2 + b\_1 \mu\_i + b\_0 \tag{19}$$

where

*b*<sup>3</sup> = *d*1*d*2*d*3(1 + *d*3*v*∗)

$$\begin{split} b\_{2} &= d\_{1}d\_{3} (1 + d\_{4}\upsilon^{\*})\upsilon^{\*} \left(\upsilon - \frac{r^{\*}\upsilon^{\*}(1 + \imath\eta\upsilon^{\*})}{\left[\upsilon^{\*} + \left(\upsilon\kappa\_{8} + \imath\eta\upsilon^{\*}\right)r^{\*} + \upsilon\iota\_{10}\right]^{2}}\right) + \frac{d\_{1}d\_{2}\upsilon\upsilon\_{12}(\upsilon\kappa\_{8} + \imath\eta\upsilon^{\*})\upsilon^{\*}r^{\*}}{\left[\upsilon^{\*} + \left(\upsilon\kappa\_{8} + \imath\eta\upsilon^{\*}\right)r^{\*} + \upsilon\iota\_{10}\right]^{2}} \\ &- \frac{d\_{1}d\_{3}d\_{4}\upsilon^{\*}\upsilon^{\*}(\upsilon^{\*} + \upsilon\iota\_{10})}{\left[\upsilon^{\*} + \left(\upsilon\kappa\_{8} + \imath\eta\upsilon^{\*}\right)r^{\*} + \upsilon\iota\_{10}\right]^{2}} + d\_{2}d\_{3}(1 + d\_{4}\upsilon^{\*})\left(\upsilon^{\*} - \frac{\imath\upsilon^{\*}\upsilon^{\*}}{\left(\upsilon + \upsilon\iota\_{4}\right)^{2}}\right) \end{split}$$

*b*<sup>1</sup> =*d*1*v*<sup>∗</sup> - *<sup>e</sup>* <sup>−</sup> *<sup>r</sup>*∗*v*∗(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 *w* − 12*v*∗*r*∗(*w*<sup>8</sup> + *w*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 *d*1*w*12*v*∗*r*∗(*v*<sup>∗</sup> + *w*10)(*w*8*r*<sup>∗</sup> + *w*10) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*r*∗) *v*<sup>∗</sup> + *w*10] <sup>4</sup> + *d*2*w*12*v*∗*r*∗(*w*<sup>8</sup> + *w*9*v*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 *<sup>u</sup>*<sup>∗</sup> <sup>−</sup> *<sup>u</sup>*∗*v*<sup>∗</sup> (*u* + *w*4)<sup>2</sup> *d*3(1 + *d*4*v*∗) - *<sup>e</sup>* <sup>−</sup> *<sup>r</sup>*∗*v*∗(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*9*r*∗) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 *<sup>u</sup>*<sup>∗</sup> <sup>−</sup> *<sup>u</sup>*∗*v*<sup>∗</sup> (*u* + *w*4)<sup>2</sup> <sup>−</sup> *<sup>d</sup>*3*d*4*r*∗*v*∗(*v*<sup>∗</sup> <sup>+</sup> *<sup>w</sup>*10) [*v*<sup>∗</sup> + (*w*<sup>8</sup> + *w*9*v*∗)*r*<sup>∗</sup> + *w*10] 2 *<sup>u</sup>*<sup>∗</sup> <sup>−</sup> *<sup>u</sup>*∗*v*<sup>∗</sup> (*u* + *w*4)<sup>2</sup> <sup>+</sup> *<sup>d</sup>*3(<sup>1</sup> <sup>+</sup> *<sup>d</sup>*4*v*∗) *<sup>w</sup>*6*w*7*u*∗*v*<sup>∗</sup> (*u* + *w*4)(*u* + *w*7)<sup>2</sup>

$$\begin{split} b\_{0} &= \left\{ u^{\*} - \frac{u^{\*}v^{\*}}{(u^{\*} + w\_{4})^{2}} \right\} \left( e - \frac{r^{\*}v^{\*}(1 + w\_{9}r^{\*})}{\left[v^{\*} + (w\_{8} + w\_{9}r^{\*})\right.r^{\*} + w\_{10}]^{2}} \right) \frac{w\_{12}v^{\*}r^{\*}(w\_{8} + w\_{9}v^{\*})}{\left[v^{\*} + (w\_{8} + w\_{9}r^{\*})\right.r^{\*} + w\_{10}]^{2}} \\ &+ \left\{ u^{\*} - \frac{u^{\*}v^{\*}}{(u^{\*} + w\_{4})^{2}} \right\} \frac{v^{\*}r^{\*}w\_{12}(v^{\*} + w\_{10})(w\_{8}r^{\*} + w\_{10})}{\left[v^{\*} + (w\_{8} + w\_{9}r^{\*})\right.r^{\*} + w\_{10}]^{4}} \\ &+ \frac{w\_{6}w\_{7}w\_{12}u^{\*}v^{\*2}r^{\*}(w\_{8} + w\_{9}v^{\*})}{\left(u + w\_{4}\right)(u + w\_{7})^{2}\left[v^{\*} + (w\_{8} + w\_{9}r^{\*})\right.r^{\*} + w\_{10}\right]^{2}} \end{split}$$

From the above expression, we clearly see that *b*<sup>0</sup> = −det(**Gu**(**u**∗)).

Let, ˜ *<sup>b</sup>*(*μ*) = *<sup>b</sup>*3*μ*<sup>3</sup> <sup>+</sup> *<sup>b</sup>*2*μ*<sup>2</sup> <sup>+</sup> *<sup>b</sup>*1*<sup>μ</sup>* <sup>−</sup> det(**Gu**(**u**∗) and let *<sup>μ</sup>*˜1, *<sup>μ</sup>*˜2, *<sup>μ</sup>*˜3 be the three roots of ˜ *<sup>b</sup>*(*μ*) = 0, such that **Re**(*μ*˜1) ≤ **Re**(*μ*˜2) ≤ **Re**(*μ*˜3). Now, *μ*˜1 *μ*˜2 *μ*˜3 = det(**Gu**(**u**∗)). As det(**Gu**(**u**∗)) < 0 due to condition specified in Equation (12), we have,

$$
\bar{\mu}\_1 \bar{\mu}\_2 \tilde{\mu}\_3 < 0
$$

As *<sup>d</sup>*<sup>1</sup> <sup>&</sup>gt; 0, *<sup>d</sup>*<sup>2</sup> <sup>&</sup>gt; 0, *<sup>d</sup>*<sup>3</sup> <sup>&</sup>gt; 0 and *<sup>d</sup>*<sup>4</sup> <sup>&</sup>gt; 0, *<sup>b</sup>*<sup>3</sup> <sup>&</sup>gt; 0. Using the theory of equations if one of *<sup>μ</sup>*21, *<sup>μ</sup>*22, *<sup>μ</sup>*2<sup>3</sup> is real and negative and the product of the other two is positive.

To obtain the values of *d*∗ <sup>3</sup>, we calculate the following limits:

$$\lim\_{d\_3 \to \infty} \frac{b\_3}{d\_3} = d\_1 d\_2 (1 + d\_4 v^\*) \stackrel{\sim}{=} \rho\_3 \tag{20}$$

(*u* + *w*4)(*u* + *w*7)<sup>2</sup>

∼= *ρ*<sup>1</sup>

$$\begin{split} \lim\_{d \to \infty} \frac{b\_2}{d\_3} &= d\_1 (1 + d\_4 v^\*) v^\* \left( c - \frac{r^\* v^\* (1 + wy^\* \, ^\*)}{\left[ v^\* + (wy + wy \, ^\*) \, r^\* + w\_{10} \right]^2} \right) - \frac{d\_1 d\_4 v^\* v^\* (v^\* + w\_{10})}{\left[ v^\* + (wy + wy \, v^\*) \, r^\* + w\_{10} \right]^2} \\ &+ d\_2 (1 + d\_4 v^\*) \left( u^\* - \frac{u^\* v^\*}{(u + w\_4)^2} \right) \cong \rho\_2 \\ \lim\_{d \to \infty} \frac{b\_1}{d\_3} &= (1 + d\_4 v^\*) v^\* \left( c - \frac{r^\* v^\* (1 + u\_3 v^\*)}{\left[ v^\* + (wy + wy \, v^\*) \, r^\* + w\_{10} \right]^2} \right) \left( u^\* - \frac{u^\* v^\*}{(u + w\_4)^2} \right) \\ &- \frac{d\_4 v^\* v^\* (v^\* + w\_{10})}{\left[ v^\* + (w\_{10} + w\_{20}) v^\* \, v^\* - w\_{10} \right]^2} \left( u^\* - \frac{u^\* v^\*}{(u + w\_4)^2} \right) + (1 + d\_4 v^\*) \left( \frac{w\_6 w \eta \, u^\* v^\*}{(u + w\_6)(u + w\_7)^2} \right) \cong \rho\_1. \end{split} \tag{22}$$

$$-\frac{\omega \bullet \left[\left(\nu^\* + (\nu \wp + \nu \wp^\*)\right)r^\* + \nu \nu\_{10}\right]^2}{\left[\nu^\* + (\nu \wp + \nu \wp^\*)r^\* + \nu \nu\_{10}\right]^2} \left(\mu^\* - \frac{\omega \bullet \nu}{(\mu + \nu \nu\_4)^2}\right) + (1 + d\_4 \nu^\*)^2$$
 
$$\text{If } \mu^\* < 0, \text{ then we then find } \nu \text{ in a non-null time ball } \mu.$$

If *ρ*<sup>1</sup> < 0, then the following inequality holds,

$$\begin{aligned} & \left\{ \varepsilon - \frac{r^\* \left( 1 + w\_{\mathfrak{z}} r^\* \right)}{\left[ v^\* + \left( w\_{\mathfrak{z}} + w\_{\mathfrak{z}} v^\* \right) r^\* + w\_{\mathfrak{z}} \right]^2} \right\} \left\{ 1 - \frac{v^\*}{\left( u^\* + w\_4 \right)^2} \right\} + \frac{w\_6 w\_{\mathfrak{z}}}{\left( u^\* + w\_4 \right) \left( u^\* + w\_{\mathfrak{z}} \right)^2} \\ & < \frac{d\_4 r^\* \left( r^\* + w\_{\mathfrak{z}} \right)}{\left( 1 + d\_4 v^\* \right) \left[ v^\* + \left( w\_8 + w\_9 v^\* \right) r^\* + w\_{\mathfrak{z}} \right]^2} \left( 1 - \frac{v^\*}{\left( u^\* + w\_4 \right)^2} \right) \end{aligned}$$

We have, *ρ*<sup>1</sup> < 0 < *ρ*3.

In addition,

$$\lim\_{d\_3 \to \infty} \frac{\tilde{b}(\mu)}{d\_3} = \rho\_3 \mu^3 + \rho\_2 \mu^2 + \rho\_1 \mu = \mu(\rho\_3 \mu^2 + \rho\_2 \mu + \rho\_1)$$

As *ρ*<sup>1</sup> < 0 < *ρ*3, which means that the equation,

$$\lim\_{d\_3 \to \infty} \frac{b(\mu)}{d\_3} = 0$$

has one positive root and one negative root. From continuity, if *<sup>d</sup>*<sup>3</sup> <sup>→</sup> <sup>∞</sup>, *<sup>μ</sup>* is real and negative. Since *<sup>μ</sup>*22*μ*2<sup>3</sup> <sup>&</sup>gt; 0, *<sup>μ</sup>*2<sup>2</sup> and *<sup>μ</sup>*2<sup>3</sup> are real and positive.

$$\lim\_{d\_2 \to \infty} \bar{\mu}\_1 = \frac{-\rho\_2 - \sqrt{\rho\_2^2 - 4\rho\_1 \rho\_3}}{2\rho\_3} < 0 \tag{23}$$

$$\lim\_{d\_2 \to \infty} \bar{\mu}\_3 = \frac{-\rho\_2 + \sqrt{\rho\_2^2 - 4\rho\_1 \rho\_3}}{2\rho\_3} > 0 = \hat{\mu} \tag{24}$$

$$\lim\_{d\_2 \to \infty} \bar{p}\_2 = 0 \tag{25}$$

Therefore, there exists *d*∗ <sup>3</sup> > 0 such that when *d*<sup>3</sup> ≥ *d*<sup>∗</sup> <sup>3</sup> then the specified criteria holds. In addition, ˜ *b*(*μ*) < 0 when *μ* ∈ (−∞, *μ*˜1) ∪ (*μ*˜2, *μ*˜3). Therefore, when 0 < *μ*<sup>2</sup> < *μ*˜, then *μ*<sup>2</sup> ∈ (*μ*˜2, *μ*˜3), it follows that ˜ *b*(*μ*2) < 0. Therefore, *B*<sup>3</sup> < 0, and the proof is complete.

Therefore, from the above theorems, we conclude that cross diffusion destabilizes the stationary uniform solution.

#### **5. Inhomogeneous Steady States**

In this section, the justification for the cross diffusion driven instability phenomenon is explained. We now prove that model system in Equation (8) admits the inhomogeneous steady state. We, consider the steady state model system in Equation (8) in the following form,

$$\begin{aligned} &1 - d\_1 \Delta u = u \left(1 - u - \frac{v}{u + w\_4} \right) \\ &1 - d\_2 \Delta v = v \left(-w\_5 + \frac{w\_6 u}{u + w\_7} - \varepsilon v - \frac{r}{v + \left(w\_8 + w\_9 v \right) r + w\_{10}} \right) \\ &1 - d\_3 \Delta \left(r + d\_4 v r \right) = r \left(-w\_{11} + \frac{w\_{12} v}{v + \left(w\_8 + w\_9 v \right) r + w\_{10}} \right) \\ &\frac{\partial u}{\partial n} = \frac{\partial v}{\partial n} = \frac{\partial r}{\partial n} = 0 \qquad \text{x, y} \in \Omega. \end{aligned} \tag{26}$$

The constants *Q*∗, *Q*, *Q*¯ depend on the domain Ω. As Ω is fixed, this dependance is not mentioned explicitly. The parameters *w*4, *w*5, *w*6, *w*7,*e*, *w*8, *w*9, *w*10, *w*11, *w*<sup>12</sup> are collectively denoted by Γ.

#### *5.1. A Priori Estimates*

We now attempt to give an a priori upper and lower bounds for the positive solutions to Equation (26). We use Harnack Inequality and Maximum Principle (for details, refer to [16,18]).

**Proposition 1.** *Harnack Inequality Let <sup>f</sup>* <sup>∈</sup> *<sup>C</sup>*2(Ω) <sup>∩</sup> *<sup>C</sup>*1(Ω¯ ) *be a positive solution to* <sup>Δ</sup>*w*(*x*) + *<sup>c</sup>*(*x*)*w*(*x*) = <sup>0</sup>*, where <sup>c</sup>* <sup>∈</sup> *<sup>C</sup>*(Ω¯ )*, satisfying the homogeneous Neumann boundary condition. Then, there exists a positive constants C*<sup>∗</sup> *that depends only on* ||*c*||<sup>∞</sup> *such that*

$$\max\_{\Omega} w \le \mathcal{C}\_\* \min\_{\Omega} w$$

**Proposition 2.** *Maximum principle Let g* <sup>∈</sup> *<sup>C</sup>*(<sup>Ω</sup> <sup>×</sup> <sup>R</sup>1) *and bij* <sup>∈</sup> *<sup>C</sup>*(Ω¯ ), *<sup>j</sup>* <sup>=</sup> 1, 2, . . . , *<sup>N</sup>*

*1. If f* <sup>∈</sup> *<sup>C</sup>*2(Ω) <sup>∩</sup> *<sup>C</sup>*1(Ω¯ ) *satisfies*

$$\begin{cases} \quad \Delta f(\mathbf{x}) + \sum\_{j=1}^{N} b\_j(\mathbf{x}) w\_{\mathbf{x}j} + \mathfrak{g}(\mathbf{x}, w(\mathbf{x})) \ge 0, & \mathbf{x} \in \Omega, \\\ \partial\_{\mathbf{n}} w(\mathbf{x}) \le 0 & \mathbf{x} \in \partial\Omega. \end{cases}$$

*and f*(*x*0) = maxΩ¯ *f , then g*(*x*0, *f*(*x*0)) ≥ 0*. 2. If f* <sup>∈</sup> *<sup>C</sup>*2(Ω) <sup>∩</sup> *<sup>C</sup>*1(Ω¯ ) *satisfies*

$$\begin{cases} \quad \Delta f(\mathbf{x}) + \sum\_{j=1}^{N} b\_j(\mathbf{x}) f\_{\mathbf{x}j} + \mathbf{g}(\mathbf{x}, f(\mathbf{x})) \le 0, & \mathbf{x} \in \Omega, \\\ \quad \partial\_n f(\mathbf{x}) \ge 0 & \mathbf{x} \in \partial \Omega. \end{cases}$$

*and f*(*x*0) = minΩ¯ *f , then g*(*x*0, *f*(*x*0)) ≤ 0*.*

**Theorem 4.** *(Upper Bound) Let δ*1, *δ*2, *δ*<sup>3</sup> *and δ*<sup>4</sup> *be fixed positive constants. Then, there exists positive constants <sup>Q</sup>*∗(Γ, *<sup>δ</sup>*1, *<sup>δ</sup>*2, *<sup>δ</sup>*3, *<sup>δ</sup>*4) *and <sup>Q</sup>*¯(Γ, *<sup>δ</sup>*1, *<sup>δ</sup>*2, *<sup>δ</sup>*3, *<sup>δ</sup>*4) *such that when <sup>d</sup>*<sup>1</sup> <sup>≥</sup> *<sup>δ</sup>*1*; <sup>d</sup>*<sup>2</sup> <sup>≥</sup> *<sup>δ</sup>*2*; <sup>d</sup>*<sup>3</sup> <sup>≥</sup> *<sup>δ</sup>*<sup>3</sup> *and <sup>d</sup>*<sup>4</sup> <sup>≤</sup> *<sup>δ</sup>*4*, the positive equilibrium solution* **u** = (*u*, *v*,*r*)*<sup>T</sup> of Equation (26) satisfies,*

$$\begin{aligned} \max\_{\Omega} (u, v, r) &\leq \bar{Q}(\Gamma, \delta\_1, \delta\_2, \delta\_3, \delta\_4) \\ \max\_{\Omega} (u, v, r) &\leq Q^\*(\Gamma, \delta\_1, \delta\_2, \delta\_3, \delta\_4) \min\_{\Omega} (u, v, r) \end{aligned} \tag{27}$$

**Proof.** We apply the maximum principle to the first part of Equation (26) and get max <sup>Ω</sup>¯ *<sup>u</sup>* <sup>≤</sup> 1. Similarly, applying the maximum principle to the second part of Equation (26) gives max <sup>Ω</sup>¯ *<sup>v</sup>* <sup>≤</sup> *<sup>w</sup>*<sup>6</sup> .

*e*(1 + *w*7) Let *<sup>ϕ</sup>*(*x*) = *<sup>d</sup>*3(*d*4*vr* <sup>+</sup> *<sup>r</sup>*). Let *<sup>x</sup>*<sup>1</sup> <sup>∈</sup> <sup>Ω</sup>¯ be such that *<sup>ϕ</sup>*(*x*1) maxΩ¯ *<sup>ϕ</sup>*. Now, applying the maximum principle to the third part of Equation (26) gives,*r*(*x*1) <sup>≤</sup> (*w*<sup>12</sup> <sup>−</sup> *<sup>w</sup>*11)*w*<sup>6</sup> <sup>−</sup> *<sup>e</sup>*(<sup>1</sup> <sup>+</sup> *<sup>w</sup>*7)*w*10*w*<sup>11</sup> *e*(1 + *w*7)*w*8*w*<sup>11</sup> + *w*6*w*9*w*<sup>11</sup> .

Defining *ϕ*(*x*) as *d*<sup>4</sup> ≤ *δ*4, we have,

$$\begin{split} \max\_{\Omega} r \leq \frac{1}{d\_3} \max\_{\Omega} \varphi(\mathbf{x}\_1) &= \frac{1}{d\_3} \varphi(\mathbf{x}\_1) = r(\mathbf{x}\_1) + d\_4 v(\mathbf{x}\_1) r(\mathbf{x}\_1) \\ &\leq r(\mathbf{x}\_1) + d\_4 \max\_{\Omega} v(\mathbf{x}\_1) r(\mathbf{x}\_1) \\ &\leq r(\mathbf{x}\_1) \left( 1 + d\_4 \max\_{\Omega} v(\mathbf{x}\_1) \right) \\ &\leq \left( 1 + \delta\_4 \frac{w\_6}{\varepsilon (1 + w\_7)} \right) \left( \frac{(w\_{12} - w\_{11}) w\_6 - \varepsilon (1 + w\_7) w\_{10} w\_{11}}{\varepsilon (1 + w\_7) w\_8 w\_{11} + w\_6 w\_9 w\_{11}} \right) \end{split}$$

Let,

$$\overline{Q}(\Gamma,d) = \max\left\{1, \frac{w\_6}{\varepsilon(1+w\gamma)}, \left(1+\delta\_4 \frac{w\_6}{\varepsilon(1+w\gamma)}\right) \left(\frac{(w\_{12}-w\_{11})w\_6 - \varepsilon(1+w\gamma)w\_{10}w\_{11}}{\varepsilon(1+w\gamma)w\_8w\_{11} + w\_6w\_9w\_{11}}\right)\right\}$$

Therefore, we have,

$$\max\_{\Omega} (u, v, r) \le \overline{\mathbb{Q}}(\Gamma, \delta\_1, \delta\_2, \delta\_3, \delta\_4)$$

Now, we show that,

$$\max\_{\vec{\Omega}} (u, v, r) \le Q^\*(\Gamma, \delta\_1, \delta\_2, \delta\_3, \delta\_4) \min\_{\vec{\Omega}} (u, v, r).$$

Using boundedness of *u*, *v* and *r*, we apply the Harnack inequality to the first and second parts of Equation (26), yielding

$$\max\_{\Omega}(\mathfrak{u}, \mathfrak{v}) \le Q^\*(\Gamma, \delta\_1, \delta\_2, \delta\_3, \delta\_4) \min\_{\Omega}(\mathfrak{u}, \mathfrak{v}),$$

Let *ϑ*(*x*) = *d*3(*d*4*vr* + *r*) and we have,

$$\begin{cases} \Delta \theta + c(x)\theta = 0, & x \in \Omega, \\\ \partial\_{\mathfrak{n}} \theta = 0, & x \in \partial \Omega. \end{cases}$$

where *<sup>c</sup>*(*x*) = (*w*<sup>12</sup> <sup>−</sup> *<sup>w</sup>*11)*<sup>v</sup>* <sup>−</sup> (*w*8*w*<sup>11</sup> <sup>+</sup> *<sup>w</sup>*9*w*11*v*)*<sup>r</sup>* <sup>−</sup> *<sup>w</sup>*10*w*<sup>11</sup> *<sup>d</sup>*3(<sup>1</sup> <sup>+</sup> *<sup>d</sup>*4*v*)(*<sup>v</sup>* + (*w*<sup>8</sup> <sup>+</sup> *<sup>w</sup>*9*v*)*<sup>r</sup>* <sup>+</sup> *<sup>w</sup>*10) is bounded. By Harnack inequality,

$$\max\_{\Omega} \theta \le Q\_3^\* \min\_{\Omega} \theta$$

for some positive constant

$$Q\_3^\* = Q\_3^\*(\Gamma, \delta\_{1\prime}\delta\_{2\prime}\delta\_{3\prime}\delta\_4).$$

Ω¯

**Lemma 1.** *Let di*,*<sup>m</sup> i* = 1, 2, 3, 4*, be positive constants, m* = 1, 2, ... , *and* **um** = (*um*, *vm*,*rm*)*<sup>T</sup> be the corresponding positive equilibrium solution of Equation (26) with di* = *di*,*m. If* **um** → *u*¯ *as m* → ∞ *and* **u¯** *is a constant vector, then* **u¯** = **u**∗*, where u*∗ *is the nontrivial solution of* **G**(**u**) = 0*.*

**Proof.** For all *m*, Ω *um g*1(**um**)*dx* = 0. If *g*1(**u¯**) > 0, then *g*1(**um**) > 0 when *m* is large. However, since *um* is positive, this is not possible. Similarly, *g*1(**u**¯) < 0 is not possible. Therefore, *g*1(**u¯**) = 0. The same argument shows that *g*2(**u¯**) = *g*3(**u¯**) = 0.

**Theorem 5.** *Lower Bound: Let δ*1*, δ*2*, δ*<sup>3</sup> *and δ*<sup>4</sup> *be fixed positive constants. There exists positive constants <sup>Q</sup>*(Γ, *<sup>δ</sup>*1, *<sup>δ</sup>*2, *<sup>δ</sup>*3, *<sup>δ</sup>*4) *such that, when <sup>d</sup>*<sup>1</sup> ≥ *<sup>δ</sup>*1*, <sup>d</sup>*<sup>2</sup> ≥ *<sup>δ</sup>*2*, <sup>d</sup>*<sup>3</sup> ≥ *<sup>δ</sup>*<sup>3</sup> *and <sup>d</sup>*<sup>4</sup> ≤ *<sup>δ</sup>*4*, the positive solution* **<sup>u</sup>** = (*u*, *<sup>v</sup>*,*r*)*<sup>T</sup> of Equation (26) satisfies,*

$$\min\_{\Omega}(u,v,r) \ge \underline{Q}(\Gamma, \delta\_1, \delta\_2, \delta\_3, \delta\_4) \tag{28}$$

**Proof.** We integrate the second part of Equation (26) in Ω and consider the inhomogeneous Neumann boundary condition, yielding

$$\int\_{\Omega} v \left( -w\_5 + \frac{w\_6 u}{u + w\_7} - \epsilon v - \frac{r}{v + (w\_8 + w\_9 v)r + w\_{10}} \right) = 0$$

therefore, there exists *<sup>x</sup>*<sup>0</sup> <sup>∈</sup> <sup>Ω</sup> such that *<sup>w</sup>*6*u*(*x*0) *u*(*x*0) + *w*<sup>7</sup> = *w*<sup>5</sup> + *ev*(*x*0) + *r*(*x*0)

*v*(*x*0)+(*w*<sup>8</sup> + *w*9*v*(*x*0))*r*(*x*0) + *w*<sup>10</sup> . As *<sup>u</sup>* <sup>≤</sup> 1, *<sup>u</sup>*(*x*0) <sup>≥</sup> *<sup>w</sup>*5*w*<sup>7</sup> *w*<sup>6</sup> . By using Harnack inequality,

$$\min\_{\Omega} u(\mathfrak{x}) \ge \frac{1}{Q^\*} \frac{w\_5 w\_7}{w\_6}$$

Now, integrating third part of Equation (26) in Ω

$$\min\_{\Omega} r(x) \ge \frac{1}{Q^\*} w\_{10} w\_{12} \dots$$

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

With reference to Equation (27), to prove Theorem 5, it is sufficient to show that,

$$\max\_{\Omega} v \ge \underline{Q\_1}(\Gamma, \delta\_1, \delta\_2, \delta\_3, \delta\_4). \tag{29}$$

Let the inequality in Equation (29) not hold. Then, there exists sequences {*d*1,*m*, *<sup>d</sup>*2,*m*, *<sup>d</sup>*3,*m*, *<sup>d</sup>*4,*m*}<sup>∞</sup> *m*=1 with *d*1,*<sup>m</sup>* ≥ *δ*1, *d*2,*<sup>m</sup>* ≥ *δ*2, *d*3,*<sup>m</sup>* ≥ *δ*3, *d*4,*<sup>m</sup>* ≤ *δ*<sup>4</sup> and the corresponding positive solutions **um** to Equation (26) such that max <sup>Ω</sup>¯ *vm* <sup>→</sup> 0, where, **um** = (*um*, *vm*,*rm*)*<sup>T</sup>* satisfies

$$\begin{aligned} -d\_{1,m}\Delta u\_m &= u\_m \left( 1 - u\_m - \frac{v\_m}{u\_m + w\_4} \right) \\ -d\_{2,m}\Delta v\_m &= v\_m \left( -w\_5 + \frac{w\_6\mu}{u\_m + w\_7} - cv\_m - \frac{r\_m}{v\_m + \left( w\_8 + w\_9 v\_m \right) r\_m + w\_{10}} \right) \\ -d\_{3,m}\Delta \left( r\_m + d\_4 v\_m r\_m \right) &= r\_m \left( -w\_{11} + \frac{w\_{12} v\_m}{v\_m + \left( w\_8 + w\_9 v\_m \right) r\_m + w\_{10}} \right) \\ \frac{\partial u\_m}{\partial n} &= \frac{\partial v\_m}{\partial n} = \frac{\partial r\_m}{\partial n} = 0. \end{aligned} \tag{30}$$

Since maxΩ¯ (*u*,*r*) > 0, we may assume that maxΩ¯ *um* → *u*¯ and maxΩ¯ *rm* → *r*¯, where *u*¯ and *r*¯ are positive constants. In addition, we claim that *um* and *rm* converge uniformly to positive constants, respectively. In fact, there are two cases of {*d*3,*m*}<sup>∞</sup> *<sup>m</sup>*=<sup>1</sup> to be considered.

**Case 1** {*d*3,*m*}<sup>∞</sup> *<sup>m</sup>*=<sup>1</sup> is bounded with respect to *m*. Set,

$$
\eta\_m = d\_{3,m}(r\_m + d\_{4,m}\upsilon\_m r\_m).
$$

Using upper bound of **um**, we have *ηm* ≤ *C* for all *m* ≥ 1. Since *η<sup>m</sup>* satisfies

$$\begin{aligned} \Delta \eta\_{\text{m}} + \frac{(w\_{12} - w\_{11})v\_{\text{m}} - (w\_{8}w\_{11} + w\_{9}w\_{11}v\_{\text{m}})r\_{\text{m}} - w\_{10}w\_{11}}{d\_{3}(1 + d\_{4}v\_{\text{m}})(v\_{\text{m}} + (w\_{8} + w\_{9}v\_{\text{m}})r\_{\text{m}} + w\_{10})} &= 0 \\ \partial\_{\text{n}} \eta\_{\text{m}} = 0 \end{aligned}$$

using *L<sup>p</sup>* estimate and the Sobolev embedding theorems, we have

⎧ ⎨ ⎩

$$||\eta\_{m}||\_{Q^{1,a}(\Omega)} \leq Q||\eta\_{m}||\_{W^{2,p}\_{2,p}(\Omega)} \leq Q\_{2}$$

Similarly, the *C*2,*α*(Ω¯ ) norms of *η<sup>m</sup>* is uniformly bounded with respect to *m*. Thus, we assume that *<sup>η</sup><sup>m</sup>* <sup>→</sup> *<sup>η</sup>* in *<sup>C</sup>*2(Ω¯ ). By the definition of *<sup>η</sup>m*, for sufficiently large *<sup>m</sup>* we have,

$$r\_m = \frac{\eta\_m}{d\_3(1 + d\_{4,m}v\_m)}$$

Since *vm* → 0 and *d*4,*<sup>m</sup>* ≤ *δ*4, we have

$$r\_m \to r \equiv \frac{\eta}{d\_3} \cdot \frac{\eta}{\cdot}$$

Hence, *η* satisfies

$$\begin{aligned} d\_3^2 \Delta \eta + \eta \left( c d\_3 - \frac{\eta}{w\_{10}} \right) &= 0 \\ \partial\_n \eta &= 0 \end{aligned}$$

Hence, *<sup>η</sup>* <sup>≡</sup> *cw*10*d*3; otherwise, *<sup>η</sup>* <sup>≡</sup> 0 on <sup>Ω</sup>¯ that implies *<sup>r</sup>* <sup>=</sup> *<sup>r</sup>*¯ <sup>=</sup> 0, which is a contradiction to *<sup>r</sup>*¯ <sup>&</sup>gt; 0, therefore *<sup>r</sup>* <sup>≡</sup> *<sup>η</sup> d*3 = *cw*10.

**Case 2** {*d*3,*m*}<sup>∞</sup> *<sup>m</sup>*=<sup>1</sup> is unbounded with respect to *m*. We may assume that *d*3,*<sup>m</sup>* → ∞. Set

$$
\psi\_{\mathfrak{m}} = r\_{\mathfrak{m}} + d\_{\mathfrak{k}} v\_{\mathfrak{m}} r\_{\mathfrak{m}}.
$$

Therefore, proceeding as above

$$-\Delta \psi = 0 \quad \text{in} \quad \Omega \qquad \qquad \qquad \partial\_n \psi = 0 \quad \text{on} \quad \partial \Omega.$$

Hence, *ψ* = *r* ≡ constant > 0. Similarly, we can prove that *u* = constant ≥ 0. The above argument shows that there are positive constants *c*1, *c*<sup>3</sup> ≥ 0, such that (*um*, *vm*,*rm*) → (*c*1, 0, *c*3). This is contradiction to Lemma 5.1, thus the proof is completed.

#### *5.2. Existence of Inhomogeneous Positive Steady States*

We, now discuss the existence of inhomogeneous positive solutions. Let **X** be defined as in Notation 1. Define

$$\begin{aligned} \mathsf{X}^+ &= \{ \mathbf{u} \in \mathsf{X} : u > 0, v > 0, r > 0 \text{ on } \bar{\Omega} \} \\ B(Q) &= \{ \mathbf{u} \in \mathsf{X} : Q^{-1} < u, v, r < Q \text{ on } \bar{\Omega} \}, \; \mathsf{C} > 0. \end{aligned}$$

Let Φ*<sup>u</sup>* = (*d*1*u*, *d*2*v*, *d*3*r*(1 + *d*4*v*))*T*. Then, Equation (26) is transformed as follows,

$$\begin{cases} -\Delta \Phi(\mathbf{u}) = G(\mathbf{u}), & \mathbf{x} \in \Omega, \\\ \partial\_n \mathbf{u} = 0, & \mathbf{x} \in \Omega. \end{cases} \tag{31}$$

As the determinant of Φ(**u**) is positive for all non-negative **u**, Φ−1(**u**) exists and det Φ−<sup>1</sup> **<sup>u</sup>** is positive. Hence, **u** is a positive solution to Equation (31) if and only if

$$F(\mathbf{u}) \cong \mathbf{u} - (I - \Delta)^{-1} \{ \Phi\_{\boldsymbol{\mu}}^{-1} [\mathbf{G}(\mathbf{u}) + \nabla \mathbf{u}. \nabla \mathbf{u}] + \mathbf{u} \} = 0 \quad \text{in} \quad \mathcal{X}^+$$

where (*<sup>I</sup>* − <sup>Δ</sup>)−<sup>1</sup> is the inverse of (*<sup>I</sup>* − <sup>Δ</sup>) in **<sup>X</sup>**, where *<sup>F</sup>*(.) is a compact perturbation of the identity operator, for any *B* = *B*(*Q*). the Leray-Schauder degree deg(*F*(.), 0, *B*) is well defined if *F*(**u**) = 0 on *∂B*.

Furthermore, we note that

$$D\_{\mathbf{u}}F(\mathbf{u}^\*) = I - (I - \Delta)^{-1} \left\{ \Phi\_{\mathbf{u}}^{-1}(\mathbf{u}^\*) \mathbf{G}\_{\mathbf{u}}(\mathbf{u}^\*) + I \right\},$$

and as *DuF*(**u**∗) is invertible, the index of *<sup>F</sup>* at **<sup>u</sup>**<sup>∗</sup> is defined as index(*F*(.), **<sup>u</sup>**∗)=(−1)*γ*, where *<sup>γ</sup>* is the number of negative eigenvalues of *DuF*(**u**∗).

For every integer *i* ≥ 1 and for each integer 1 ≤ *j* ≤ dim*E*(*μi*), **Xij** is invariant under *DuF*(**u**∗), and *λ* is an eigenvalue of *DuF*(**u**∗) on **Xij** if and only if it is an eigenvalue of the matrix,

$$I - \frac{1}{1 + \mu\_i} \left[ \Phi\_{\boldsymbol{\mu}}^{-1} \mathbf{u}^\* \mathbf{G} \mathbf{u} (\mathbf{u}^\*) \right] = \frac{1}{1 + \mu\_i} \left[ \mu\_i I - \Phi\_{\boldsymbol{\mu}}^{-1} \mathbf{u}^\* \mathbf{G} \mathbf{u} (\mathbf{u}^\*) \right]$$

Thus, *DuF*(**u**∗) is invertible if and only if, for all *<sup>i</sup>* <sup>≥</sup> 1, the matrix *<sup>I</sup>* <sup>−</sup> <sup>1</sup> 1 + *μ<sup>i</sup>* 3 Φ−<sup>1</sup> *<sup>u</sup>* **u**∗**Gu**(**u**∗) 4 is non-singular. We have,

$$H(\mu) = H(\mathbf{u}^\*; \mu) \cong \det \left\{ \mu I - \Phi\_\mu^{-1} \mathbf{u}^\* \mathbf{G}\_\mathbf{u} \left( \mathbf{u}^\* \right) \right\},$$

Further, if *H*(*μi*) = 0, then for each 1 ≤ *j* ≤ dim*E*(*μi*), the number of negative eigenvalues of *DuF*(**u**∗) on **Xij** is odd if and only if *H*(*μi*) < 0. We have results as follows.

**Proposition 3.** *Suppose that for all, i* <sup>≥</sup> <sup>1</sup>*, the matrix <sup>μ</sup>iI* <sup>−</sup> <sup>Φ</sup>−<sup>1</sup> *<sup>i</sup>* (**u**∗)*Gu*(**u**∗) *is nonsingular. Then,*

$$index(F(.), \mathbf{u}^\*) = (-1)^\tau, \quad where \quad \tau = \sum\_{\mathbf{i} \ge \mathbf{1}, \mathbf{H}(\tau\_\mathbf{i}) < \mathbf{0}} dim E(\mu\_\mathbf{i}) .$$

To compute *index*(*F*(.), **u**∗), we consider the sign of *H*(*μi*). As we want to study the existence of positive solutions of Equation (31) with respect to *d*2, we focus on dependence of *H*(*μi*) on *d*2, and consider *d*1, *d*<sup>3</sup> and *d*<sup>4</sup> to be fixed. We denote,

$$H(\mu) = \det \left\{ \Phi\_{\boldsymbol{\mu}}^{-1}(\mathbf{u}^\*) \right\} \det \left\{ \mu \Phi\_{\boldsymbol{\mu}}(\mathbf{u}^\*) - \mathbf{G}\_{\mathbf{u}}(\mathbf{u}^\*) \right\}.$$

As det ' Φ−<sup>1</sup> *<sup>u</sup>* (**u**∗) ( <sup>&</sup>gt; 0, we only need to consider det {*μ*Φ*u*(**u**∗) <sup>−</sup> **Gu**(**u**∗)}. In fact *<sup>B</sup>*<sup>3</sup> <sup>=</sup> det {*μ*Φ*u*(**u**∗) − **Gu**(**u**∗)}, and *B*<sup>3</sup> < 0 by Theorem 3. Therefore, we make the following proposition.

**Proposition 4.** *Assuming the condition specified in Theorem 3 holds, then there exists a positive number d*∗ <sup>3</sup> *such that, for all d*<sup>3</sup> > *d*<sup>∗</sup> <sup>3</sup>*, the roots μ*˜1, *μ*˜2, *μ*˜3 *of* det {*μ*Φ*u*(**u**∗) − **Gu**(**u**∗)} = 0 *are all real and satisfy Equations (23)–(25). Moreover, for all d*<sup>3</sup> ≥ *d*<sup>∗</sup> 3*.*

$$\begin{cases} \begin{aligned} & -\infty < \tilde{\mu}\_{1} < 0 < \tilde{\mu}\_{2} < \tilde{\mu}\_{3}, \\ & \det\left\{\mu \Phi\_{\boldsymbol{\upmu}}(\mathbf{u}^{\ast}) - \mathbf{G}\_{\mathbf{u}}(\mathbf{u}^{\ast})\right\} < 0, \quad \text{when} \quad \boldsymbol{\upmu} \in (\infty, \tilde{\mu}\_{1}) \cup (\tilde{\mu}\_{2}, \tilde{\mu}\_{3}), \\ & \det\left\{\mu \Phi\_{\boldsymbol{\upmu}}(\mathbf{u}^{\ast}) - \mathbf{G}\_{\mathbf{u}}(\mathbf{u}^{\ast})\right\} > 0, \quad \text{when} \quad \boldsymbol{\upmu} \in (\bar{\mu}\_{1}, \bar{\mu}\_{2}) \cup (\bar{\mu}\_{3}, \infty). \end{aligned} \end{cases} \tag{32}$$

To discuss the effect of cross diffusion on the existence of inhomogeneous positive solution to Equation (31), we first deduce a non-existence result when the cross diffusion term is absent.

**Theorem 6.** *Suppose that*

$$d\_1 \ge \frac{1}{\mu\_2}, \qquad \qquad d\_3 \ge \frac{w\_{11}}{w\_{12}\mu\_2}$$

*where μ*<sup>2</sup> *is given previously. Then, there exists a positive constants D*<sup>2</sup> *when d*<sup>2</sup> ≥ *D*<sup>2</sup> *without cross diffusion term d*4*, Equation (26) has no non-constant positive solution. Furthermore,*

$$\begin{aligned} (u, v, r) &= (\vec{u}, \vec{v}, \vec{r}), \\ \text{where} \quad \vec{u} &\cong \frac{1}{\text{measure}(\Omega)} \int\_{\Omega} u; \quad \vec{v} \cong \frac{1}{\text{measure}(\Omega)} \int\_{\Omega} v; \quad \vec{r} \cong \frac{1}{\text{measure}(\Omega)} \int\_{\Omega} r. \end{aligned}$$

**Proof.** Assume **u** = (*u*, *v*,*r*) is a positive solution of Equation (26) with *d*<sup>4</sup> = 0. Let

$$\text{all} \cong \frac{1}{\text{measure}(\Omega)} \int\_{\Omega} u; \quad \vartheta \cong \frac{1}{\text{measure}(\Omega)} \int\_{\Omega} v; \quad \text{\textdegree} \quad \text{\textdegree} \quad \frac{1}{\text{measure}(\Omega)} \int\_{\Omega} v$$

Multiplying the first part of Equation (26) by (*u* − *u*¯), the second part by (*v* − *v*¯) and the third part by (*r* − *r*¯) (taking *d*<sup>4</sup> = 0) and integrating by parts over Ω, we have

*d*1 Ω |∇*u*| <sup>2</sup> = Ω (*u* − *u*¯)(*ug*1(*u*, *v*) − *ug*¯ <sup>1</sup>(*u*¯, *v*¯)) = Ω (*<sup>u</sup>* <sup>−</sup> *<sup>u</sup>*¯)<sup>2</sup> <sup>1</sup> <sup>−</sup> *<sup>u</sup>* <sup>−</sup> *<sup>u</sup>*¯ <sup>−</sup> *<sup>w</sup>*4*<sup>v</sup>* (*u* + *w*4)(*u*¯ + *w*4) − Ω (*u* − *u*¯)(*v* − *v*¯)(*w*4*u*¯ + *uu*¯) *d*2 Ω |∇*v*| <sup>2</sup> = Ω (*v* − *v*¯)(*vg*2(*u*, *v*,*r*) − *vg*¯ <sup>2</sup>(*u*¯, *v*¯,*r*¯)) = Ω (*<sup>v</sup>* <sup>−</sup> *<sup>v</sup>*¯)2<sup>×</sup> *w*6*uu*¯ (*<sup>u</sup>* <sup>+</sup> *<sup>w</sup>*7)(*u*¯ <sup>+</sup> *<sup>w</sup>*7) <sup>−</sup> *<sup>w</sup>*<sup>5</sup> <sup>−</sup> *<sup>e</sup>*(*<sup>v</sup>* <sup>+</sup> *<sup>v</sup>*¯) <sup>−</sup> *<sup>w</sup>*8*rr*¯ [*v* + (*w*<sup>8</sup> + *w*9*r*) *v* + *w*10] [*v*¯ + (*w*<sup>8</sup> + *w*9*r*¯) *v*¯ + *w*10] + Ω *w*6*w*7(*vu* − *v*¯*u*¯)(*v* − *v*¯) (*<sup>u</sup>* <sup>+</sup> *<sup>w</sup>*7)(*u*¯ <sup>+</sup> *<sup>w</sup>*7) <sup>−</sup> Ω *w*10(*v*¯*r*¯ − *rv*)(*v* − *v*¯) (*u* + *w*7)(*u*¯ + *w*7) + Ω *vv*¯(*v* − *v*¯)(*r* − *r*¯) [*v* + (*w*<sup>8</sup> + *w*9*r*) *v* + *w*10] [*v*¯ + (*w*<sup>8</sup> + *w*9*r*¯) *v*¯ + *w*10]

$$\begin{split} d\_3 \int\_{\Omega} |\nabla r|^2 &= \int\_{\Omega} (r-\mathfrak{r}) (r\mathfrak{g}\_3(\mathfrak{v},r) - r\mathfrak{g}\_3(\mathfrak{v},\mathfrak{r})) \\ &= \int\_{\Omega} (r-\mathfrak{r})^2 \left[ \frac{w\_{12}\upsilon\mathfrak{v}}{\left[\upsilon + (w\_8 + w\_9\mathfrak{r})\upsilon + w\_{10}\right] \left[\vartheta + (w\_8 + w\_9\mathfrak{r})\vartheta + w\_{10}\right]} \right] \\ &+ \int\_{\Omega} \frac{w\_{12}\upsilon\mathfrak{v}(r+\mathfrak{r})}{\left[\upsilon + (w\_8 + w\_9\mathfrak{r})\upsilon + w\_{10}\right] \left[\vartheta + (w\_8 + w\_9\mathfrak{r})\vartheta + w\_{10}\right]} \\ &+ \int\_{\Omega} \frac{w\_{10}w\_{12}(r\upsilon - \mathfrak{r}\mathfrak{r})(r-\mathfrak{r})}{\left[\upsilon + (w\_8 + w\_9\mathfrak{r})\upsilon + w\_{10}\right] \left[\vartheta + (w\_8 + w\_9\mathfrak{r})\vartheta + w\_{10}\right]} \end{split}$$

By using Cauchy inequality, we have *d*<sup>1</sup> Ω |∇*u*| <sup>2</sup> + *d*<sup>2</sup> Ω |∇*v*| <sup>2</sup> + *d*<sup>3</sup> Ω |∇*r*| 2 ≤ Ω (*<sup>u</sup>* <sup>−</sup> *<sup>u</sup>*¯)<sup>2</sup> <sup>1</sup> <sup>−</sup> *<sup>u</sup>* <sup>−</sup> *<sup>u</sup>*¯ <sup>−</sup> *<sup>w</sup>*4*<sup>v</sup>* (*<sup>u</sup>* <sup>+</sup> *<sup>w</sup>*4)(*u*¯ <sup>+</sup> *<sup>w</sup>*4) <sup>+</sup>  + Ω (*<sup>v</sup>* <sup>−</sup> *<sup>v</sup>*¯)<sup>2</sup> <sup>×</sup> *w*6*uu*¯ (*<sup>u</sup>* <sup>+</sup> *<sup>w</sup>*7)(*u*¯ <sup>+</sup> *<sup>w</sup>*7) <sup>−</sup> *<sup>w</sup>*<sup>5</sup> <sup>−</sup> *<sup>e</sup>*(*<sup>v</sup>* <sup>+</sup> *<sup>v</sup>*¯) <sup>−</sup> *<sup>w</sup>*8*rr*¯ [*v* + (*w*<sup>8</sup> + *w*9*r*) *v* + *w*10] [*v*¯ + (*w*<sup>8</sup> + *w*9*r*¯) *v*¯ + *w*10] + Ω (*<sup>r</sup>* <sup>−</sup> *<sup>r</sup>*¯)<sup>2</sup> *w*12*vv*¯ [*v* + (*w*<sup>8</sup> + *w*9*r*) *v* + *w*10] [*v*¯ + (*w*<sup>8</sup> + *w*9*r*¯) *v*¯ + *w*10] + 1 4 Ω *w*6*w*7(*vu* − *v*¯*u*¯)(*v* − *v*¯) (*u* + *w*7)(*u*¯ + *w*7) <sup>+</sup>*vv*¯(*<sup>r</sup>* <sup>−</sup> *<sup>r</sup>*¯)(*<sup>v</sup>* <sup>−</sup> *<sup>v</sup>*¯) + *<sup>w</sup>*12*vv*¯(*<sup>r</sup>* <sup>+</sup> *<sup>r</sup>*¯) + *<sup>w</sup>*10*w*12(*rv* <sup>−</sup> *<sup>r</sup>*¯*v*¯)(*<sup>r</sup>* <sup>−</sup> *<sup>r</sup>*¯) <sup>−</sup> *<sup>w</sup>*10(*r*¯*v*¯ <sup>−</sup> *rv*)(*<sup>v</sup>* <sup>−</sup> *<sup>v</sup>*¯) [*v* + (*w*<sup>8</sup> + *w*9*r*) *v* + *w*10] [*v*¯ + (*w*<sup>8</sup> + *w*9*r*¯) *v*¯ + *w*10] −(*u* − *u*¯)(*v* − *v*¯)(*uw*¯ <sup>4</sup> + *uu*¯) ,

where  is any arbitrary positive constant.

By Poincare inequality we have,

$$\begin{aligned} &d\_1 \int\_{\Omega} |\nabla u|^2 + d\_2 \int\_{\Omega} |\nabla v|^2 + d\_3 \int\_{\Omega} |\nabla r|^2 \\ &\ge \int\_{\Omega} d\_1 \mu\_2 (u - \bar{u})^2 + \int\_{\Omega} d\_2 \mu\_2 (v - \bar{v})^2 + \int\_{\Omega} d\_3 \mu\_2 (r - \bar{r})^2. \text{ By Theorem 5.1, we can choose a sufficiently small } \epsilon \text{ such that,} \\ &\text{ sufficiently small } \epsilon\_0 \text{ such that,} \end{aligned}$$

$$\begin{aligned} d\_1 \mu\_2 &> 1 - \mathfrak{u} - \overline{\mathfrak{u}} - \frac{w\_4 \upsilon}{\left(\mathfrak{u} + w\_4\right) \left(\overline{\mathfrak{u}} + w\_4\right)} + \mathfrak{e}\_0\\ d\_3 \mu\_2 &> \frac{w\_{12} \upsilon \overline{\upsilon}}{\left[\upsilon + \left(w \mathfrak{s} + w \upsilon r\right) \upsilon + w\_{10}\right] \left[\bar{\upsilon} + \left(w \mathfrak{s} + w \upsilon \bar{r}\right) \bar{\upsilon} + w\_{10}\right]} - \mathfrak{w}\_{11} + \mathfrak{e}\_1 \end{aligned}$$

Lastly, by taking

*D*¯ <sup>2</sup> > 1 *μ*2 *w*6*uu*¯ (*<sup>u</sup>* <sup>+</sup> *<sup>w</sup>*7)(*u*¯ <sup>+</sup> *<sup>w</sup>*7) <sup>−</sup> *<sup>w</sup>*<sup>5</sup> <sup>−</sup> *<sup>e</sup>*(*<sup>v</sup>* <sup>+</sup> *<sup>v</sup>*¯) <sup>−</sup> *<sup>w</sup>*8*rr*¯ [*v* + (*w*<sup>8</sup> + *w*9*r*) *v* + *w*10] [*v*¯ + (*w*<sup>8</sup> + *w*9*r*¯) *v*¯ + *w*10] + 1 4 *w*6*w*7(*uv* − *u*¯*v*¯)(*v* − *v*¯) (*<sup>u</sup>* <sup>+</sup> *<sup>w</sup>*7)(*u*¯ <sup>+</sup> *<sup>w</sup>*7) <sup>+</sup> *vv*¯(*<sup>r</sup>* <sup>−</sup> *<sup>r</sup>*¯)(*<sup>v</sup>* <sup>−</sup> *<sup>v</sup>*¯) [*v* + (*w*<sup>8</sup> + *w*9*r*) *v* + *w*10] [*v*¯ + (*w*<sup>8</sup> + *w*9*r*¯) *v*¯ + *w*10] <sup>+</sup> *<sup>w</sup>*12*vv*¯(*<sup>r</sup>* <sup>+</sup> *<sup>r</sup>*¯) + *<sup>w</sup>*10(*rv* <sup>−</sup> *<sup>r</sup>*¯*v*¯)(*w*<sup>12</sup> + (*<sup>v</sup>* <sup>−</sup> *<sup>v</sup>*¯)) [*v* + (*w*<sup>8</sup> + *w*9*r*) *v* + +*w*10] [*v*¯ + (*w*<sup>8</sup> + *w*9*r*¯) *v*¯ + *w*10] + (*v* − *v*¯)(*u*¯ − *u*)(*uw*¯ <sup>4</sup> + *uu*¯) 

We conclude that (*u*, *v*,*r*)=(*u*¯, *v*¯,*r*¯), which completes the proof. Therefore, we conclude that the model system has no non-constant positive solution. Now, we prove that in the presence of cross diffusion inhomogeneous solution is generated. The following theorem discusses the global existence of of inhomogeneous solution of Equation (8) with respect to *d*<sup>3</sup> as the other parameters *d*1, *d*<sup>2</sup> and *d*<sup>4</sup> are fixed.

**Theorem 7.** *Let d*1, *d*2, *d*<sup>4</sup> *be fixed and satisfy conditions of Theorems 2 and 3. Let μ*˜ *be defined in Equations (23)–(25). If μ*˜ ∈ (*μn*, *μn*+1) *for some n* ≥ 1*, and the sum σ<sup>n</sup>* = *n* ∑ *i*=2 *dimEμ<sup>i</sup> is odd. Then, there exist a positive number d*∗ <sup>3</sup> *such that, if d*<sup>3</sup> > *d*<sup>∗</sup> <sup>3</sup>*, Equation (8) has at least one inhomogeneous positive steady state solution.*

**Proof.** Let *δ*1, *δ*2, *δ*<sup>3</sup> and *δ*<sup>4</sup> be positive constants and satisfy *δ*<sup>1</sup> < *d*1, *δ*<sup>2</sup> < *d*2, *δ*<sup>3</sup> < *d*<sup>3</sup> and *δ*<sup>4</sup> > *d*4. By Equations (23)–(25) and Proposition 4, there exists a positive constant *d*∗ <sup>3</sup> such that, when *δ*<sup>3</sup> > *d*<sup>∗</sup> 3, Equation (32) holds and

$$0 = \mu\_1 < \overline{\mu}\_2 < \mu\_2, \qquad \overline{\mu}\_3 \in \left(\mu\_n, \mu\_{n+1}\right) \tag{33}$$

Now, we prove that for, *δ*<sup>3</sup> > *d*<sup>∗</sup> <sup>3</sup>, Equation (26) has at least one inhomogeneous positive solution. The statement has been proved by contradiction and is based on homotopy invariance of topological degree. On the contrary, we assume that the assertion is not true. For *t* ∈ [0, 1], define

$$\Phi(t; \mathbf{u}) = (\hat{d}u + t(d - 1 - \hat{d}\_1)u, \hat{d}\_2v + t(d\_2 - \hat{d}\_2)v \times (\hat{d}\_3 + t(\hat{d}\_3 - d\_3))(r + td\_4vr))^T$$

where, <sup>ˆ</sup>*d*<sup>2</sup> <sup>=</sup> *<sup>D</sup>*¯ 2, <sup>ˆ</sup>*d*<sup>1</sup> <sup>=</sup> <sup>1</sup> *μ*2 , <sup>ˆ</sup>*d*<sup>3</sup> <sup>=</sup> *<sup>c</sup> μ*2 , and consider

$$\begin{cases} -\Delta\Phi(t,\mathbf{u}) = \mathbf{G}(\mathbf{u}) & \text{in}\Omega, \quad 0 \le t \le 1\\ \partial\_{\mathbf{u}}\mathbf{u} = 0 & \text{on} \quad \partial\Omega \end{cases} \tag{34}$$

Then, **u** is a positive inhomogeneous solution of Equation (26) if and only if it is a solution of Equation (34) for *t* = 1. **u**<sup>∗</sup> is the unique constant solution of Equation (34) for any 0 ≤ *t* ≤ 1. Now, we have for any 0 ≤ *t* ≤ 1, **u** is a unique positive solution of Equation (34) if and only if,

$$F(t; \mathbf{u}) \cong \mathbf{u} - (I - \Delta)^{-1} \left\{ \Phi\_{\boldsymbol{\mu}}^{-1} [\mathbf{G}(\mathbf{u}) + \nabla \mathbf{u}. \Phi\_{\boldsymbol{\mu}\mathbf{u}}(t; \mathbf{u}). \nabla \mathbf{u}] + \mathbf{u} \right\} = 0 \quad \text{in} \quad \mathbb{X}^+$$

As *F*(1; **u**) = *F*(**u**), by Theorem 4.2, **F**(**0**; **u**) = 0 has only positive solution **u**<sup>∗</sup> in *X*+, therefore we have,

$$\mathcal{D}\_{\mathfrak{u}}F(0,\mathfrak{u}^\*) = I - (I - \Delta)^{-1} \left\{ \mathcal{D}^{-1} \mathbf{G}\_{\mathfrak{u}}(\mathfrak{u}^\*) + 1 \right\}$$

$$\mathcal{D}\_{\mathfrak{u}}F(1,\mathfrak{u}^\*) = I - (I - \Delta)^{-1} \left\{ \Phi\_{\mathfrak{u}}^{-1} \mathbf{G}\_{\mathfrak{u}}(\mathfrak{u}^\*) + I \right\} = \mathcal{D}\_{\mathfrak{u}}F(\mathfrak{u}^\*).$$

where, <sup>D</sup> <sup>=</sup> diag( <sup>ˆ</sup>*d*1, <sup>ˆ</sup>*d*2, <sup>ˆ</sup>*d*3). Using Proposition 4 and Equation (33), it follows,

$$\begin{cases} \begin{array}{ll} H(\mu\_1) = H(0) > 0, \\ H(\mu\_i) < 0, & 2 \le i \le n\_i, \\ H(\mu\_i) > 0. & i \ge n+1. \end{array} \end{cases} $$

Therefore, 0 is not an eigenvalue of the matrix *<sup>μ</sup>iI* <sup>−</sup> <sup>Φ</sup>−<sup>1</sup> *<sup>u</sup>* (**u**∗)**Gu**(**u**∗) for all *i* ≥ 1, and

$$\sum\_{i \ge 1, H(\mu\_i) < 0} \dim E(\mu\_i) = \sum\_{i=2}^n \dim E(\mu\_i) = \sigma\_{n\nu} \qquad \text{which is odd.}$$

By Proposition 3,

$$\text{index}(F(1;.), \mathfrak{u}^\*) = (-1)^\tau = (-1)^{\sigma\_n} = -1 \tag{35}$$

Similarly,

$$\text{index}(F(0;.), \mathbf{u}^\*) = \left(-1\right)^0 = 1\tag{36}$$

By Theorems 3 and 4, there exists a positive constants *C* such that, for all 0 ≤ *t* ≤ 1, the positive solutions of Equation (34) satisfying <sup>1</sup> *<sup>C</sup>* <sup>≤</sup> *<sup>u</sup>*, *<sup>v</sup>*,*<sup>r</sup>* <sup>≤</sup> *<sup>C</sup>*. Therefore, *<sup>F</sup>*(*t*; **<sup>u</sup>**) <sup>=</sup> 0 on *<sup>∂</sup><sup>B</sup>* for all 0 <sup>≤</sup> *<sup>t</sup>* <sup>≤</sup> 1. By the homotopy invariance of the topological degree,

$$\deg(F(1;.),0,B(\mathbb{C})) = \deg(F(0,.),0,B(\mathbb{C})).\tag{37}$$

However, by our deduction, both equations *F*(1; **u**) = 0 and *F*(0; **u**) = 0 have only the positive solution **u**∗ in *B*(*C*), and hence by Equations (35) and (36)

$$\deg(F(0;.),0,B(\mathbb{C})) = \text{index}(F(0;.),\mathbf{u}^\*) = 1$$

$$\deg(F(1;.),0,B(\mathbb{C})) = \text{index}(F(1;.),\mathbf{u}^\*) = -1$$

which contradicts Equation (37) and we complete the proof.

#### **6. Numerical Simulation**

The numerical simulation was carried out to understand the spatiotemporal dynamics of top predator *r* under the influence of cross diffusion. For this purpose, in this section, we present a detailed investigation of the patterns in the top level predator for different diffusivity coefficients. The system of partial differential equations given in Equation (8) was numerically solved using semi implicit finite difference technique. Forward difference scheme was used for the reaction terms and standard five point explicit finite difference scheme was used for diffusion term. Turing patterns were obtained from the effect of nonlinear diffusion term for the top level predator *r* as in [18–20].

To discretize the system, we considered Taylor's series expansion about the non-trivial equilibrium point (*u*∗, *v*∗,*r*∗) of the cross diffusive term of the third equation describing the dynamics of *r* for the model system in Equation (8). We obtained a system of the following form:

$$\begin{aligned} \frac{\partial u}{\partial t} - d\_1 \Delta u &= u \left( 1 - u - \frac{v}{u + w\_4} \right) \\ \frac{\partial v}{\partial t} - d\_2 \Delta v &= v \left( -w\_5 + \frac{w\_6 u}{u + w\_7} - cv - \frac{r}{v + (w\_8 + w\_9 v)} r + w\_{10} \right) \\ \frac{\partial r}{\partial t} - d\_3 d\_4 r^\* \Delta v - d\_3 (1 + d\_4 r^\*) \Delta r &= r \left( -w\_{11} + \frac{w\_{12} v}{v + (w\_8 + w\_9 v)} r + w\_{10} \right) \\ \frac{\partial u}{\partial n} = \frac{\partial v}{\partial n} = \frac{\partial r}{\partial n} = 0 & \quad (x, y) \in \partial \Omega. \end{aligned} \tag{38}$$

The space step size and time step size were chosen appropriately to ensure the convergence of the scheme.

We used the standard five-point approximation for the two-dimensional Laplacian with the zero-flux boundary conditions. Initially, the entire system was at the stationary state (*u*∗, *v*∗,*r*∗). The perturbation introduced in the initial condition was of the order 5 × <sup>10</sup><sup>−</sup>4, as given in Equation (39):

$$\begin{aligned} u(\mathbf{x}, y) &= u^\* + \epsilon\_1 \sin\left(\frac{2\pi(\mathbf{x} - \mathbf{x}\_0)}{0.2}\right) + \epsilon\_2 \sin\left(\frac{2\pi(y - y\_0)}{0.2}\right) \\ v(\mathbf{x}, y) &= v^\* + \epsilon\_1 \sin\left(\frac{2\pi(\mathbf{x} - \mathbf{x}\_0)}{0.2}\right) + \epsilon\_2 \sin\left(\frac{2\pi(y - y\_0)}{0.2}\right) \\ r(\mathbf{x}, y) &= r^\* + \epsilon\_1 \sin\left(\frac{2\pi(\mathbf{x} - \mathbf{x}\_0)}{0.2}\right) + \epsilon\_2 \sin\left(\frac{2\pi(y - y\_0)}{0.2}\right) \end{aligned} \tag{39}$$

where <sup>1</sup> = <sup>2</sup> = <sup>5</sup> × <sup>10</sup><sup>−</sup>4, *<sup>x</sup>*<sup>0</sup> = *<sup>y</sup>*<sup>0</sup> = 0.1.

The set of parameter values at which the system would yield a locally asymptotically stable solution is:

$$\begin{aligned} w\_4 &= 0.25, & w\_5 &= 0.25, & w\_6 &= 0.8, & w\_7 &= 0.25, & w\_8 &= 0.01, \\ w\_9 &= 0.1, & w\_{10} &= 0.28, & w\_{11} &= 0.25, & w\_{12} &= 0.78, & \varepsilon &= 2. \end{aligned}$$

To analyze the role of cross diffusion on *r*, we considered the following set of parameters:

$$\begin{aligned} w4 &= 17.25, & w5 &= 17.25, & w6 &= 17.25, & w7 &= 17.25, & w8 &= 0.25, \\ w9 &= 0.25, & w10 &= 0.28, & w11 &= 18.26, & w12 &= 3.05, & \varepsilon &= 22. \end{aligned}$$

The parameter values of diffusivity coefficients are presented in Table 1. We performed the simulation on a 50 × 50 grid with spatial step size 0.5 and time step size 0.1. To investigate the role of cross diffusion and self diffusion in the pattern formation of the top predator *r*, we performed simulations for a wide range of self diffusive coefficient *d*<sup>3</sup> and cross diffusivity coefficient *d*4. The different values of self and cross diffusive coefficients used in numerical experiments of top level predator *r* are presented in Table 1. We carried out all the numerical simulation at time level *t* = 10,000 for the model system given in Equation (38).

**Table 1.** Values of diffusivity coefficients *d*1, *d*2, *d*<sup>3</sup> and *d*<sup>4</sup> used in the simulations.


**Figure 1.** No Patterns for top predator of the model system in Equation (38) were obtained at time level *t* = 10,000 in the absence of cross diffusion, i.e., self diffusion and cross diffusion coefficients being *d*<sup>3</sup> = 10, *d*<sup>4</sup> = 0, respectively.

**Figure 2.** (**a**) A mix of hot spot and labyrinth turing patterns obtained at *d*<sup>3</sup> = 1.7 and *d*<sup>4</sup> = 2.3. (**b**) Floral Turing patterns obtained at *d*<sup>3</sup> = 1 and *d*<sup>4</sup> = 5.

**Figure 3.** (**a**) Pentagonal Turing patterns obtained at *d*<sup>3</sup> = 0.8 and *d*<sup>4</sup> = 6. (**b**) Floral Turing patterns obtained at *d*<sup>3</sup> = 0.3 and *d*<sup>4</sup> = 8.

**Figure 4.** (**a**) Another hot spot turing patterns obtained at *d*<sup>3</sup> = 0.5 and *d*<sup>4</sup> = 10. (**b**) A mix of Hot Spot and Labyrinth Turing patterns obtained at *d*<sup>3</sup> = 0.1 and *d*<sup>4</sup> = 15.

**Figure 5.** (**a**) Hot Spot patterns obtained at *d*<sup>3</sup> = 0.05 and *d*<sup>4</sup> = 85. (**b**) Hexagonal Spot Turing patterns obtained at *d*<sup>3</sup> = 0.06 and *d*<sup>4</sup> = 125.

**Figure 6.** (**a**) A mixture of Hot Spot and Labyrinth Turing patterns obtained at *d*<sup>3</sup> = *d*<sup>4</sup> = 0.5. (**b**) A mix of Hot Spot and Labyrinth Turing patterns obtained at *d*<sup>3</sup> = *d*<sup>4</sup> = 1.

**Figure 7.** (**a**) A mixture of hot spot and labyrinth turing patterns obtained at *d*<sup>3</sup> = *d*<sup>4</sup> = 1.5. (**b**) A mixture of hot spot and labyrinth turing patterns obtained at *d*<sup>3</sup> = *d*<sup>4</sup> = 2.

Here, we performed numerical simulation to investigate the role of self and cross diffusion on pattern formation. In Figure 1, we notice that, in the absence of cross diffusion, i.e., *d*<sup>4</sup> = 0, no patterns were formed for the top level predator *r*, of the system in Equation (38). The behavior persists even for the higher values of self diffusivity coefficients *d*3. This numerical experiment shows that in the absence of cross diffusion there is no destabilizing effect on the system even when there is self diffusion in the system.

In our second numerical experiment, we introduced the cross diffusion by gradually increasing the cross diffusivity coefficient *d*<sup>4</sup> and along with that we gradually decreased the self diffusivity coefficient *d*3. The goal was to see the impact of cross diffusivity on the dynamics when there is very less role of self diffusion in the system. To this end, in Figure 2a, we increase *d*<sup>4</sup> from 0 to 2.3 and decrease *d*<sup>3</sup> from 10 to 1.7. An increase in cross diffusion coefficient results in destabilization in the top predator *r* dynamics leading to a mixture of hot spots and labyrinth patterns, as seen in Figure 2a. We further increased *d*<sup>4</sup> to 5 and decreased *d*<sup>3</sup> to 1 to observe floral turing patterns (Figure 2b).

On further increasing the cross diffusivity coefficient at *d*<sup>4</sup> = 6 and taking self diffusion coefficient at *d*<sup>3</sup> = 0.8, we observed pentagonal Turing patterns, as shown in Figure 3a. Further, at *d*<sup>4</sup> = 8 and *d*<sup>3</sup> = 0.3, we obtained Floral Turing Pattern.

Increasing the value of *d*<sup>4</sup> to 10 and fixing *d*<sup>3</sup> = 0.5, we observed hot spothl turing patterns, as given in Figure 4a at *d*<sup>4</sup> = 15. At *d*<sup>3</sup> = 0.1, a mix of hot spot and Labyrinth Turing patterns were observed, as shown in Figure 4b.

As shown in Figure 5a,b, we observed that, for a very high value of cross diffusion coefficients *d*<sup>4</sup> = 85 and *d*<sup>4</sup> = 125, and very low values of *d*<sup>3</sup> = 0.05 and *d*<sup>3</sup> = 0.06, hot spot and Hexagonal Spot Turing patterns are obtained, respectively.

From these experiments, we can conclude that the role of cross diffusion in destabilizing the dynamics is much more significant than self diffusion.

In our third numerical experiment, we started with equal values of both *d*<sup>3</sup> and *d*<sup>4</sup> fixed at 0.5. The hot spot patterns began to appear. Further, when we increased diffusivity coefficient *d*<sup>3</sup> and *d*<sup>4</sup> to 1, a mix of hot spot and Labyrinth Turing Patterns were seen. We further increased both *d*<sup>3</sup> and *d*<sup>4</sup> with the restriction that both self and cross diffusion coefficient must remain the same. We observed that the patterns became more conspicuous, leading to a mix of hot spot and Labyrinth Turing Patterns (see Figures 6a,b and 7a,b).

Through extensive simulations over a wide range of diffusivity coefficients, we observed that cross diffusion coefficient *d*<sup>4</sup> plays an important role to realize the phenomenon of pattern formation. In the absence of cross diffusion, no patterns were formed, even at higher values of *d*3, as shown in Figure 1. Moreover we also determined the Turing patterns in the cases where cross diffusion coefficient equals to the self diffusion coefficient, i.e., *d*<sup>3</sup> = *d*4. Table 1 describes the various self and cross diffusivity coefficients used for obtaining the spatial patterns of the system in Equation (38). Table 1 indicates a negative correlation between self diffusion coefficient and cross diffusion coefficients, which results in pattern formation. The above simulations results are coherent to our analytical results that non-constant positive solutions exist only in the presence of cross diffusion hence the pattern formation.

#### **7. Conclusions**

In this study, we concluded a three-dimensional continuous time ecological system, modeling a tritrophic food chain based on a hybrid type of Holling Type II and Crowley-Martin functional response. We also took into account inhomogeneous spatial distribution of the species involved and dependence of movements of the species on the population density due to the mutual interference between the individuals of the same species and included nonlinear self diffusion and cross diffusion in our discussion. We established that, under certain parametric conditions, the model is unstable only in the presence of cross diffusion. This phenomenon can be regarded as an extension to diffusion driven instability and can be referred to as cross diffusion driven instability. The phenomenon of pattern formation in spatially extended models including cross diffusion as a destabilizing mechanism was studied with varying the self diffusion and cross diffusion coefficients over a wide range, which serves as an extension to Turing pattern formations exhibited by reaction diffusion system. In addition, we analytically proved the existence of inhomogeneous steady states and gave a priori upper and lower bounds for the positive solutions of the system under steady state. We proved that no non-constant positive solutions exist in the presence of only self diffusion subject to certain conditions. Through simulation, we concluded that cross diffusion is necessary for pattern formation of models in which species interactions are governed by Crowley-Martin and Holling Type II functional response. This establishes the fact that Crowley-Martin functional response governed species interaction is stable in spatial temporal aspects when the species diffuse because of their own characteristic behavior but when under the influence of other species the system tends towards instability. Based on our analysis and simulation, we concluded that cross diffusion induces at least a inhomogeneous stationary solution and mathematically explains the segregation phenomenon of ecosystem under the influence of predator dependent functional response.

**Author Contributions:** The authors contributed equally to this work.

**Funding:** The research of the first author was funded by Science and Engineering Research Board (SERB), grant number EMR/2017/005203.

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

#### **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/).

*Article*
