Turing patterns induced by cross-diffusion in a predator-prey system in presence of habitat complexity

Chaos, Solitons and Fractals 91 (2016) 421–429

Contents lists available at ScienceDirect

Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos

Turing patterns induced by cross-diffusion in a predator-prey system in presence of habitat complexity Santu Ghorai, Swarup Poria∗ Department of Applied Mathematics, University of Calcutta, 92 APC Road, Kolkata-700009, India

a r t i c l e

i n f o

Article history: Received 10 March 2016 Revised 1 July 2016 Accepted 8 July 2016

Keywords: Habitat complexity Cross diffusion Two parameter bifurcation Turing pattern formation Control of pattern

a b s t r a c t In this paper, we have investigated the phenomena of Turing pattern formation in a predator-prey model with habitat complexity in presence of cross diffusion. Using the linear stability analysis, the conditions for the existence of stationary pattern and the existence of Hopf bifurcation are obtained. It is shown analytically that the presence of cross diffusion in the system supports the formation of Turing pattern. Two parameter bifurcation analysis are done analytically and corresponding bifurcation diagrams are presented numerically. A series of simulation results are plotted for different biologically meaningful parameter values. Effects of variation of habitat complexity and the predator mortality rate and birth rate of prey on pattern formation are also reported. It is shown that cross-diffusion can lead to a wide variety of spatial and spatiotemporal pattern formation. It is found that the model exhibits spot and stripe pattern, and coexistence of both spot and strip patterns under the zero flux boundary condition. It is observed that cross-diffusion, habitat complexity, birth rate of prey and predator’s mortality rate play a significant role in the pattern formation of a distributed population system of predator-prey type. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Predator-prey dynamics have been continuing as one of the dominant themes in mathematical ecology since the pioneer work of Lotka [1] and Volterra [2]. The main goal of mathematical ecology is to understand the relationship between different species and their living environment. The crucial components of a mathematical model of the predator prey system are the growth function of prey species, the mortality function of predator species and the functional response. Several types of functional responses such as Holling type I -IV [1–5], Beddington–Deangelis [6] etc. have been used for modeling predator prey systems. The existence of habitat complexity in an eco-system reduces the probability of capturing a prey by reducing the searching efficiency of predator and habitat complexity affects the attack coefficient [7] as a result the functional response will be modified. Habitat complexity is the structural complexity of habitats. Habitat complexity can strongly mediate predator-prey interactions, affecting not only total predation rates, but also modifying selectivities for different prey species [8– 11]. Pennings [12] and Grabowski [13] reported that habitat complexity reduces encounter rates of predators with prey. Aquatic habitat becomes structurally complex in presence of submerged

∗

Corresponding author. E-mail address: [email protected] (S. Poria).

http://dx.doi.org/10.1016/j.chaos.2016.07.003 0960-0779/© 2016 Elsevier Ltd. All rights reserved.

vegetation or aquatic weeds. It is observed that structural complexity of the habitat stabilizes the predator-prey interaction between piscivorous perch (predator) and juvenile perch and roach (prey) by reducing predator foraging efficiency. Luckinbill prolonged the coexistence of Paramecium aurelia (prey) and Didinium nasutum (predator) in laboratory system by increasing strength of habitat complexity using methyl cellulose in the Cerophyl medium (nutrient)[14]. Therefore, it is important to incorporate the effect of habitat complexity in predator-prey functional response for theoretical models. Reaction-diffusion (RD) systems have attracted increasing attention from the mathematical biologists in recent years to seek insights into the fascinating patterns that occur in living organisms and in ecological systems. Turing instability constitutes a universal paradigm for the spontaneous generation of spatially organized patterns. Alan Turing [15] proposed a dynamical mechanism, which has been extensively used to explain how Turing patterns are formed and is now known as Turing bifurcation or Turing instability. Turing pointed out that to generate spatial patterns, a reaction diffusion system should contain at least two reactive species that diffuse at very different rates: one slowly diffusing substance and other rapidly diffusing substance. Based on the pioneer work of Turing [15], Segel and Jackson [16] first introduced the reaction diffusion system in ecology. Under some ecological settings, diffusion should be thought of as dispersal of population density and often be considered as a stabilizing process, thus it is the diffusion

422

S. Ghorai, S. Poria / Chaos, Solitons and Fractals 91 (2016) 421–429

results in spatial patterns [17–20]. Existence of non-Turing shown the pattern under the self diffusion was reported [21–23] when the homogeneous state is unstable but all those investigations are based on the consequence of Turing concept. Effects of migration (directed movement) and diffusion (random movement) in a predator-prey system with Holling-Tanner form was reported by Sun et al. [24]. Recently an epidemic model with spatial diffusion is considered by Li [25] to understand the patterns in disease spreading. To reveal the relationship between inducible defenses and herbivore outbreak, a reaction diffusion model was proposed and analyzed by Sun et al. [26]. Periodic solutions in a herbivore-plant system with time delay and spatial diffusion was first reported by Li et al. [27]. Relationship between pattern structures and ecosystems collapse was investigated by Sun et al. [28] formulating a reaction diffusion model. Systems with cross-diffusion are widespread in nature and play an important role, in biophysical and biomedical situations. The escaping-chasing phenomena was described of cross-diffusion by Kerner (1959) [29], and first applied in competitive systems by Shigesada et al. (1979) [30]. An important example of crossdiffusion is chemotaxis, in which a field (e.g., an insect population) diffuses towards regions of high concentrations of another field (e.g., a chemically desirable substance). The escaping velocity of preys may be considered as dispersive velocity of the predators and the chase velocity of predators may be considered as dispersive velocity of prey. The presence of cross-diffusion in reaction process yields a large variety of pattern-forming instabilities. In case of cross-diffusion, the gradient in the concentration of one species induces a flux of another species which is very significant in generating spatial structures. After worked of Kerner [29], Shigesada [30], the study of in ecological models with cross diffusion has attracted attention of biologists [31–35]. The resulting governing system is a reaction-diffusion equation of the linear form can be written in the following form :

∂U = G ( U ) + D∇ 2 U. ∂t Where U is the n dimensional vector with components ui (x, t), (i=1,2, , n) as the species- densities. D is the n × n matrix of the diffusion coefficients, where the diagonal elements are called the self-diffusion coefficient and the off-diagonal elements are called cross-diffusion coefficients, G is the reaction term indicating the interaction between the involved species and ∇ 2 be Laplacian operator. Recently cross-diffusion driven instabilities have gained a considerable attention in population dynamics [32–36], mainly due to their ability to predict some important features in the study of the spatial distribution of species in ecological systems [30,37,38], in epidemic system [39]. Despite recognition of the significance of habitat complexity in community dynamics, current reaction diffusion models of theoretical ecology have rarely considered its effects. In this paper, we consider a mathematical model of predator prey system in presence of habitat complexity [40]. The effects of cross diffusion term together with self diffusion term in the pattern forming instability are investigated both analytically and numerically. The effects of variation of habitat complexity parameter, birth rate of prey, mortality rate of predator and cross diffusion term are investigated. The goal of this paper is to show that cross-diffusion can drive the emergence of Turing pattern and the impact of variation of cross-diffusion coefficients in the spatially inhomogeneous distribution of population density. Our analyze reveal that the crossdiffusion may induce spatial patterns and different cross-diffusion coefficient may lead to transition between spot and stripe patterns. The paper is organized as follows. In Section 2, we introduce a diffusive predator-prey model with habitat complexity in presence of the linear cross diffusion term. In Section 3, the analysis of

the local model is carried out. In Section 4, the bifurcation is done. The possibility of existence of Turing pattern in presence of cross diffusion is shown analytically. The result of Turing pattern formation through extensive numerical simulations are presented in Section 5. In Section 6 the controllability aspect of spatiotemporal pattern are discussed. Finally, conclusions are drawn in Section 7. 2. Mathematical model Habitat complexity exists in every terrestrial or aquatic ecological systems. In our model, we assume the logistic growth of prey linear mortality rate of predator. The existence of habitat complexity reduces the probability of capturing a prey by reducing the searching efficiency of predator and habitat complexity affects the attack coefficient [41]. Now following the Jana and Bairagi [40], Sahoo and Poria [42] and introducing a dimension less parameter c(0 < c < 1) that measures the strength of habitat complexity we obtain the following model,





X A ( 1 − c )X Y dX = RX 1 − − dt K B + ( 1 − c )X dY ξ A ( 1 − c )X Y = − δY. dt B + ( 1 − c )X

(1)

where X(t) and Y(t) stand for the prey and predator densities at time t, respectively. The parameter R is the growth rate of the prey, K is its carrying capacity, ξ is the conversion rate of prey into the predator and δ is the death rate of the predator in the absence of prey. In absence of habitat complexity, i.e. when c = 0 the system (4) becomes the well known Rosenzweig–MacArthur model [43]. Now re-defining

Xnew =

Xold , K

tnew = R told ,

Ynew =

AYold BR

the model (1) can be transformed to the following nondimensional form,

XY dX = X (1 − X ) − dt a + bX dY ξ XY = − δY dt a + bX

(2)

where

a=

1 , 1−c

ξnew =

AK ξold , BR

δnew =

δold R

,

b=

K . B

The new parameter a is the new habitat complexity dependent parameter and notice that the model has habitat complicity if a > 1 and in case of no habitat complexity a = 1. In presence of spatial diffusion of species the model (2) takes the following form,

XY ∂X = X (1 − X ) − + D11 ∇ 2 X ∂t a + bX ∂Y ξ XY = − δY + D22 ∇ 2Y. ∂t a + bX

(3)

where D11 > 0 and D22 > 0 are the diffusion coefficients for prey and predator respectively and ∇ 2 = ∂∂x2 + ∂∂y2 is the Laplacian operator. We now focus on the following predator-prey model considering the presence of linear cross-diffusion in our model :

XY ∂X = X (1 − X ) − + D11 ∇ 2 X + D12 ∇ 2Y ∂t a + bX ∂Y ξ XY = − δY + D21 ∇ 2 X + D22 ∇ 2Y. ∂t a + bX

(4)

The coefficient D12 and D21 are called the cross-diffusion coefficient describe the respective population fluxes of preys and predators resulting from the presence of the other species, respectively,

S. Ghorai, S. Poria / Chaos, Solitons and Fractals 91 (2016) 421–429

i.e, the values of D12 imply that the prey species approaches towards the lower concentration of the predator species. The values of D21 imply predator species tends to diffuse in the direction of higher concentration of the prey species. In many cases, the predator prefers to avoid group defense by a huge number of prey and chooses to catch its prey from a smaller concentration group which unable to sufficiently resist [44]. For this reason, it can be assume that D12 can take any value, but D21 is always positive [45]. In this investigation we assume that the cross diffusion coefficient D12 > 0 and D21 > 0. From principle of thermodynamic aspects of diffusion [46], we assume

D11 D22 − D12 D21 > 0

(5)

which indicates that the flow of the respective densities in the spatial domain depends more strongly on their own density than on the others [45]. We investigate the model (4) under the following positive initial conditions

− → − → − → X ( x , 0 ) > 0, Y ( x , 0 ) > 0, x = (x, y ) ∈ 

(6)

with zero-flux boundary conditions

Now we do linear stability analysis of the equilibrium point E∗ of the system (8). We evaluate the Jacobian matrix A0 of the system at E∗ as

⎛ A0 =⎝

3. Local analysis The non-diffusive model can be written in the following form,

dX = f (X, Y ) dt dY = g(X, Y ) dt

(8)

⎞ X a + bX ⎠ = ξX −δ a + BX E∗ −



a11 a21



a12 . a22 (12)

δ

a11 = ξ (ξ −bδ ) [ξ (b − a ) − bδ (a + b)]; a12 = − ξδ ; a21 =

where,

ξ − aδ − bδ and a22 = 0. The characteristic equation for A0 is λ2 − tr (A0 )λ + 0 = 0, and the roots are given by,



λ1,2 (0 ) =

1 tr (A0 ) ± 2

(tr (A0 ) )2 − 4





0

where,

δ [ξ (b − a ) − bδ (a + b)] and ξ (ξ − bδ ) δ = det (A0 ) = (ξ − aδ − bδ ). ξ

tr (A0 ) =

(7)

where  is the 2-D system size and ν is the outward unit normal vector of the boundary ∂ . The Eq. 7 imply that no external input is imposed from outside i.e, species in  can’t leave out side of the domain  and out side species can not enter in the domain  .

aY (a + bX )2 aξ Y (a + bX )2

1 − 2X −



∂ X ∂Y = = 0, On ∂  × [0, t]. ∂ν ∂ν

423

0

(13)

Notice that at E∗ the stability of the non-diffusive model (8) only depends only the trace of the Jacobian matrix A0 and due to absence of the a22 term the stability nature only depends on the sign of the a11 term. With the condition b ≤ a, the local system (8) globally asymptotically stable. From the Eq. 13 it is clear that if b > a then the non-diffusive model (8) is asymptotically stable. In general, if b > a then the non-diffusive model (8) is stable if δ > δ H . Where δ H is the Hopf bifurcation parameter, defined by the following expression

δH =

ξ b



b−a b+a

 (14)

where

f (X, Y ) = X (1 − X ) −

XY a + bX

and

g(X, Y ) =

ξ XY − δY. a + bX

4. Bifurcation analysis

Lemma 1. All the solutions of system (8) which initiate in R+ are 2 uniformly bounded within the region , where

  Y 1  = (X, Y ) : 0 ≤ X + ≤ 1 + ξ 4δ

(9)

Proof : Let us consider the function Z (t ) = X (t ) + Y ξ(t ) Then after straight forward calculation corresponding to the trajectories of the system (8). We obtain

1 dZ (t ) + δ Z (t ) = X (1 − X ) + δ ≤ + δ dt 4

(10)

Now using the comparison theorem of ordinary differential equation we obtain ∀ T ∈ [0 t]

0≤1+





1 1 − 1+ − Z (T ) e−(t−T ) 4δ 4δ

Hence we have lim sup Z (t ) ≤ 1 + t→∞

1 . 4δ

(11)

Proof is complete.

Equilibria : The model (8) has three equilibrium points. The equilibrium points are E0 ≡ (0, 0), E1 ≡ (1, 0) and E∗ ≡ (X∗ , Y∗ ) where

X∗ =

aδ aξ (ξ − aδ − bδ ) , Y∗ = ξ − bδ (ξ − bδ )2

From the biological point of view we are only interested in the stability behavior of the positive equilibrium point. Note that the non-trivial equilibrium point is positive if ξ > δ (a + b).

In this section we derive cross-diffusion-driven instability conditions and show that these are a generalization of the classical diffusion-driven instability conditions in the absence of crossdiffusion. One of the most revealing results is that, in contrast to the classical reaction-diffusion systems without cross-diffusion, it is no longer necessary to enforce that one of the species diffuse much faster than the other. Furthermore, it is no longer necessary to have an activator-inhibitor mechanism for pattern formation, activator-activator, inhibitor-inhibitor reaction kinetics as well as short-range inhibition and long-range activation all have the potential of forming patterns due to cross-diffusion-driven instability. We linearize the system (4) around E∗ for small space and time dependent fluctuations and expand in Fourier space as − →− → x

− → X( x ,t) ∼

X ∗ eλt ei k

− → Y ( x ,t) ∼

Y ∗ eλt ei k

− →− → x

where k, λ are the wave number vector and frequency, respectively. Substituting the above relations in (4) and neglecting the higher order terms we obtain the characteristic equation,

|Ak − λI| = 0 where Ak = A0 − k2 D and D =



D11 D21

D12 D22



the diffusion matrix.  The roots of the characteristic equation λ2 − tr (Ak )λ + k = 0

424

S. Ghorai, S. Poria / Chaos, Solitons and Fractals 91 (2016) 421–429

are



λ1,2 (k ) =

1 tr (Ak ) ± 2

(tr (Ak ))2 − 4

In view of relation (5) and (22) we have





D11 D22 D12 a21 − D22 a11 > D21 > > 0. D12 −a12

k

where tr (Ak ) = tr (A0 ) − k (D11 + D22 ) 2

and



k

= det (Ak ) = (D11 D22 − D21 D12 )k4 − (D22 a11 + D11 a22 − D21 a12 − D12 a21 )k2 +

 0

.

(15) Turing bifurcation : Turing bifurcation occurs when the equilibrium state that is stable in absence of diffusion but it becomes unstable in presence of diffusion. The matrix Ak is dependent on the wave number k. The system (4) will be stable at the equilibrium point E∗ if tr(Ak ) < 0 and k > 0. The stability of the local system (8) guarantee that at the equilibrium point E∗ , tr(Ak ) < 0, since

tr (Ak ) = tr (A0 ) − k2 (D1 + D2 ) < tr (A0 ) < 0. Thus, the only way for the E∗ to become an unstable point of the diffusive system (4) is that k < 0. The determinant k is the following quadratic polynomial in k2 .  h(k2 ) = k4 − k2 + . (16) 0

Here, = D22 a11 + D11 a22 − D21 a12 − D12 a21 and = D11 D22 − D12 D21 > 0. After a straight forward calculation h(k2 ) can be written in the following form



h ( k2 ) =

√

k2 − √ 2

2



+

 0

−



2 . 4

(17)

(23)

Since the equilibrium point of the local system (8) is stable the conditions (21a) and (21b) hold. Next condition (21c) holds when the criteria (22) satisfied which can be done by choosing the diffusion parameter suitable. Finally the condition (21d) should be satisfied for obtaining Turing patterns. The condition (21d) can be written as

 δ δ ξ + bδ b−a D22 − (ξ − bδ − aδ )D12 + D21 ξ ξ − bδ ξ  δ >2 (ξ − aδ − bδ )(D11 D22 − D12 D21 ). ξ

(24)

Taking δ as a bifurcation parameter the critical parameter value for Turing bifurcation δ T can be obtained by solving the inequation (24). Note that in absence of cross diffusion (i.e, for D12 = 0,D21 = 0) the condition for Turing instability is D22 a11 + D11 a22 >  2 D11 D22 0 . However this condition does not hold due to the absence of a22 term. Hence Turing pattern does not occur in the model (4) in absence of cross diffusion term. It is also important to note that in presence of cross diffusion the diffusion parameters D11 and D22 need not necessary be unequal. This is the one of the most important feature for obtaining Turing instability in presence of cross diffusion.

The minimum value of h occurs when the first square term vanish in the above equation i.e, at

k2min =

5. Numerical simulation results

. 2

(18)

and the corresponding minimum value of h(k2 ) is given by the relation,

h(k2min ) =

 0

−

2 . 4

(19)

Hence the Turing instability condition is h(k2min ) < 0 and the critical wave number is k2cr =



0 . Corresponding the critical wave

length is λT = 2kπ . Again from the equation h(k2 ) = 0, it can be cr easily observed that the range of the wave numbers for which the 2 h(k ) < 0 is given by k1 < k < k2 where

1

k1,2 = ∓ 2 2



− 4 2

 0

.

(20)

The Turing bifurcation conditions in presence of cross diffusion term can be written as the following :

a11 + a22 < 0  0

The non-uniform stationary state of a prey-predator model, which corresponds to the spatio-temporal patterns, cannot be found analytically. Therefore, computer simulations are necessary to understand the behavior of the dynamics of the model. In this section we perform a series of numerical simulations of the model (4).

(21a)

= a11 a22 − a12 a21 > 0

(21b)

5.1. Numerical scheme We simulate our model in the certain bounded spatial domain

 ⊆ R2 . We assume here zero flux boundary conditions. Mathematically, it is equivalent to the Neumann boundary condition. The main reason for choosing such boundary conditions is that we are interested in the self-organization of patterns and zero flux boundary conditions imply no external input. We simulate the model (4) for all (x, y) ∈  with positive initial condition around the point E∗ in particularly we assume the initial condition (6) takes the following form

X (x, y, 0 ) = X ∗ − 1 × ψ1

= D22 a11 + D11 a22 − D21 a12 − D12 a21 > 0

(21c)



= D22 a11 + D11 a22 − D21 a12 − D12 a21 > 2

 0

.

(21d)

Since k2 is a non negative number so for getting the nonnegative wave number from Eq. (18) we must have

D21 >

D12 a21 − D22 a11 . −a12

(22)

Y (x, y, 0 ) = Y ∗ − 2 × ψ2 . Where 1 = 2 × 10−4 , 2 = 6 × 10−4 and ψi (i = 1, 2 ) are random numbers in between −0.5 – 0.5. With zero flux boundary condi∂ [X, Y ] = 0 on ∂  × [0, t] where ∂  be the boundary of the tion ∂ν domain and ν the outward normal to ∂ . Without loss of generality we assume ∂  ∈ C2 . To solve the system (4) we transform the continuous model to a finite dimensional model, i.e., discretized the model in time and

S. Ghorai, S. Poria / Chaos, Solitons and Fractals 91 (2016) 421–429

425

Fig. 1. ( A ) : Bifurcation diagram in a − δ space. Taking b = 1.5, ξ = 0.1, D11 = 1, D12 = 1, D21 = 4 and D22 = 5. (B) : Bifurcation diagram in D21 − δ space. Taking a = 1.1, b = 1.5, ξ = 0.1, D11 = 1, D12 = 1 and D22 = 5.

space. We use here FTCS scheme in the following way

→n ∂− X = ∂ t i, j

− →n − → X i, j − X ni,−1 j

δt − →n − → − → − → − → X i+1, j + X ni−1, j + X ni, j+1 + X ni, j−1 − 4 X ni, j − → 2 n ∇ X = . h2

Where δ t is the time step and h = δ x = δ y is the step length of both space direction. Substituting k ∼ i + j (N + 1 ) ; i, j = 0, 1, · · · , N where i, j are space discretization index we convert the two dimensional vector as an one dimensional vector. Then after rearranging the term with respect to the time label the scheme can be written in the following matrix form,

− →n+1 − → X k = L X nk .

(25)

We choose  = [0 200] × [0 200] time step t = 0.01 and both space step h = 0.5. As we early proved that the model (4) can not produce the Turing pattern in absence of cross diffusion term (D12 = 0, D21 = 0 ). In this present investigation we are only interested the Turing pattern formation for the model (4). We perform a series of numerical simulation in the Turing domain by varying the δ , a and D21 . We fixed the other parameters at b = 1.5, ξ = 0.1, D11 = D12 = 1, D22 = 5. In the series of numerical calculation we used the MATLAB ‘Jet’ color map for presenting the pattern. It is well known that prey and predator always gives the same type of pattern. This means that there is a tendency that the prey stays away from the predators and the predators will get closer to the prey, which reflects the effect of the cross-diffusion. So we only presented the dynamical behavior of prey specie. We present the dynamical behavior of prey when the pattern independent on time. We also demonstrated two examples to show that the type of pattern for prey and predator remains same. 5.2. Effect of variation of δ According to our model the non-dimensional parameter δ is the ratio of the death rate of predator and growth rate prey. Increase of this ratio implies increase of death of predator and decrease of this ratio implies increase of birth rate of prey. Predators death rate and growth rate of prey are the most important parameters in any food chain model. For this reason the dimensionless parameter δ is considered as a control parameter. In this section we first fixed the control parameter D21 = 4 and draw the bifurcation diagram in a − δ space, which is presented in the Fig. 1 A. Please note that

from the biological point of view the bifurcation diagram is only 1 valid for δ < 15+10 a , which is coming from the positivity condition of equilibrium point (ξ > δ (a + b)). Now from the bifurcation diagram we fixed the parameter value a = 1.1. Hence according to the choice of parameter set the Hopf bifurcation parameter δ H ∼ 0.010256. Since in this case b > a the local dynamics (8) is stable if δ > δ H ∼ 0.010256. To investigate the patterns for various choice of δ we disturb the parameter δ around the δ H and δ T , where the δ T is the appropriate root of δ which can be determined from the expression (24). Now we choose δ around the point δ H and move towards δ T by proper choices. We plotted the behavior of k and the corresponding dispersion relation, which is shown in Fig. 2 B for some selected value of δ . Fig. 2 strongly recommended the existence of Turing region for the model (4). We follow the exact same procedure for investigating the Turing pattern formation for the other cases. Choosing δ = 0.011 we find that the prey system (4) has cold spot pattern which is presented in Fig. 3 A. For choice of δ = 0.013 we found the coexistence of cold spots and stripes, which is plotted in Fig. 3 B. For δ = 0.015 we found the almost cold stripes pattern with some spots, which is presented in Fig. 3 C. On the other hand the hot stripe appears at δ = 0.017, which is depicted in Fig. 3 D. Further increasing of δ = 0.018 we found mixture of hot spots and hot stripes , which is depicted in Fig. 3 E. For the choice of δ = 0.02 the system (4) rise to the Hot spot patterns, which is presented in the Fig. 3 F. Next we present some Turing pattern to investigate the effects of the habitat complexity parameter a. 5.3. Effect of variation of a 1 The parameter a = 1−c measure the complexity of the habitat in our model. In absence of habitat complexity a = 1(c = 0 ), the model is nothing but the well known Rosenzweig–MacArthur model. Habitat complexity is another most important biological parameter which is also choosen as control parameter. In the previous simulation, we concentrated on the parameter δ for fixing the other two control parameter a and D21 and of the other parameters remain fixed. We constructed the bifurcation diagram in a, δ space by fixing the parameter D21 = 4. From the bifurcation diagram 1 A we can choose the appropriate value of a for fixed value of δ in Turing region. Now we investigate the pattern formation for different choice of the parameter a. Let us fix δ = 0.015. Before investigating the effect of the complexity parameter we first look at the dynamical behavior of the

426

S. Ghorai, S. Poria / Chaos, Solitons and Fractals 91 (2016) 421–429

Fig. 2. (A) : Plot of k against k2 for different values of δ . Taking a = 1.1, b = 1.5, ξ = 0.1, D11 = 1, D12 = 1, D21 = 4 and D22 = 5. (B) : The corresponding dispersion relation.

Fig. 3. Evaluation of the prey of the model (4) after the time t = 50 0 0 for a = 1.1 and D21 = 4. ( A ) : δ = 0.011, ( B ) : δ = 0.013, ( C ) : δ = 0.015, ( D ) : δ = 0.017, ( E ) : δ = 0.018 and (F) : δ = 0.02. The other parameters are fixed at b = 1.5, ξ = 0.1, D11 = 1, D12 = 1 and D22 = 5.

system having no habitat complexity which can be found by setting the parameter a = 1 (c = 0 ). We perform the numerical simulation and found that the prey population gives the mixture of clod spot and cold stripes patterns, which is presented in Fig. 4 A. Now if we go back to the previous result shown in Fig. 3 C for a = 1.1 and compare the Fig. 3 C with the Fig. 4 A. Then it is observable that the system (4) having no habitat complexity (a=1) end up with the mixture of cold spot and clod stripe pattern, but due to small (c ∼ 0.09091) present of habitat complexity the pattern turns in to almost cold stripe pattern. Notice that for small changes of the complexity parameter there is an observable change in the density scale. Now set a = 1.15 we found pure stripe which is presented in Fig. 4 B. Further increase of a = 1.2 we get mixture of spots and hot stripes which is depicted in Fig. 4 C. 5.4. Effect of variation of D21 Next we investigate the effect of variation of cross diffusion parameter D21 . On pattern formation fixed a = 1.1 and draw the bifurcation diagram in D21 − δ space which is depicted in Fig. 1 B. In

this section we choose two values of δ which are taken from the Turing region of bifurcation diagram presented in Fig. 1 B. Fixing δ = 0.015 and taking D21 = 3.7. We found stripe pattern which is quite similar to our previous result depicted in Fig. 3 C. Now if we increase D21 and set D21 = 4.5 we found the mixture of cold spots and cold stripes. A Further increasing D21 to 4.8 we found almost similar type of mixture pattern with different quality. Please notice that for all the above three choices of D21 the density scale has observable changes. Hence, due to the change of reaction rate the different variety of mixture patterns are observe for the system (4). One can investigate the different kind of spots and stripes etc by varying the other parameter also. For another evidence of variation of pattern due to the cross diffusion term we consider one more example. Let us change the δ and set δ = 0.02 and the other parameter a = 1.1 fixed as before. Now choose D21 = 4.5 after numerical simulation we found spot pattern as like in Fig. 3 E but notice that there is a more spots in Fig. 6 A when compare with Fig. 3 E. Now setting D21 = 4.8 and 4.9, we found different variety of the mixture of spot and stripe patterns, which is depicted in Figs. 6 B, C. From the Fig. 6 it is clear that

S. Ghorai, S. Poria / Chaos, Solitons and Fractals 91 (2016) 421–429

427

Fig. 4. Evaluation of the prey of the model (4) after the time t = 50 0 0 for δ = 0.015 and D21 = 4. ( A ) : a = 1, ( B ) : a = 1.15, ( C ) : a = 1.2. The other parameters are fixed at b = 1.5, ξ = 0.1, D11 = 1, D12 = 1 and D22 = 5.

Fig. 5. Evaluation of the prey of the model (4) after the time t = 50 0 0 for δ = 0.015 and a = 1.1. ( A ) : D21 = 3.7, ( B ) : D21 = 4.5, ( C ) : D21 = 4.8. The other parameters are fixed at b = 1.5, ξ = 0.1, D11 = 1, D12 = 1 and D22 = 5 .

Fig. 6. Evaluation of the prey (upper row) and predator (lower row) of the model (4) for a = 1.1 and δ = 0.02. ( A, D ) : D21 = 4.5 at t=5,0 0 0, ( B, E ) : D21 = 4.8 at t=3,0 0 0 ( C, F ) : D21 = 4.9 at t=1,0 0 0. The other parameters are fixed at b = 1.5, ξ = 0.1, D11 = 1, D12 = 1 and D22 = 5.

for the choice of different cross diffusion parameter the dynamical behavior of prey and predator are always of the same type. This means that there is a tendency that the prey stays away from the predators and the predators will get closer to the prey. Now we give one more example for which we strongly establish the fact that change of the local parameter for the model (4) under the fixed cross diffusion part the dynamical behavior of prey and

predator remain unchanged. Let us fixed a = 1.2 and D21 = 4.8 and take δ as the control parameter. Setting δ = 0.008, 0.011 and 0.02 and simulated the model (4) for a long time. It is found that in different instant of time the dynamics of prey and predator always same type, which is presented in Fig. 7.

428

S. Ghorai, S. Poria / Chaos, Solitons and Fractals 91 (2016) 421–429

Fig. 7. Evaluation of the prey (upper row) and predator (lower row) of the model (4) after the time t = 20 0 0 for a = 1.2 and D21 = 4.8. ( A, D ) : δ = 0.008, ( B, E ) : δ = 0.011, ( C, F ) : δ = 0.02. The other parameters are fixed at b = 1.5, ξ = 0.1, D11 = 1, D12 = 1 and D22 = 5.

Fig. 8. Surface plot of the prey of the model (4) after the time t = 10 0 0 for a = 1.1 and D21 = 4. ( A ) : δ = 0.006, ( B ) : δ = 0.011, ( C ) : δ = 0.03. The other parameters are fixed at b = 1.5, ξ = 0.1, D11 = 1, D12 = 1 and D22 = 5.

6. Controllability Controlling an ecological system is an challenging task to mathematical biologists for many years. By appropriate setting the biological parameter a model can be stabilized forever. Consider Fig. 3 we found different types of spots, stripes or their mixture by setting different values of δ . Now if we set δ > δ T ( ∼ 0.0202505) then we can not found the pattern since the simulation region fall in the stable region. This also mean that the with respect to parameter set a = 1.1, b = 1.5, ξ = 0.1, D11 = 1, D12 = 1, D21 = 4 and D22 = 5, the appropriate choice of δ (0.0202505 < δ < 0.0384615) can control the system (4) to fully stable system. For the strong evidence of control we present the 3-dimensional view in 8 for the each prey node (X(i, j)) after time t = 1, 0 0 0 for the model (4). From the bifurcation diagram Fig. 1 A we set δ = 0.06 which corresponds to give the Hopf-Turing region. For choice of δ = 0.006 we found that the each point of the prey situated in unsymmetrical order due to both Turing as well as Hopf oscillation, which is present in Fig. 8 A. For choice δ = 0.011 which is in the Turing region, the each point of prey are situated in symmetrical order, which is depicted in Fig. 8 B. For the stable region we choose δ = 0.03 which is shown in the Fig. 8 C. The above discussion reflect that the spatiotemporal dynamics can be control though the parameter δ when the other parameter

are fixed. Again from the bifurcation diagram Fig. 1 we can control the patterns through the variation of a and D21 . In general the stability region can be defined when the condition (21a) is hold and (21d) is violated together with the condition > δ (a + b) (cf Fig. 1). Also note that the condition (21a) always hold if a > b and the condition (21d) always violated if (21c) is violated. 7. Conclusion We have investigated the effects of habitat complexity on the spatio-temporal pattern of a predator-prey model in presence of linear cross diffusion. It is shown that the said model (4) does not produce the Turing pattern in presence of self diffusion terms only (in absence of cross diffusion term). The effects of variation of the three control parameters : the mortality rate of predator, the cross diffusion term and habitat complexity on the pattern formation are reported. Habitat complexity exists in almost every terrestrial or aquatic ecological systems and the effects of habitat complexity 1 can be observed varying the parameter a(= 1−c ). Both the effects of variation of predator’s death rate or predator migration rate on the predator prey system dynamics can be incorporated varying the death rate of predator δ only. Lower the death rate of predator may cause higher death rate of prey via feeding. The dimensionless parameter δ is the ratio between death rate of predator and

S. Ghorai, S. Poria / Chaos, Solitons and Fractals 91 (2016) 421–429

growth rate of prey. Biologically low value of δ means that either growth rate of prey is high or death rate of predator is low. Moreover, the feeding rate of predators may be the determining factor in producing different type of spatial patterns which can be visualize through the variation of predators death rate. We have done two parameter bifurcation analysis of the model in the a − δ parameter space and corresponding figure is plotted in Fig. 1 A. The bifurcation analysis is also done in the D21 − δ parameter space and the corresponding bifurcation diagram is presented in Fig. 1 B. The bifurcation diagrams in Fig. 1 help us to choose the parameters for obtaining suitable patterns as well as for controlling the spatiotemporal patterns. It is clear from the Fig. 1 A that varying the parameter δ only the transition from pure Hopf → Hopf-Turing → Turing → stable patterns are possible. Similar kind of transitions are possible in case of other control parameters also. The effect of variation of δ in the Turing region are presented in Fig. 3. It is observed that transition from cold spots → mixture of cold spots and cold stripes → cold stripes → hot stripes → mixture of hot spots and hot stripes → hot spots are possible in the Turing region with the variation of δ only keeping all other parameters fixed. Our results show that the death rate of predator or growth rate of prey can play a vital role in the spatial predator-prey model (4). Similarly the habitat complexity parameter a (cf Fig. 4) also plays an important role for making the several kind of patterns as well as to control spatiotemporal patterns of the model (4). The cross diffusion parameter D21 plays the crucial role for making the Tuing pattern for the said model ( cf Figs. 5 and 4). The cross diffusion parameter D21 can also have important role for large variation of the Turing patterns (cf Figs. 6 and 7). From the Figs. 6 and 7 we found that the dynamical behavior of prey and predator are always similar type. Choosing the appropriate parameters form the bifurcation diagrams Fig. 1 we can control the spatio-temporal patterns ( cf Fig. 8). One can investigate the effects of variation of other parameters such as b, ξ , D11 , D12 and D22 in the pattern formation for further investigation. Of course, the methods and the results in the present paper can be applied to the pattern formation in the diffusive predator-prey model or other reaction-diffusion systems. We have seen that addition of linear cross diffusion term enables the models to produce much richer dynamics than would otherwise be possible. We hope that our work could be useful to study the effects of cross-diffusion on the pattern formation in a complex population dynamics. Acknowledgments We are thankful to reviewers for their valuable comments and suggestions. References [1] Lotka AJ. Elements of physical biology. Baltimore, Maryland USA: Williams & Wilkins; 1925. [2] Volterra V. Variazione e fluttuazini del numero d’individui in specie animali conviventi. Mem R Accad Naz dei Lincei 1926;2:31–113. [3] Holling CS. The components of predation as revealed by a study of small mammal predation of the european pine sawfiy. Canadian Entomologist 1959;91:293–320. [4] Nagano S, Maeda Y. Phase transitions in predator-prey systems. Phys Rev E 2012;85:011915. [5] Berryman AA. The origins and evolution of predator-prey theory. Ecology 1992:1530–5. [6] Beddington JR. Mutual interference between parasites or predators and its effect on searching efficiency. J Anim Ecol 1975;44(1):331–40. [7] Winfield IJ. The influence of simulated aquatic macrophytes on the zooplankton consumption rate of juvenile roach, rutilus rutilus, rudd, scardinius erythrophthalmus, and perch, perca fluviatilis. J Fish Biol 1986;29(sA):37–48. [8] Anderson O. Optimal foraging by largemouth bass in structured environments. Ecology 1984;65:851–61. [9] Ryer CH. Pipefish foraging: effects of fish size, prey size and altered habitat complexity. Marine Ecol Progr Ser 1988;48(1):37–45.

429

[10] Bell S, McCoy E, Mushinsky H. Habitat structure: the physical arrangement of objects in space. Springer Science & Business Media; 2012. [11] Lassau S, Hochuli D. Effects of habitat complexity on ant assemblages. Ecography 2004;27(2):157–64. [12] Pennings SC. Predator-prey interactions in opisthobranch gastropods: effects of prey body size and habitat complexity. Marine Ecol-Progress Series 1990;62:95–101. [13] Grabowski JH. Habitat complexity disrupts predator-prey interactions but not the trophic cascade on oyster reefs. Ecology 20 04;85(4):995–10 04. [14] Luckinbill LS. Coexistence in laboratory populations of paramecium aurelia and its predator didinium nasutum. Ecology 1973;54(6):1320–7. [15] Turing AM. On the chemical basis of morphogenesis. Philos Trans Roy Soc London B 1952;237:37–72. [16] Segel L, Jackson JL. Dissipative structure: an explanation and an ecological example. J Theor Biol 1972;37(3):545–59. [17] Wang W, Liu Q, Jin Z. Spatiotemporal complexity of a ratio-dependent predator-prey system. Phys Rev E 2007;75(5):051913. [18] Abdusalam H, Fahmy E. Cross-diffusional effect in a telegraph reaction diffusion lotka-volterra two competitive system. Chaos, Solitons Fractals 2003;18(2):259–66. [19] Tang X, Song Y. Bifurcation analysis and turing instability in a diffusive predator-prey model with herd behavior and hyperbolic mortality. Chaos, Solitons Fractals 2015;81:303–14. [20] Yuan S, Xu C, Zhang T. Spatial dynamics in a predator-prey model with herd behavior. Chaos 2013;23(3):033102. [21] Medvinsky A, Petrovskii S, Tikhonova I, Malchow H, Li B. Spatiotemporal complexity of plankton and fish dynamics. SIAM Rev 2002;44(3):311–70. [22] Hu G, Li X, Wang Y. Pattern formation and spatiotemporal chaos in a reaction-diffusion predator-prey system. Nonlinear Dynam 2015;81(1):265–75. [23] Ghorai S, Poria S. Pattern formation and control of spatiotemporal chaos in a reaction diffusion prey-predator system supplying additional food. Chaos, Solitons Fractals 2016;85:57–67. [24] Sun GQ, Zhang J, Song LP, Li BL. Pattern formation of a spatial predator-prey system. Appl Math Comput 2012;218(22):11151–62. [25] Li L. Patch invasion in a spatial epidemic model. Appl Math Comput 2015;258:342–9. [26] Sun GQ, Wang S, Ren Q, Jin Z, Wu Y. Effects of time delay and space on herbivore dynamics: linking inducible defenses of plants to herbivore outbreak. Scientific Reports, 2015; 2015. 5 [27] Li L, Jin Z, Li J. Periodic solutions in a herbivore-plant system with time delay and spatial diffusion. Appl Math Model 2016;40(7):4765–77. [28] Sun G, Wu Z, Wang Z, Jin Z. Influence of isolation degree of spatial patterns on persistence of populations. Nonlinear Dynam 2016;83(1–2):811–19. [29] Kerner E. Further considerations on the statistical mechanics of biological associations. Bulletin Math Biophys 1959;21(2):217–55. [30] Shigesada N, Kawasaki K, Teramoto E. Spatial segregation of interacting species. J Theor Biol 1979;79(1):83–99. [31] Tang X, Song Y. Cross-diffusion induced spatiotemporal patterns in a predator-prey model with herd behavior. Nonlinear Anal 2015;24:36–49. [32] Vanag VK, Epstein I. Cross-diffusion and pattern formation in reaction-diffusion systems. Phys Chem Ch Phys 2009;11(6):897–912. [33] Wang W, Lin Y, Zhang L, Rao F, Tan Y. Complex patterns in a predator-prey model with self and cross-diffusion. Commun Nonlinear Sci Numer Simul 2011;16(4):2006–15. [34] Sun G, Jin Z, Liu Q, Li L. Pattern formation induced by cross-diffusion in a predator-prey system. Chinese Phys B 2008;17(11):3936–41. [35] Gambino G, Lombardo M, Sammartino M. Turing instability and traveling fronts for a nonlinear reaction-diffusion system with cross-diffusion. Math Comput Simulat 2012;82(6):1112–32. [36] Sun G, Jin Z, Zhao Y, Liu Q, Li L. Spatial pattern in a predator-prey system with both self-and cross-diffusion. Int J Mod Phys C 2009;20(01):71–84. [37] Murray JD. Mathematical biology II: spatial models and biomedical applications. third ed. New York: Springer-Verlag; 2003. [38] Sun G, Jin Z, Li L, Haque M, Li B. Spatial patterns of a predator-prey model with cross diffusion. Nonlinear Dynam 2012;69(1–2):1631–8. [39] Sun G, Jin Z, QX L, Li L. Spatial pattern in an epidemic system with cross-diffusion of the susceptible. J Biol Syst 2009;17(1):1–12. [40] Jana D, Bairagi N. Habitat complexity, dispersal and metapopulations: macroscopic study of a predator-prey system. Ecol Complexity 2014;17:131–9. [41] Kot M. Elements of mathematical ecology. Cambridge: Cambridge University Press; 2001. [42] Sahoo B, Poria S. Effects of additional food on an ecoepidemic model with time delay on infection. Appl Math Comput 2014;245:17–35. [43] Rosenzweig M, MacArthur R. Graphical representation and stability conditions of predator-prey interaction. Am Nat 1963;97:209–23. [44] Dubey B, Das B, Hussain J. A predator-prey interaction model with self and cross-diffusion. Ecol Model 2001;141(1):67–76. [45] Baek H. Bifurcation analysis of a predator-prey system with self and cross-diffusion and constant harvesting rate. Electron J Qual Theor Diff Equ 2014;29:1–14. [46] Vitagliano V. Some phenomenological and thermodynamic aspects of diffusion in multicomponent systems. Pure Appl Chem 1991;63(10):1441–8.