Table Of Content

Cross-diffusion induced Turing instability in two-prey ✩ one-predator system Yahui Chena, Canrong Tianb, Zhi Linga aSchool of Mathematical Science, Yangzhou University, Yangzhou 225002, P.R. China 5 1 bDepartment of Basic Science, Yancheng Institute of Technology, 0 Yancheng 224003, P. R. China. 2 n a J 3 2 Abstract ] P In this paper, we study a strongly coupled two-prey one-predator system. We first prove A . the unique positive equilibrium solution is globally asymptotically stable for the cor- h t a responding kinetic system (the system without diffusion) and remains locally linearly m [ stable for the reaction-diffusion system without cross-diffusion, hence it does not belong 1 to the classical Turing instability scheme. Moreover we prove that the positive equilib- v 8 rium solution is globally asymptotically stable for the reaction-diffusion system without 0 7 cross-diffusion. But it becomes linear unstable only when cross-diffusion also plays a role 5 0 . in the reaction-diffusion system, thus it is a cross-diffusion induced instability. Finally, 1 0 the corresponding numerical simulations are also demonstrated and we obtain the spatial 5 1 patterns. : v i Keywords: Prey-predator system, Cross-diffusion, Turing instability. X r 2008 MSC: 35K60, 35R35 a ✩ The work supported by PRC grant NSFC (61472343,11201406)and NSF of YZU (2014CXJ002). Email address: [email protected](Zhi Linga) Preprint submitted to Physic A January 26, 2015 1. Introduction In 2009, Elettreby considered the following prey-predator model [3] u′ = au (1−u )−u u := u f (u ,u ), 1 1 1 1 3 1 1 1 3   u′2 = bu2(1−u2)−u2u3 := u2f2(u2,u3), (1.1)  u′ = −cu2 +(du +eu )u := u f (u ,u ,u ), 3 3 1 2 3 3 3 1 2 3    where u ,u and uare the population densities of three species. This system models the 1 2 3 dynamic of two-prey one-predator ecosystem, i.e. the third species preys on the second and the first one. In the absence of any predation, each term of preys grows logistically. The effect of the predation is to reduce the prey growth rate. In the absence of any prey for sustenance, the predator’s death rate results in inverse decay, which is the term −cu2. 3 The prey’s contribution to growth rate of the predators are respectively du u and eu u . 1 3 2 3 They studyed the global stability and persistence of the model. However, in reality, individual organisms are distributed in space. We can use the reaction-diffusion equations to establish spatio-temporal dynamical system which can modelthepursuit-evasionphenomenon(predatorspursuing preyandpreyescaping predators) in the prey-predator system. Therefore, in present paper we further investigate the following reaction-diffusion model with cross-diffusion: u −∆[(k +k u )u ] = au (1−u )−u u , in Ω×(0,∞), 1t 11 13 3 1 1 1 1 3   u −∆[(k +k u )u ] = bu (1−u )−u u , in Ω×(0,∞),  2t 22 23 3 2 2 2 2 3    u −∆[(k u +k u +k )u ] = −cu2 +(du +eu )u , in Ω×(0,∞), (1.2)  3t 31 1 32 2 33 3 3 1 2 3  ∂u1 = ∂u2 = ∂u3 = 0, on ∂Ω×(0,∞), ∂η ∂η ∂η     u1(x,0) = u10(x),u2(x,0) = u20(x),u3(x,0) = u30(x), in Ω,    where Ω is a bounded domain in RN with smooth boundary ∂Ω. η is the unit outward normal to ∂Ω. The homogeneous Neumann boundary condition indicates that there is zero population flux across the boundary. The parameters a,b,c,d,e and k (1 ≤ i, ij j ≤ 3) are all positive constants. k is the diffusion rate of i-th species. This diffusion ii term represent simple Brownian type motion of particle dispersal. k (i 6= j) is the cross- ij diffusion rate of i-th species. It is necessary to note that the cross-diffusion coefficient 2 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. On the contrary, the negative cross-diffusion coefficient denotes the population flux of one species in the direction of higher concentration of another species. Here the cross-diffusion term presents the tendency of predators to avoid the group defense by a large number of prey, i.e. the predator diffuses in the direction of lower concentration of the prey species. More biological background can be found in [1, 14, 17]. As we know, the problem of cross-diffusion was proposed first by Kerner [9] and first applied to competitive population systems by Shigesada et al. [20]. Since then, the role of cross-diffusion in the models of many physical, chemical and biological processes has been extensively studied. In the field of population dynamics some models of multispecies population are described by reaction-diffusion systems. Jorne [8] examined the effect of cross diffusion on the diffusive Lotka-Volterra system. They found that the cross-diffusion may give rise to instability in the system, although this situation seems quite rare from an ecological point of view. Gurtin [6] developed some mathematical models for population dynamics with the inclusion of cross-diffusion as well as self-diffusion and showed that the effect ofcross-diffusion maygiverisetothesegregationoftwo species. Someconditionsfor the existence of global solutions have been given by several authors, for example, Deuring [2], Kim [10], Pozio and Tesei [19], Yamada [25]. Moreover, due to a most interesting qualitative feature: pattern formation induced by cross-diffusion effect there are some works on the diffusion driven instability (Turing instability [22]) and the existence of a non-constant stationary solution, please refer to [11, 12, 13, 15, 16, 18, 23] and the references cited therein. The main purpose of this paper is to study the Turing instability which is driven solely from the effect of cross-diffusion by using mathematical analysis and numerical simulations. The rest of this paper is organized as follows. In section 2 we show that the unique positive equilibrium of the ODE system (1.1) is globally asymptotically stable. In the section 3 we show that the positive equilibrium remains linearly stable in the presence of self-diffusion. It becomes linearly unstable with theinclusion ofsome appropriatecross- 3 diffusion influences. The Turing instability occurs only when the cross-diffusion rates k 23 and k are large. The resulting patterns are computed by a numerical method and also 32 we devoted to some conclusions in section 4. 2. Stability of the positive equilibrium solution of the ODE system In this section, we consider the stability of the positive equilibrium solution of the system (1.1). It is easy to know that if abc > max{e(b−a),d(a−b)} (2.3) the ODE system (1.1) has a unique positive equilibrium u¯ = (u¯ ,u¯ ,u¯ ) which is given by 1 2 3 abc+ae−be abc+bd−ad ab(d+e) u¯ = , u¯ = , u¯ = . (2.4) 1 2 3 abc+bd+ae abc+bd+ae abc+bd+ae We have the following result: Theorem 2.1. The unique positive equilibrium u¯ is globally asymptotically stable for the ODE system (1.1). Proof. In order to prove the theorem, we need construct a Lyapunov function for the system (1.1). u u u V(u(t)) = d(u −u¯ −u¯ ln 1)+e(u −u¯ −u¯ ln 2)+(u −u¯ −u¯ ln 3). (2.5) 1 1 1 2 2 2 3 3 3 u¯ u¯ u¯ 1 2 3 Then V(u¯) = 0 and V(u) > 0 if u 6= u¯. By using (1.1), we compute dV u¯ u¯ u¯ = d(1− 1)u′ +e(1− 2)u′ +(1− 3)u′ dt u 1 u 2 u 3 1 2 3 = d(u −u¯ )[−a(u −u¯ )−(u −u¯ )]+e(u −u¯ )[−b(u −u¯ )−(u −u¯ )] 1 1 1 1 3 3 2 2 2 2 3 3 +(u −u¯ )[−c(u −u¯ )+d(u −u¯ )+e(u −u¯ )] 3 3 3 3 1 1 2 2 = −ad(u −u¯ )2 −be(u −u¯ )2 −c(u −u¯ )2 < 0 1 1 2 2 3 3 for all u 6= u¯. By the Lyapunov-LaSalle invariance principle [7], u¯ given by (2.4) is globally asymptotically stable for the kinetic system (1.1). 4 Theorem 2.2. The unique positive equilibrium u¯ is globally asymptotically stable for the reaction-diffusion system (1.2) without cross-diffusion, i.e. k = 0 for i 6= j. ij Proof. To study the global behavior of system (1.2), we introduce the following Lyapunov functional W(t) = V(u(x,t))dx (2.6) Z Ω where V(u(x,t)) is given by (2.5). By direct computation, we have dW ∂u = grad V · dx u dt Z ∂t Ω u¯ u¯ u¯ 1 2 3 = d(1− ),e(1− ),(1− ) ·(k ∆u +u f ,k ∆u +u f , 11 1 1 1 22 2 2 2 Z (cid:18) u u u (cid:19) Ω 1 2 3 k ∆u +u f )dx 33 3 3 3 u¯ u¯ 1 2 = d k (1− )∆u dx+ e k (1− )∆u dx 11 1 22 2 Z (cid:18) u (cid:19) Z (cid:18) u (cid:19) Ω 1 Ω 2 u¯ dV 3 + k (1− )∆u dx+ dx. 33 3 Z (cid:18) u (cid:19) Z dt Ω 3 Ω From Green’s identity, it follows that u¯ u¯ ∂u u i i i i k (1− )∆u dx = k (1− ) dS − k ∇ (1− )·∇ u dx ii i ii ii x x i Z (cid:18) u (cid:19) Z u ∂n Z u Ω i ∂Ω i Ω i = − k u¯ u−2|∇ u |2dx ≤ 0. ii i i x i Z Ω Since dV ≤ 0, dV ≤ 0. Thus, dW < 0 for all u 6= u¯. By the Lyapunov-LaSalle dt Ω dt dt R invarianceprinciple[7], u¯ isgloballyasymptotically stableforthereaction-diffusionsystem (1.2) without cross-diffusion. 3. Effects of cross-diffusion on Turing instability For simplicity, we denote (k +k u )u au (1−u )−u u 11 13 3 1 1 1 1 3     K(u) = (k +k u )u , F(u) = bu (1−u )−u u . 22 23 3 2 2 2 2 3      (k31u1 +k32u2 +k33)u3   −cu23 +(du1+eu2)u3      5 Then the reaction-diffusion system (1.2) can be rewritten in matrix notation as: ∂u −∆K(u) = F(u) in Ω×(0,∞),  ∂t  ∂u = 0 on Ω×(0,∞), (3.1)  ∂η u(x,0) = (u (x),u (x),u (x))T in Ω.  10 20 30     Linearizingthereaction-diffusionsystem(3.1)aboutthepositiveequilibriumu¯ = (u¯ ,u¯ ,u¯ ), 1 2 3 we have ∂Ψ −Ku(u¯)∆Ψ = Gu(u¯)Ψ (3.2) ∂t where Ψ = (ψ ,ψ ,ψ )T and 1 2 3 k +k u¯ 0 k u¯ 11 13 3 13 1   Ku(u) = 0 k +k u¯ k u¯ , 22 23 3 23 2      k31u¯3 k32u¯3 k33 +k31u¯1 +k32u¯2    −au¯ 0 −u¯ 1 1   Gu(u) = 0 −bu¯ −u¯ . 2 2      du¯3 eu¯3 −cu¯3    Let 0 = µ < µ < µ < · · · be the eigenvalues of the operator −∆ on Ω with the 1 2 3 homogeneous Neumann boundary condition, and E(µ ) be the eigenspace corresponding i to µ in C2(Ω). Let X = {u ∈ [C1(Ω¯)]3| ∂u = 0 on ∂Ω}, {φ } be an i ∂η ij j=1,2,...,dimE(µi) orthonormal basis of E(µ ), and X = {cφ | c ∈ R3}. Then i ij ij ∞ dimE(µi) X = X and X = X . i i ij Mi=1 Mj=1 For each i ≥ 1, Xi is invariant under the operator Ku(u¯)∆+Gu(u¯). Then problem (3.2) has a non-trivial solution of the form Ψ = cφexp(λt) if and only if (λ,c) is an eigenpair for the matrix −µiKu(u¯)+Gu(u¯), where c is a constant vector. Then the equilibrium u¯ is unstable if at least one eigenvalue λ has a positive real part for some µ . i 6 The characteristic polynomial of −µiKu(u¯)+Gu(u¯) is given by ρ (λ) = λ3 +A λ2 +A λ+A , (3.3) i 2i 1i 0i where A = (k +k +k +k u¯ +k u¯ +k u¯ +k u¯ )u +au¯ +bu¯ +cu¯ (3.4) 2i 11 22 33 13 3 23 3 31 1 32 2 i 1 2 3 A = [(k +k u¯ )(k +k u¯ +k +k u¯ +k u¯ )+(k +k u¯ ) 1i 11 13 3 22 23 3 33 13 1 32 2 22 23 3 (k +k u¯ +k u¯ )−k k u¯ u¯ −k k u¯ u¯ ]µ2 33 31 1 32 2 23 32 2 3 13 31 1 3 i +[(bu¯ +cu¯ )(k +k u¯ )+au¯ (k +k +k u¯ +k u¯ +k u¯ ) 2 3 11 13 3 1 22 33 23 3 31 1 32 2 +bu¯ (k +k u¯ +k u¯ )+cu¯ (k +k u¯ )+u¯ u¯ (k e−k ) 2 33 31 1 32 2 3 22 23 3 2 3 23 32 +u¯ u¯ (k d−k )]µ +au¯ (bu¯ +cu¯ )+bcu¯ u¯ +eu¯ u¯ +du¯ u¯ (3.5) 1 3 13 31 i 1 2 3 2 3 2 3 1 3 A = [(k k +k k u¯ +k k u¯ +k u¯ u¯ +k k u¯ +k k u¯ u¯ )k 0i 33 11 33 13 3 31 11 1 31 1 3 32 11 2 32 13 2 3 22 +(k k +k k u¯ +k k u¯ )k ]u3 +[cu¯ (k +k u¯ )(k +k u¯ ) 33 11 33 13 3 31 11 1 23 i 3 11 13 3 22 23 3 +(k +k +k u¯ )[au¯ (k +k u¯ )+bu¯ (k +k u¯ )] 33 31 32 2 1 22 23 3 2 11 13 3 +(k +k u¯ )(ek u¯ u¯ )−ak k u¯ u¯ u¯ +dk u¯ u¯ (k +k u¯ ) 11 13 3 23 2 3 23 32 1 2 3 13 1 3 22 23 3 −k k bu¯ u¯ u¯ −k k u¯ u¯ −k k u¯ u¯2]µ2 +[acu¯ u¯ (k +k u¯ ) 13 31 1 2 3 31 22 1 2 31 32 1 3 i 1 3 22 23 3 +bcu¯ u¯ (k +k u¯ )+abu¯ u¯ (k +k u¯ +k u¯ )+eu¯ u¯ (k 2 3 11 13 3 1 2 33 31 1 32 2 2 3 11 +k u¯ )+au¯ (ek u¯ u¯ −k u¯ u¯ )+k u¯ u¯ u¯ bd 13 3 1 23 2 3 23 2 3 11 1 2 3 +du¯ u¯ (k +k u¯ −k bu¯ u¯ u¯ )]u +(abc+ae+bd)u¯ u¯ u¯ . (3.6) 1 3 22 23 3 31 1 2 3 i 1 2 3 Let λ ,λ ,λ be the three roots of (3.3). In order to obtain the stability of u¯, we need 1i 2i 3i to show that three exists a positive constant δ such that Re{λ },Re{λ },Re{λ } < −δ, for all i ≥ 1. (3.7) 1i 2i 3i The aim of the following theorem is to prove that the diffusion alone (with out cross- diffusion, i.e k = k = k = k = 0) can not drive instability for this model. 31 13 32 23 7 Theorem 3.1. Suppose that (2.3) holds and k = k = k = k = 0. Then the positive 13 31 23 32 equilibrium u¯ of (3.1) is linearly stable. Proof: Substituting k =k =k =k = 0 into (3.4), (3.5) and (3.6) we have 13 31 23 32 A = au¯ +bu¯ +cu¯ +(k +k +k )µ > 0 2i 1 2 3 11 22 33 i A = (k k +k k +k k )µ2 +[a(k +k )u¯ +b(k +k )u¯ +c(k +k )u¯ ]µ 1i 11 22 11 33 22 33 i 22 33 1 11 33 2 11 22 3 i +abu¯ u¯ +acu¯ u¯ +du¯ u¯ +bcu¯ u¯ +eu¯ u¯ > 0 1 2 1 3 1 3 2 3 2 3 A = k k k µ3 +(k k cu¯ +k k bu¯ +k k au¯ )µ2 0i 11 22 33 i 11 22 3 11 33 2 22 33 1 i +(abk u¯ u¯ +dk u¯ u¯ +ack u¯ u¯ +bck u¯ u¯ +ek u¯ u¯ +bdk u¯ u¯ u¯ )µ 33 1 2 22 1 2 22 1 3 11 2 3 11 2 3 11 1 2 3 i +(ae+bd+abc)u¯ u¯ u¯ > 0. 1 2 3 A direct calculation shows that A A − A > 0 for all i ≥ 1. It follows from Routh- 2i 1i 0i Hurwize criterion that all the three roots λ ,λ ,λ of ρ (λ) = 0 have negative real parts 1i 2i 3i i for each i ≥ 1. Let λ = µ ε, then i ρ (λ) = µ3ξ3 +A µ2ξ2 +A µ ξ +A = ρ˜(ξ). i i 2i i 1i i 0i i Since µ → ∞, as i → ∞, we have i ρ˜(ξ) ρ¯(ξ) = lim i = ξ3 +(k +k +k )ξ2 +(k k +k k +k k )ξ +k k k . i→∞ µ3 11 22 33 11 22 22 33 11 33 11 22 33 i ApplyingtheRouth-Hurwitzcriterionitfollowsthatthethreerootsξ ,ξ ,ξ ofρ¯(ξ) = 0all 1 2 3 ¯ havenegativerealparts. Thus, thereexistsapositiveconstantδ suchthatRe{ξ },Re{ξ }, 1 2 ¯ Re{ξ } ≤ −2δ. By continuity, we see that there exists i ≥ 1, such that µ > 1 and the 3 0 i0 ¯ ¯ ¯ three roots ξ ,ξ ,ξ of ρ˜(ξ) = 0 satisfy Re{ξ },Re{ξ }, Re{ξ } ≤ −µ δ ≤ −µ δ ≤ −δ i1 i2 i3 i i1 i2 i3 i i0 ˜ ˜ ¯ for any i ≥ i , Let −δ = max {Re{λ },Re{λ },Re{λ }} and δ = min{δ,δ}, Then 0 1≤i≤i0 i1 i2 i3 (3.7) holds. Consequently the equilibrium u¯ is linearly stable. Note that A > 0,A > 0,A > 0, and A A − A > 0 if k = k = 0 since 2i 1i 0i 2i 1i 0i 31 32 the possible negative terms all involve either k or k . By the same arguments as in 31 32 Theorem 3.1, we have 8 Theorem 3.2. Suppose that (2.3) holds and k = k = 0, Then the positive equilibrium 31 32 u¯ of (1.2) is linearly stable. Next we consider the Turing instability i.e. the stability of the positive equilibrium u¯ = (u¯ ,u¯ ,u¯ ) changing from stable, for the ODE dynamics (1.1), to unstable for the 1 2 3 PDE dynamics (1.2). Here we give sufficient conditions on cross-diffusion which drives the instability. k and k are chosen as variation parameters. 31 32 Theorem 3.3. (1) Suppose that au¯ −u¯ < 0. Consider k as the variation parameter, 1 3 31 then there exists a positive constant δ such that when k > δ , the equilibrium u¯ is 31 31 31 linearly unstable for some domain Ω. (2) Suppose that bu¯ −u¯ < 0. Consider k as the variation parameter, then there exists 2 3 32 a positive constant δ such that when k > δ , the equilibrium u¯ is linearly unstable for 32 32 32 some domain Ω. Proof: Denote A(µ) = −(C µ3 +C µ2 +C µ+C ) (3.8) 3 2 1 0 where C = [(k k +k k u¯ +k k u¯ +k u¯ u¯ +k k u¯ +k k u¯ u¯ )k 3 33 11 33 13 3 31 11 1 31 1 3 32 11 2 32 13 2 3 22 +(k k +k k u¯ +k k u¯ )k ] 33 11 33 13 3 31 11 1 23 C = [cu¯ (k +k u¯ )(k +k u¯ )+(k +k +k u¯ )[au¯ (k +k u¯ ) 2 3 11 13 3 22 23 3 33 31 32 2 1 22 23 3 +bu¯ (k +k u¯ )]+(k +k u¯ )(ek u¯ u¯ )−ak k u¯ u¯ u¯ 2 11 13 3 11 13 3 23 2 3 23 32 1 2 3 +dk u¯ u¯ (k +k u¯ )−k k bu¯ u¯ u¯ −k k u¯ u¯ −k k u¯ u¯2] 13 1 3 22 23 3 13 31 1 2 3 31 22 1 2 31 32 1 3 C = [abu¯ u¯ k (au¯ −u¯ )+au¯ u¯ k (bu¯ −u¯ )+acu¯ u¯ (k +k u¯ )+ 1 1 2 31 1 3 1 2 32 2 3 1 3 22 23 3 bcu¯ u¯ (k +k u¯ )+abu¯ u¯ k +eu¯ u¯ (k +k u¯ )+aek u¯ u¯ u¯ 2 3 11 13 3 1 2 33 2 3 11 13 3 23 1 2 3 +k bdu¯ u¯ u¯ +du¯ u¯ (k +k u¯ ) 11 1 2 3 1 3 22 23 3 C = (abc+ae+bd)u¯ u¯ u¯ . 0 1 2 3 9 Case 1: k is the variation parameter. 31 We assume that au¯ − u¯ < 0. The following arguments by continuation are based on 1 3 the fact that each root of the algebraic equation (3.8) is a continuous function of the variation parameter k . It is easy to prove that equation (3.8) has three real roots 31 (i) (i) (1) µ = µ (k ), i = 1,2,3 when k goes to infinity and they are lim µ (k ) < 0, 1 1 31 31 k31→∞ 1 31 (2) (3) lim µ (k ) = 0 andlim µ (k ) > 0. By continuation, there exists a positive k31→∞ 1 31 k31→∞ 1 31 constant δ such that when k > δ , C > 0 anddet(A(µ)) has three real roots. Because 31 31 31 1 C > 0 and C > 0, the mumber of sign changes of (3.8) is exactly two. Therefore by 3 0 Descartes’rule, the three real roots have the following properties: (1) (2) (3) (1) −∞ < µ < 0 < µ < µ < ∞, 1 1 1 (1) (2) (3) (2) det(A(µ)) > 0 if µ ∈ (−∞,µ )∪(µ ,µ ), 1 1 1 (1) (2) (3) (3) det(A(µ)) < 0 if µ ∈ (µ ,µ )∪(µ ,∞). 1 1 1 (2) (3) If µ ∈ (µ ,µ ) for some i, then det(A(µ )) > 0 by (2), and consequently A = i 1 1 i 0i −det(A ) < 0. The number of sign of changes of the characteristic polynomial (3.3) (i) ρ (λ) = λ3+A λ2+A λ+A is either one or three. By Descartes’rule, the characteristic i 2i 1i 0i polynomial (3.3) has at least one positive eigenvalue. Hence, the equilibrium u¯ of (1.2) is linearly unstable for any domain Ω on which at least one eigenvalue µ of −∆ is in the i (2) (3) interval (µ ,µ ). 1 1 Case 2: k is the variation parameter. 32 We assume that bu¯ − u¯ < 0. The following arguments by continuation are based on 2 3 the fact that each root of the equation (3.8) is a continuous function of the variation k . 32 (i) (i) It is easy to prove that equation (3.8) has three real roots µ = µ (k ), i = 1,2,3 2 2 32 (1) (2) when k goes to infinity and they are lim µ (k ) < 0, lim µ (k ) = 0 and 32 k32→∞ 2 32 k32→∞ 2 32 (3) lim µ (k ) > 0. By continuation, there exists a positive constant δ such that k32→∞ 2 32 32 when k > δ , C > 0 and det(A(µ)) has three real roots. Because C > 0 and C > 0, 32 32 1 3 0 the mumber of sign changes of (3.8) is exactly two. Therefore by Descartes’rule,the three real roots have the folloeing properties: (1) (2) (3) (1) −∞ < µ < 0 < µ < µ < ∞, 2 2 2 (2) (2) (3) (2) det(A(µ)) > 0 if µ ∈ (−∞,µ ) (µ ,µ ), 2 2 2 S 10