- Email: [email protected]
Turing patterns in a reaction–diffusion system modeling hunting cooperation
Accepted Manuscript Turing patterns in a reaction–diffusion system modeling hunting cooperation F. Capone, M.F. Carfora, R. De Luca, I. Torcicollo
PII: DOI: Reference:
S0378-4754(19)30106-5 https://doi.org/10.1016/j.matcom.2019.03.010 MATCOM 4766
To appear in:
Mathematics and Computers in Simulation
Received date : 20 November 2018 Revised date : 26 February 2019 Accepted date : 17 March 2019 Please cite this article as: F. Capone, M.F. Carfora, R. De Luca et al., Turing patterns in a reaction–diffusion system modeling hunting cooperation, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.03.010 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Manuscript Click here to view linked References
Turing patterns in a reaction-diffusion system modeling hunting cooperation F. Caponea , M.F. Carforab , R. De Lucaa , I. Torcicollob,∗ a
Department of Mathematics and Applications “R.Caccioppoli”, University of Naples Federico II, Via Cintia, Naples, Italy b Istituto per le Applicazioni del Calcolo “Mauro Picone”, Via P. Castellino 111, CNR, Naples, Italy
Abstract A reaction-diffusion system governing the prey-predator interaction with hunting cooperation is investigated. Definitive boundedness of solutions is proved via the existence of positive invariants and attractive sets. Linear stability of the coexistence equilibria is performed and conditions guaranteeing the occurrence of Turing instability are found. Numerical simulations on the obtained results are provided. Keywords: Reaction-diffusion, Hunting cooperation, Predator-prey, Turing, Routh-Hurwitz, Stability 1. Introduction Population dynamics is a very active research field aimed to analyze how different dynamical behaviors may affect the ecosystem evolution. In predator-prey interactions, a crucial role is played by the cooperation which can help the survival of individuals. In fact, hunting cooperation, observed in many species (carnivores, birds, aquatic organisms, spiders), implies, as a direct consequence, Allee effect in predators which allows these last to persist even if the prey population does not sustain them in the absence of hunting cooperation [1]. In many papers addressing population dynamics it is ∗
Corresponding author Email addresses: [email protected] (F. Capone), [email protected] (M.F. Carfora), [email protected] (R. De Luca), [email protected] (I. Torcicollo) Preprint submitted to Mathematics and Computers in Simulation
February 26, 2019
assumed that population is homogeneously distributed in the environment, leading to the formulation of ODEs systems [1, 2, 3, 6, 11, 23]. However, in several biological contexts, the formation of spatial patterns which evolve in time is observed [25] and hence, in order to model this phenomenon, it is necessary that the species diffuse in the environment leading to a PDE system. From a mathematical point of view, in order to explore the formation of patterns, one has to look for conditions guaranteeing that a stable equilibrium in the absence of diffusion becomes unstable in the presence of diffusion (Turing instability or diffusion-driven instability). Calling “activator” the species promoting growth and “inhibitor” the one repressing growth, it is well known that a necessary condition for the occurrence of Turing instability is that the inhibitor has to diffuse faster than the activator [16]. Turing patterns in predator-prey models have been deeply investigated in literature. As examples, in [5] the case of a predator-prey system with linear diffusion and Robin boundary conditions has been analyzed; in [22] the herd behavior and cross-diffusion have been considered; in [27] a nonlinear diffusion has been introduced; in [28] the influence of self and cross-diffusion and Beddington-DeAngelis functional response has been investigated. In the recent years, there has been a growing interest in fractional diffusion models, obtaining by replacing the time derivative by a non-integer order derivative (i.e., an integro-differential operator). Indeed, memory effects in dynamical systems are well described by nonlinear reaction diffusion models with time fractional derivatives (see [12] for an extended survey). Several authors investigated pattern formation in such systems ([9], [10]), mainly by computer simulations, due to the lack of a completely developed stability theory for these models. More recently, Yin and Wen [26] show that fractional-derivative systems can form spatial patterns even though its firstderivative counterpart can’t exhibit any steady pattern. Such mechanism of pattern formation induced by anomalous diffusion suggests new directions of research. In this paper, however, we limit our investigation to usual diffusion and explore the Turing patterns formation in a predator-prey system with hunting cooperation on introducing a suitable reaction-diffusion model and look for conditions guaranteeing that a biological meaningful equilibrium, stable in the absence of diffusion, becomes unstable in the presence of diffusion. We claim that such a study provides challenges and ideas in many other fields of applied mathematics such as structural engineering, ecology, aerospace science and economics in which nonlinear mathematical models having a 2
similar structure are considered (see for instance [5, 6, 7, 8, 18, 19, 20, 21, 23, 24], and references therein). The plan of the paper is as follows. Section 2 introduces some preliminaries results on the mathematical model. Linear stability analysis is performed in Section 3 and Turing instability is investigated in Section 4. After Section 5 in which numerical simulations are shown, the paper ends with some concluding remarks. 2. Mathematical model We extend the classical Lotka-Volterra model with logistic growth for the prey including hunting cooperation [1] by assuming that both preys and predators can linearly spread in the environment, i.e. N ∂N = rN 1 − − (λ + aP )N P + d1 ∆N, ∂τ K (1) ∂P = e(λ + aP )N P − mP + d2 ∆P, ∂τ
where N (X, τ ) and P (X, τ ) denote, respectively, the prey and predator densities; (X, τ ) ∈ Ω × R+ with Ω a bounded, connected, open subset of Rn (n = 1, 2, 3) having the internal cone property; r(> 0) is the per capita intrinsic growth rate of prey; K(> 0) is the carrying capacity of prey; e(> 0) is the food conversion efficiency; m(> 0) is the per capita mortality rate of predators; λ(> 0) is the attack rate; a(> 0) is the predator cooperation in hunting rate; di (> 0) are the diffusion coefficients. To (1) we append homogeneous Neumann (no-flux) boundary conditions ∇N · n = 0,
∇P · n = 0
on ∂Ω × R+ ,
(2)
being n the outward unit normal to ∂Ω and smooth positive initial data N (X, 0) = N0 (X),
P (X, 0) = P0 (X),
X ∈ Ω.
Using the non-dimensional variables adopted in [1] and setting X = A the diameter of Ω), model (1)-(3) becomes ∂n = σn 1 − n − (1 + αp)np + γ1 ∆n, ∂t k ∂p = (1 + αp)np − p + γ2 ∆p, ∂t 3
(3) x , (with A
(4)
( ∇n · n = 0, ∇p · n = 0, on ∂Ω × R+ . n(x, 0) = n0 (x), p(x, 0) = p0 (x), on Ω.
(5)
with
r eλK am di , k= , α = 2 , γi = 2 , i = 1, 2. (6) m m λ Am As observed in [1], k comprises the dimensional carrying capacity, the conversion efficiency as well as the per-capita predator mortality and attack rate; then it can be interpreted as the predator basic reproduction number. Denoting by k·k and k·k∞ the L2 (Ω)-norm and the L∞ (Ω)−norm, respectively; by T > 0 an arbitrary fixed time, by ΩT = Ω × (0, T ] the parabolic ¯ × [0, T ]), by ΓT = ∂Ω × [0, T ) cylinder (ΩT being the parabolic interior of Ω ˜ T = ΓT ∪ {Ω × {t = 0}}, the following theorems hold. and by Γ σ=
¯ T )]2 be the non negative solution of Theorem 1. Let (n, p) ∈ [C12 (ΩT ) ∩ C(Ω (4)-(5). Then n(x, t) ≤ max k, max n0 (x) := M1 , p(x, t) ≤ C∞ (p0 (x)), a.e. in Ω ¯ Ω
(7)
where C∞ is a positive constant depending on the initial data. Proof. Setting max n(x, t) = n(x1 , t1 ), ¯T Ω
(8)
we have to distinguish three cases. 1) If (x1 , t1 ) belongs to the interior of ΩT , then (4)1 implies that ∂n n − σn 1 − − γ1 ∆n < 0. ∂t k (x1 ,t1 )
∂n Since ∂t
(x1 ,t1 )
= 0, [∆n](x1 ,t1 ) < 0, (9) can hold if n(x1 , t1 ) < k.
2) If (x1 , t1 ) ∈ Ω × {t = 0} then n(x1 , t1 ) ≤ max n0 (x). ¯ Ω
4
(9)
3) If (x1 , t1 ) ∈ ΓT , in view of the regularity of the domain Ω, since Ω satisfies in any point x0 ∈ ∂Ω the interior ball condition, there exists an open ball B ∗ ⊂ Ω with x0 ∈ ∂B ∗ . If n(x1 , t1 ) > k, on choosing the radius of B ∗ sufficiently small, it follows that γ1 ∆n −
∂n > 0, ∂t
in B ∗ ,
and by virtue of Hopf’s Lemma [17], one obtains that dn > 0, dn (x1 ,t1 ) which is not possible in view of (5)1 . In order to prove (7)2 , in view of (7)1 and (4)2 , it turns out that p(x, t) is a sub-solution of the problem ∂U ∂t − γ2 ∆U = M1 (1 + αU )U, (10) ∇U · n = 0, on ∂Ω × R+ , U (x, 0) = U0 (x) = max p0 (x). ¯ Ω
Since
1 M2 M1 (1 + αU )U ≤ αM1 + U2 + 1 , 2 2
(11)
in view of Theorem 1 of [13], one obtains that, denoting by τ (U0 ) the maximal existence time of the solution U (x, t) of (10), since – from the continuous dependence on the initial data – there exists a positive constant C(U0 ) such that kU (·, t)k ≤ C(U0 ), ∀t ∈ (0, τ (U0 )), (12) the solution U (x, t) exists for all time and there exists a positive constant C∞ such that kU (·, t)k∞ ≤ C∞ (U0 ), ∀t > 0 (13) and the thesis is proved. Theorem 2. ∀ε > 0 the manifold n a ¯o Σε= (n, p) ∈ [R+ ]2 : knk2 +kpk2 ≤ (1 + ε) 2 5
(14)
with
2 a ¯ = 2|Ω| (σ + 1)M12 + M1 C∞ (1 + αC∞ ) ,
(15)
and |Ω| being the measure of Ω, is an absorbing set for system (4). Proof. Multiplying (4)1 by n, (4)2 by p, adding the resulting equations and integrating over Ω, by virtue of the divergence theorem, the boundary conditions (5)1 and (7), it turns out that 1 d(knk2 + kpk2 ) 2 ≤ σ knk2 − kpk2 + M1 C∞ |Ω|(1 + αC∞ ). 2 dt
(16)
Then, setting E = knk2 + kpk2 , it follows that dE ≤ −2E + a ¯. dt
(17)
Following the procedure in [14], one can prove that Σ is an absorbing set. The coexistence equilibria E ∗ = (n∗ , p∗ ) are given by n∗ =
1 1 + αp∗
(18)
and p∗ (positive) solution of f (x) := kα2 x3 + 2αkx2 + k(1 − σα)x + σ(1 − k) = 0. (19) √ −1 + 1 + 3σα 1 < k ≤ 1, σ > there exIf k > 1, E ∗ is unique, while if ασ α ist two coexistence equilibria. In the all other cases, no coexistence equilibria are admissible. Introducing the perturbation fields {X = n − n∗ , Y = p − p∗ }, (4) becomes ∂ X a11 a12 X γ1 ∆X F1 = + + , (20) a21 a22 Y γ2 ∆Y F2 ∂t Y with
σ 1 + 2αp∗ a = − , a = − , 11 12 k(1 + αp∗ ) 1 + αp∗ ∗ a = (1 + αp∗ )p∗ , a = αp , 21 22 1 + αp∗ 2 σX − (1 + αp∗ )XY − αY (XY + n∗ Y + p∗ X), F 1 = − k F2 = (1 + αp∗ )XY + αY (XY + n∗ Y + p∗ X). 6
(21)
To (20) we associate the initial-boundary conditions ( ∇X · n = ∇Y · n = 0, on ∂Ω × R+ . X(x, 0) = X0 (x), Y (x, 0) = Y0 (x), x ∈ Ω.
(22)
Denoting by W ∗ (Ω) the functional space defined by dϕ ∗ 1,2 1,2 + , W (Ω) = ϕ ∈ W (Ω) ∩ W (∂Ω) : = 0 on ∂Ω×R dn our aim is to study the stability of (n∗ , p∗ ) with respect to the perturbations (X, Y ) ∈ [W ∗ (Ω)]2 . 3. Linear stability of the coexistence equilibria Let 0 = α0 < α1 < α2 < · · · < αn < · · ·
(23)
∂X = LX, ∂t
(24)
be the eigenvalues of the operator −∆ on Ω with homogeneous Neumann boundary conditions. The linearized version of (20) is
where L=
a11 + γ1 ∆ a12 a21 a22 + γ2 ∆
,
X=
X Y
.
(25)
For each i = 0, 1, 2, 3, ..., λi is an eigenvalue of L if and only if λi is an eigenvalue of the matrix (see [4]) a11 − αi γ1 a12 ˜ . (26) Li = a21 a22 − αi γ2 The characteristic equation of L˜i is λ2i − I1i λi + I2i = 0, where Iji , (j = 1, 2), are the principal L˜i −invariants given by I1i = trL˜i = I01 − αi (γ1 + γ2 ), σγ − kαp∗ γ1 I2i = det L˜i == γ1 γ2 αi2 + 2 αi + I02 , ∗ k(1 + αp ) 7
(27)
(28)
I0j , (j = 1, 2), being the principal invariants in absence of diffusion, i.e. kαp∗ − σ , I01 = a11 + a22 = k(1 + αp∗ ) (29) σα ∗ 0 ∗ + 1 + 2αp . I2 = a11 a22 − a12 a21 = p − k(1 + αp∗ )2
It is to recall that, in the absence of diffusion, the necessary and sufficient condition guaranteeing the linear stability of E ∗ is {I01 < 0, I02 > 0} [15]. Simple calculation shows that E ∗ is linearly stable in the absence of diffusion if and only if k(1 + αp∗ )2 (1 + 2αp∗ ) kαp∗ < σ < , (30) α holds. Let us underline that, when there exist two coexistence equilibria, (30) implies that only the one with larger predator density is linearly stable. Moreover, (30) is a necessary condition for linear stability in the presence of diffusion too. Theorem 3. If either k(1 + αp∗ )2 (1 + 2αp∗ ) kαp∗ < σ < , α ∗ γ kαp 2 ≥ , γ1 σ
(31)
or
k(1 + αp∗ )2 (1 + 2αp∗ ) kαp∗ < σ < , α 2 ∗ ∗ kαp (kγ1 αp −σγ2 ) γ 2< , < 4kp∗ k(1+2αp∗ )(1+αp∗ )2 −σα , γ1 σ γ1 γ2
(32)
holds, then the coexistence equilibrium (when it exists) is linearly stable. Proof. When (30) holds, I01 < 0 and hence I1i < 0, ∀i. In view of (28)2 , it turns out that either (31) or (32) guarantees that I2i > 0, ∀i. 4. Turing instability In this Section we look for conditions guaranteeing that E ∗ – stable in the absence of diffusion – becomes unstable in the presence of diffusion (Turing 8
instability). Since I01 < 0 ⇒ I1i < 0, ∀i, Turing instability may occur only if I02 > 0 together with I2¯i < 0 for at least one index ¯i. Moreover, in view of (28)2 , a necessary condition for the occurrence of Turing instability is kαp∗ γ2 < . γ1 σ
(33)
It should be noted that in the proposed model, prey is the activator for the predator population (a21 > 0) but it is self-inhibiting (a11 < 0). Predators are inhibitors for the prey (a12 < 0), but they are also self-activating (a22 > 0). In this context, a necessary condition for diffusion-driven instability is that prey disperses faster than predators. In fact, in view of (30), it turns out kαp∗ < 1 and (33) implies necessarily that that σ γ1 > γ2 .
(34)
As illustrated in details by Murray in [16] (see Sect. 2.3, pp. 87-89), biologically there are two types of interaction with consequent qualitatively different reactions: predators spread faster than prey or vice versa. Although a large proportion of the bibliography falls in the first case, the other one is not an uncommon situation and it is precisely the object of the present study. Theorem 4. If (30) and (33) hold together with either αp∗ , γ < 2 α1 (1 + αp∗ ) k(1 + αp∗ )I02 + σα1 γ2 γ1 > kα1 [αp∗ − α1 γ2 (1 + αp∗ )]
(35)
or
(kγ1 αp∗ − σγ2 )2 > 4kp∗ k(1 + 2αp∗ )(1 + αp∗ )2 − σα , γ1 γ2 then Turing instability occurs.
(36)
Proof. If (30) holds, then there is stability in the absence of diffusion. In view of (28)2 , it follows that αp∗ σα1 γ2 + I02 k(1 + αp∗ ) α γ + I21 = α1 γ2 − , (37) 1 1 1 + αp∗ k(1 + αp∗ ) that is negative if (33) and (35), or (33) and (36) hold. 9
10
10
9
0.708
9
1.808
8
0.706
8
1.806
7
0.704
7
1.804
6
0.702
6
1.802
5
0.7
5
1.8
4
0.698
4
1.798
3
0.696
3
1.796
2
0.694
2
1.794
1
0.692
1
1.792
0
0 0
2
4
6
8
10
0
2
4
6
8
10
Figure 1: Plot of the initial state of system (4) in the 2D setting. The initial values are fixed as (n0 , p0 ) = (0.7,1.8) plus a random perturbation of amplitude 0.01. Left panel: initial spatial distribution of prey; right panel: initial spatial distribution of predators.
5. Numerical simulations In this paper, a theoretical analysis of an evolutionary process that involves the interaction of two spatially distributed population, predators and prey, with local diffusion, has been performed. In order to highlight the theoretical results, in this section we provide some numerical simulations which show the impact of different choices for the diffusion coefficients on the dynamics of the two populations. First, we consider a configuration in the parameters space that leads to a single stable coexistence equilibrium in absence of diffusion, according to (30). Specifically, we choose α = 0.3 (weak hunting cooperation among predators), σ = 3 (per capita growth rate of prey) and k = 5 (high predator basic reproduction number for homogeneously mixed populations) and the unique stable coexistence equilibrium is given by E ∗ ≈ (0.66,1.72). In the PDE model (4) we assume an initial condition very close to this equilibrium value plus a small random perturbation. As well known [16], the size of the space domain has to be chosen sufficiently large with respect to the wavelength of the unstable modes, in order to see the instability patterns in the solution. Then we set as space domain the interval [0,10] for 1D simulations, or the square [0,10] × [0,10] for 2D simulations, with a discretization step h = 0.05. Figure 1 shows the chosen initial condition of the system for the prey (left panel) and predator (right panel) populations in the 2D setting. Now, if the diffusion coefficients are chosen, according to Theorem 4, in the range given by (33)-(35), the conditions for instability to emerge are 10
10
10
9
0.6595
1.760
9 1.750
8
8 0.6590
7
7
1.740
0.6585
6
6
5
0.6580
4
1.730
5 4
0.6575
3
1.720
3 0.6570
2 1
0.6565
2
1.710
1 1.700
0
0 0
2
4
6
8
10
0
2
4
6
8
10
Figure 2: Plot of the state of the system for T = 100 in presence of (unbalanced) diffusion: γ1 = 1.5; γ2 = 0.01 with parameters choice α=0.3, σ=3, k=5. Note that both populations (prey, in the left panel; predators in the right panel) start showing spatial patterns.
satisfied. Specifically, we fix γ1 = 1.5 and γ2 = 0.01. Then, the coexistence equilibrium loses stability due to the Turing effect: as shown in Figure 2, spatial patterns emerge in both populations; these patterns persist and increase in amplitude with time (Figure 3). However, it is worth noting that the amplitude of such patterns remains confined within the absorbing set given by (14). This result can be very helpful in analyzing and also forecasting the spatial redistribution of populations in real world applications. 6. Conclusions Ecologists increasingly emphasize how ecosystems reveal spatial self-organization, where large-scale ordered spatial patterns emerge from disordered initial conditions through local interactions. Indeed, organisms concentrate resources in their local environment, enabling persistence even when mean field resource levels become too low for their survival. Ecosystems with regular patterns may be more resilient and resistant to the environmental change compared to homogeneous ecosystems. This process is the key to understand ecological stability and diversity. In view of these real world applications, we have presented a theoretical analysis of an evolutionary process that involves the interaction of two spatially distributed population, predators and prey, with local diffusion. Our analysis, also illustrated by numerical simulations, reveals the impact of diffusion and hunting cooperation on the overall dynamics. We have investigated the stability of the coexistence equilibria and found the conditions guaranteeing the occurrence of Turing instability, lead11
10
10 0.7500
9
6.000
9 0.7000
8 7
0.6500
6
8
5.000
7 4.000
6 0.6000
5
5 3.000
4
0.5500
3
4 3
2.000
0.5000
2
2
1
0.4500
0
1.000
1 0
0
2
4
6
8
10
0
2
4
6
8
10
Figure 3: Plot of the state of the system for T =500 in presence of (unbalanced) diffusion: γ1 = 1.5; γ2 = 0.01 and parameters choice as in the previous Figure. Note that both populations (prey, in the left panel; predators in the right panel) show spatial patterns of great amplitude.
ing to the formation of isolated group. Specifically, we identified a scenario in which weak hunting cooperation allows predators and prey population to spatially redistribute to form isolated (spotted) groups, where predators density rises well above the level sustained by the same prey population in the absence of diffusion. All these results can have important applications for global environmental change, ecosystem adaptation and biological control using predators. Acknowlegements This paper has been performed under the auspices of the G.N.F.M. of INdAM. References [1] M. Alves, F. Hilker, Hunting cooperation and Allee effects in predators, J. Theor. Biol. 419 (2017) 13–22. [2] L. Berec, Impacts of foraging facilitation among predators on predatorprey dynamics, Bull. Math. Biol. 72 (2010) 94–121. [3] D. Boukal, M. Sabelis, L. Berec, How predator functional responses and Allee effects in prey affect the paradox of enrichment and population collapses, Theor. Popul. Biol. 72 (2007) 136–147. 12
[4] R. Cantrell, C. Cosner, Spatial ecology via reaction-diffusion equations, John Wiley and Sons, Ltd, Chichester, UK, 2003. [5] F. Capone, On the dynamics of predator-prey models with the Beddington-De Angelis functional response, under robin boundary conditions, Ricerche Mat. 57 (2008) 137–157. [6] F. Capone, M. Carfora, R. De Luca, I. Torcicollo, On the dynamics of an intraguild predator–prey model, Math. Comput. Simulation 149 (2018) 17–31. [7] F. Capone, V. De Cataldis, R. De Luca, On the stability of a SEIR reaction diffusion model for infections under Neumann boundary conditions, Acta Appl. Math. 132 (2014) 165–314. [8] F. Capone, R. De Luca, On the nonlinear dynamics of an ecoepidemic reaction-diffusion model, Int. J. Nonlinear Mech. 95 (2017) 307–314. [9] B. Datsko, V. Gafiychuk, Complex nonlinear dynamics in subdiffusive activator–inhibitor systems., Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 1673–1680. [10] B. Datsko, Y. Luchko, V. Gafychuk, Pattern formation in fractional reaction difusion systems with multiple homogeneous states, Int. J. Bifurcation Chaos. 22, (2012). 22 (2012) 1250087. [11] R. De Luca, On the long-time dynamics of nonautonomous predatorprey models with mutual interference, Ricerche Mat. 61 (2012) 275–290. [12] M. Du, Z. Wang, H. Hu, Measuring memory with the order of fractional derivative, Scientific Reports 3 (2013) 3431. [13] L. Dung, Dissipativity and global attractors for a class of quasilinear parabolic systems, Commun. Partial Diff. Eqns. 22 (1997) 213–433. [14] J. Flavin, S. Rionero, Qualitative Estimates for Partial Differential Equations: An Introduction, volume 8111, CRC Press Inc., Boca Raton, 1996. [15] D. Merkin, Introduction to the Theory of Stability, volume 24, Text in Applied Mathematic, Springer, 1997. 13
[16] J. Murray, Mathematical Biology II: Spatial Models and Biomedical Applications, 3rd ed., Springer-Verlag, Berlin Heidelberg, 1993. [17] M. Protter, H. Weinberger, Maximum Principles in Differential Equations, Springer, New York, 1999. [18] S. Rionero, A peculiar Lyapunov functional for ternary reaction-diffusion dynamical systems, Boll. U.M.I. 9 (2011) 393–407. [19] S. Rionero, Stability of ternary reaction-diffusion dynamical systems, Rend. Lincei Mat. Appl. 22 (2011) 245–268. [20] S. Rionero, I. Torcicollo, Stability of a continuous reaction-diffusion Cournot-Kopel duopoly game model, Acta Appl. Math. 132 (2014) 505– 513. [21] S. Rionero, I. Torcicollo, On the dynamics of a nonlinear reaction– diffusion duopoly model, Int. J. Non-Linear Mech. 99 (2018) 105–111. [22] X. Tang, Y. Song, T. Zhang, Turing-Hopf bifurcation analysis of a predator-prey model with herd behavior and cross-diffusion, Nonlinear Dyn. 86 (2016) 73–89. [23] A. Terry, Predator-prey models with component Allee effect for predator reproduction, J. Math. Biol. 71 (2015) 1325–1352. [24] I. Torcicollo, On the non-linear stability of a continuous duopoly model with constant conjectural variation, Int. J Non-Linear Mech. 81 (2016) 268–273. [25] A. Turing, The Chemical Basis for Morphogenesis, Philos. Trans. R. Soc. Lond. Ser. B, Biol. Sci. 237 (1952) 37–72. [26] H. Yin, X. Wen, Pattern formation through temporal fractional derivatives, Scientific Reports 8 (2018) 5070. [27] H. Yin, X. Xiao, X. Wen, Turing patterns in a predator-prey system with self-diffusion, Abstract and Applied Analysis 2013 (2013). [28] J. Zhang, W. Li, Y. Wang, Turing patterns of a strongly coupled predator-prey system with diffusion effects, Nonlinear Anal. 74 (2011) 847–858. 14
PII: DOI: Reference:
S0378-4754(19)30106-5 https://doi.org/10.1016/j.matcom.2019.03.010 MATCOM 4766
To appear in:
Mathematics and Computers in Simulation
Received date : 20 November 2018 Revised date : 26 February 2019 Accepted date : 17 March 2019 Please cite this article as: F. Capone, M.F. Carfora, R. De Luca et al., Turing patterns in a reaction–diffusion system modeling hunting cooperation, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.03.010 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Manuscript Click here to view linked References
Turing patterns in a reaction-diffusion system modeling hunting cooperation F. Caponea , M.F. Carforab , R. De Lucaa , I. Torcicollob,∗ a
Department of Mathematics and Applications “R.Caccioppoli”, University of Naples Federico II, Via Cintia, Naples, Italy b Istituto per le Applicazioni del Calcolo “Mauro Picone”, Via P. Castellino 111, CNR, Naples, Italy
Abstract A reaction-diffusion system governing the prey-predator interaction with hunting cooperation is investigated. Definitive boundedness of solutions is proved via the existence of positive invariants and attractive sets. Linear stability of the coexistence equilibria is performed and conditions guaranteeing the occurrence of Turing instability are found. Numerical simulations on the obtained results are provided. Keywords: Reaction-diffusion, Hunting cooperation, Predator-prey, Turing, Routh-Hurwitz, Stability 1. Introduction Population dynamics is a very active research field aimed to analyze how different dynamical behaviors may affect the ecosystem evolution. In predator-prey interactions, a crucial role is played by the cooperation which can help the survival of individuals. In fact, hunting cooperation, observed in many species (carnivores, birds, aquatic organisms, spiders), implies, as a direct consequence, Allee effect in predators which allows these last to persist even if the prey population does not sustain them in the absence of hunting cooperation [1]. In many papers addressing population dynamics it is ∗
Corresponding author Email addresses: [email protected] (F. Capone), [email protected] (M.F. Carfora), [email protected] (R. De Luca), [email protected] (I. Torcicollo) Preprint submitted to Mathematics and Computers in Simulation
February 26, 2019
assumed that population is homogeneously distributed in the environment, leading to the formulation of ODEs systems [1, 2, 3, 6, 11, 23]. However, in several biological contexts, the formation of spatial patterns which evolve in time is observed [25] and hence, in order to model this phenomenon, it is necessary that the species diffuse in the environment leading to a PDE system. From a mathematical point of view, in order to explore the formation of patterns, one has to look for conditions guaranteeing that a stable equilibrium in the absence of diffusion becomes unstable in the presence of diffusion (Turing instability or diffusion-driven instability). Calling “activator” the species promoting growth and “inhibitor” the one repressing growth, it is well known that a necessary condition for the occurrence of Turing instability is that the inhibitor has to diffuse faster than the activator [16]. Turing patterns in predator-prey models have been deeply investigated in literature. As examples, in [5] the case of a predator-prey system with linear diffusion and Robin boundary conditions has been analyzed; in [22] the herd behavior and cross-diffusion have been considered; in [27] a nonlinear diffusion has been introduced; in [28] the influence of self and cross-diffusion and Beddington-DeAngelis functional response has been investigated. In the recent years, there has been a growing interest in fractional diffusion models, obtaining by replacing the time derivative by a non-integer order derivative (i.e., an integro-differential operator). Indeed, memory effects in dynamical systems are well described by nonlinear reaction diffusion models with time fractional derivatives (see [12] for an extended survey). Several authors investigated pattern formation in such systems ([9], [10]), mainly by computer simulations, due to the lack of a completely developed stability theory for these models. More recently, Yin and Wen [26] show that fractional-derivative systems can form spatial patterns even though its firstderivative counterpart can’t exhibit any steady pattern. Such mechanism of pattern formation induced by anomalous diffusion suggests new directions of research. In this paper, however, we limit our investigation to usual diffusion and explore the Turing patterns formation in a predator-prey system with hunting cooperation on introducing a suitable reaction-diffusion model and look for conditions guaranteeing that a biological meaningful equilibrium, stable in the absence of diffusion, becomes unstable in the presence of diffusion. We claim that such a study provides challenges and ideas in many other fields of applied mathematics such as structural engineering, ecology, aerospace science and economics in which nonlinear mathematical models having a 2
similar structure are considered (see for instance [5, 6, 7, 8, 18, 19, 20, 21, 23, 24], and references therein). The plan of the paper is as follows. Section 2 introduces some preliminaries results on the mathematical model. Linear stability analysis is performed in Section 3 and Turing instability is investigated in Section 4. After Section 5 in which numerical simulations are shown, the paper ends with some concluding remarks. 2. Mathematical model We extend the classical Lotka-Volterra model with logistic growth for the prey including hunting cooperation [1] by assuming that both preys and predators can linearly spread in the environment, i.e. N ∂N = rN 1 − − (λ + aP )N P + d1 ∆N, ∂τ K (1) ∂P = e(λ + aP )N P − mP + d2 ∆P, ∂τ
where N (X, τ ) and P (X, τ ) denote, respectively, the prey and predator densities; (X, τ ) ∈ Ω × R+ with Ω a bounded, connected, open subset of Rn (n = 1, 2, 3) having the internal cone property; r(> 0) is the per capita intrinsic growth rate of prey; K(> 0) is the carrying capacity of prey; e(> 0) is the food conversion efficiency; m(> 0) is the per capita mortality rate of predators; λ(> 0) is the attack rate; a(> 0) is the predator cooperation in hunting rate; di (> 0) are the diffusion coefficients. To (1) we append homogeneous Neumann (no-flux) boundary conditions ∇N · n = 0,
∇P · n = 0
on ∂Ω × R+ ,
(2)
being n the outward unit normal to ∂Ω and smooth positive initial data N (X, 0) = N0 (X),
P (X, 0) = P0 (X),
X ∈ Ω.
Using the non-dimensional variables adopted in [1] and setting X = A the diameter of Ω), model (1)-(3) becomes ∂n = σn 1 − n − (1 + αp)np + γ1 ∆n, ∂t k ∂p = (1 + αp)np − p + γ2 ∆p, ∂t 3
(3) x , (with A
(4)
( ∇n · n = 0, ∇p · n = 0, on ∂Ω × R+ . n(x, 0) = n0 (x), p(x, 0) = p0 (x), on Ω.
(5)
with
r eλK am di , k= , α = 2 , γi = 2 , i = 1, 2. (6) m m λ Am As observed in [1], k comprises the dimensional carrying capacity, the conversion efficiency as well as the per-capita predator mortality and attack rate; then it can be interpreted as the predator basic reproduction number. Denoting by k·k and k·k∞ the L2 (Ω)-norm and the L∞ (Ω)−norm, respectively; by T > 0 an arbitrary fixed time, by ΩT = Ω × (0, T ] the parabolic ¯ × [0, T ]), by ΓT = ∂Ω × [0, T ) cylinder (ΩT being the parabolic interior of Ω ˜ T = ΓT ∪ {Ω × {t = 0}}, the following theorems hold. and by Γ σ=
¯ T )]2 be the non negative solution of Theorem 1. Let (n, p) ∈ [C12 (ΩT ) ∩ C(Ω (4)-(5). Then n(x, t) ≤ max k, max n0 (x) := M1 , p(x, t) ≤ C∞ (p0 (x)), a.e. in Ω ¯ Ω
(7)
where C∞ is a positive constant depending on the initial data. Proof. Setting max n(x, t) = n(x1 , t1 ), ¯T Ω
(8)
we have to distinguish three cases. 1) If (x1 , t1 ) belongs to the interior of ΩT , then (4)1 implies that ∂n n − σn 1 − − γ1 ∆n < 0. ∂t k (x1 ,t1 )
∂n Since ∂t
(x1 ,t1 )
= 0, [∆n](x1 ,t1 ) < 0, (9) can hold if n(x1 , t1 ) < k.
2) If (x1 , t1 ) ∈ Ω × {t = 0} then n(x1 , t1 ) ≤ max n0 (x). ¯ Ω
4
(9)
3) If (x1 , t1 ) ∈ ΓT , in view of the regularity of the domain Ω, since Ω satisfies in any point x0 ∈ ∂Ω the interior ball condition, there exists an open ball B ∗ ⊂ Ω with x0 ∈ ∂B ∗ . If n(x1 , t1 ) > k, on choosing the radius of B ∗ sufficiently small, it follows that γ1 ∆n −
∂n > 0, ∂t
in B ∗ ,
and by virtue of Hopf’s Lemma [17], one obtains that dn > 0, dn (x1 ,t1 ) which is not possible in view of (5)1 . In order to prove (7)2 , in view of (7)1 and (4)2 , it turns out that p(x, t) is a sub-solution of the problem ∂U ∂t − γ2 ∆U = M1 (1 + αU )U, (10) ∇U · n = 0, on ∂Ω × R+ , U (x, 0) = U0 (x) = max p0 (x). ¯ Ω
Since
1 M2 M1 (1 + αU )U ≤ αM1 + U2 + 1 , 2 2
(11)
in view of Theorem 1 of [13], one obtains that, denoting by τ (U0 ) the maximal existence time of the solution U (x, t) of (10), since – from the continuous dependence on the initial data – there exists a positive constant C(U0 ) such that kU (·, t)k ≤ C(U0 ), ∀t ∈ (0, τ (U0 )), (12) the solution U (x, t) exists for all time and there exists a positive constant C∞ such that kU (·, t)k∞ ≤ C∞ (U0 ), ∀t > 0 (13) and the thesis is proved. Theorem 2. ∀ε > 0 the manifold n a ¯o Σε= (n, p) ∈ [R+ ]2 : knk2 +kpk2 ≤ (1 + ε) 2 5
(14)
with
2 a ¯ = 2|Ω| (σ + 1)M12 + M1 C∞ (1 + αC∞ ) ,
(15)
and |Ω| being the measure of Ω, is an absorbing set for system (4). Proof. Multiplying (4)1 by n, (4)2 by p, adding the resulting equations and integrating over Ω, by virtue of the divergence theorem, the boundary conditions (5)1 and (7), it turns out that 1 d(knk2 + kpk2 ) 2 ≤ σ knk2 − kpk2 + M1 C∞ |Ω|(1 + αC∞ ). 2 dt
(16)
Then, setting E = knk2 + kpk2 , it follows that dE ≤ −2E + a ¯. dt
(17)
Following the procedure in [14], one can prove that Σ is an absorbing set. The coexistence equilibria E ∗ = (n∗ , p∗ ) are given by n∗ =
1 1 + αp∗
(18)
and p∗ (positive) solution of f (x) := kα2 x3 + 2αkx2 + k(1 − σα)x + σ(1 − k) = 0. (19) √ −1 + 1 + 3σα 1 < k ≤ 1, σ > there exIf k > 1, E ∗ is unique, while if ασ α ist two coexistence equilibria. In the all other cases, no coexistence equilibria are admissible. Introducing the perturbation fields {X = n − n∗ , Y = p − p∗ }, (4) becomes ∂ X a11 a12 X γ1 ∆X F1 = + + , (20) a21 a22 Y γ2 ∆Y F2 ∂t Y with
σ 1 + 2αp∗ a = − , a = − , 11 12 k(1 + αp∗ ) 1 + αp∗ ∗ a = (1 + αp∗ )p∗ , a = αp , 21 22 1 + αp∗ 2 σX − (1 + αp∗ )XY − αY (XY + n∗ Y + p∗ X), F 1 = − k F2 = (1 + αp∗ )XY + αY (XY + n∗ Y + p∗ X). 6
(21)
To (20) we associate the initial-boundary conditions ( ∇X · n = ∇Y · n = 0, on ∂Ω × R+ . X(x, 0) = X0 (x), Y (x, 0) = Y0 (x), x ∈ Ω.
(22)
Denoting by W ∗ (Ω) the functional space defined by dϕ ∗ 1,2 1,2 + , W (Ω) = ϕ ∈ W (Ω) ∩ W (∂Ω) : = 0 on ∂Ω×R dn our aim is to study the stability of (n∗ , p∗ ) with respect to the perturbations (X, Y ) ∈ [W ∗ (Ω)]2 . 3. Linear stability of the coexistence equilibria Let 0 = α0 < α1 < α2 < · · · < αn < · · ·
(23)
∂X = LX, ∂t
(24)
be the eigenvalues of the operator −∆ on Ω with homogeneous Neumann boundary conditions. The linearized version of (20) is
where L=
a11 + γ1 ∆ a12 a21 a22 + γ2 ∆
,
X=
X Y
.
(25)
For each i = 0, 1, 2, 3, ..., λi is an eigenvalue of L if and only if λi is an eigenvalue of the matrix (see [4]) a11 − αi γ1 a12 ˜ . (26) Li = a21 a22 − αi γ2 The characteristic equation of L˜i is λ2i − I1i λi + I2i = 0, where Iji , (j = 1, 2), are the principal L˜i −invariants given by I1i = trL˜i = I01 − αi (γ1 + γ2 ), σγ − kαp∗ γ1 I2i = det L˜i == γ1 γ2 αi2 + 2 αi + I02 , ∗ k(1 + αp ) 7
(27)
(28)
I0j , (j = 1, 2), being the principal invariants in absence of diffusion, i.e. kαp∗ − σ , I01 = a11 + a22 = k(1 + αp∗ ) (29) σα ∗ 0 ∗ + 1 + 2αp . I2 = a11 a22 − a12 a21 = p − k(1 + αp∗ )2
It is to recall that, in the absence of diffusion, the necessary and sufficient condition guaranteeing the linear stability of E ∗ is {I01 < 0, I02 > 0} [15]. Simple calculation shows that E ∗ is linearly stable in the absence of diffusion if and only if k(1 + αp∗ )2 (1 + 2αp∗ ) kαp∗ < σ < , (30) α holds. Let us underline that, when there exist two coexistence equilibria, (30) implies that only the one with larger predator density is linearly stable. Moreover, (30) is a necessary condition for linear stability in the presence of diffusion too. Theorem 3. If either k(1 + αp∗ )2 (1 + 2αp∗ ) kαp∗ < σ < , α ∗ γ kαp 2 ≥ , γ1 σ
(31)
or
k(1 + αp∗ )2 (1 + 2αp∗ ) kαp∗ < σ < , α 2 ∗ ∗ kαp (kγ1 αp −σγ2 ) γ 2< , < 4kp∗ k(1+2αp∗ )(1+αp∗ )2 −σα , γ1 σ γ1 γ2
(32)
holds, then the coexistence equilibrium (when it exists) is linearly stable. Proof. When (30) holds, I01 < 0 and hence I1i < 0, ∀i. In view of (28)2 , it turns out that either (31) or (32) guarantees that I2i > 0, ∀i. 4. Turing instability In this Section we look for conditions guaranteeing that E ∗ – stable in the absence of diffusion – becomes unstable in the presence of diffusion (Turing 8
instability). Since I01 < 0 ⇒ I1i < 0, ∀i, Turing instability may occur only if I02 > 0 together with I2¯i < 0 for at least one index ¯i. Moreover, in view of (28)2 , a necessary condition for the occurrence of Turing instability is kαp∗ γ2 < . γ1 σ
(33)
It should be noted that in the proposed model, prey is the activator for the predator population (a21 > 0) but it is self-inhibiting (a11 < 0). Predators are inhibitors for the prey (a12 < 0), but they are also self-activating (a22 > 0). In this context, a necessary condition for diffusion-driven instability is that prey disperses faster than predators. In fact, in view of (30), it turns out kαp∗ < 1 and (33) implies necessarily that that σ γ1 > γ2 .
(34)
As illustrated in details by Murray in [16] (see Sect. 2.3, pp. 87-89), biologically there are two types of interaction with consequent qualitatively different reactions: predators spread faster than prey or vice versa. Although a large proportion of the bibliography falls in the first case, the other one is not an uncommon situation and it is precisely the object of the present study. Theorem 4. If (30) and (33) hold together with either αp∗ , γ < 2 α1 (1 + αp∗ ) k(1 + αp∗ )I02 + σα1 γ2 γ1 > kα1 [αp∗ − α1 γ2 (1 + αp∗ )]
(35)
or
(kγ1 αp∗ − σγ2 )2 > 4kp∗ k(1 + 2αp∗ )(1 + αp∗ )2 − σα , γ1 γ2 then Turing instability occurs.
(36)
Proof. If (30) holds, then there is stability in the absence of diffusion. In view of (28)2 , it follows that αp∗ σα1 γ2 + I02 k(1 + αp∗ ) α γ + I21 = α1 γ2 − , (37) 1 1 1 + αp∗ k(1 + αp∗ ) that is negative if (33) and (35), or (33) and (36) hold. 9
10
10
9
0.708
9
1.808
8
0.706
8
1.806
7
0.704
7
1.804
6
0.702
6
1.802
5
0.7
5
1.8
4
0.698
4
1.798
3
0.696
3
1.796
2
0.694
2
1.794
1
0.692
1
1.792
0
0 0
2
4
6
8
10
0
2
4
6
8
10
Figure 1: Plot of the initial state of system (4) in the 2D setting. The initial values are fixed as (n0 , p0 ) = (0.7,1.8) plus a random perturbation of amplitude 0.01. Left panel: initial spatial distribution of prey; right panel: initial spatial distribution of predators.
5. Numerical simulations In this paper, a theoretical analysis of an evolutionary process that involves the interaction of two spatially distributed population, predators and prey, with local diffusion, has been performed. In order to highlight the theoretical results, in this section we provide some numerical simulations which show the impact of different choices for the diffusion coefficients on the dynamics of the two populations. First, we consider a configuration in the parameters space that leads to a single stable coexistence equilibrium in absence of diffusion, according to (30). Specifically, we choose α = 0.3 (weak hunting cooperation among predators), σ = 3 (per capita growth rate of prey) and k = 5 (high predator basic reproduction number for homogeneously mixed populations) and the unique stable coexistence equilibrium is given by E ∗ ≈ (0.66,1.72). In the PDE model (4) we assume an initial condition very close to this equilibrium value plus a small random perturbation. As well known [16], the size of the space domain has to be chosen sufficiently large with respect to the wavelength of the unstable modes, in order to see the instability patterns in the solution. Then we set as space domain the interval [0,10] for 1D simulations, or the square [0,10] × [0,10] for 2D simulations, with a discretization step h = 0.05. Figure 1 shows the chosen initial condition of the system for the prey (left panel) and predator (right panel) populations in the 2D setting. Now, if the diffusion coefficients are chosen, according to Theorem 4, in the range given by (33)-(35), the conditions for instability to emerge are 10
10
10
9
0.6595
1.760
9 1.750
8
8 0.6590
7
7
1.740
0.6585
6
6
5
0.6580
4
1.730
5 4
0.6575
3
1.720
3 0.6570
2 1
0.6565
2
1.710
1 1.700
0
0 0
2
4
6
8
10
0
2
4
6
8
10
Figure 2: Plot of the state of the system for T = 100 in presence of (unbalanced) diffusion: γ1 = 1.5; γ2 = 0.01 with parameters choice α=0.3, σ=3, k=5. Note that both populations (prey, in the left panel; predators in the right panel) start showing spatial patterns.
satisfied. Specifically, we fix γ1 = 1.5 and γ2 = 0.01. Then, the coexistence equilibrium loses stability due to the Turing effect: as shown in Figure 2, spatial patterns emerge in both populations; these patterns persist and increase in amplitude with time (Figure 3). However, it is worth noting that the amplitude of such patterns remains confined within the absorbing set given by (14). This result can be very helpful in analyzing and also forecasting the spatial redistribution of populations in real world applications. 6. Conclusions Ecologists increasingly emphasize how ecosystems reveal spatial self-organization, where large-scale ordered spatial patterns emerge from disordered initial conditions through local interactions. Indeed, organisms concentrate resources in their local environment, enabling persistence even when mean field resource levels become too low for their survival. Ecosystems with regular patterns may be more resilient and resistant to the environmental change compared to homogeneous ecosystems. This process is the key to understand ecological stability and diversity. In view of these real world applications, we have presented a theoretical analysis of an evolutionary process that involves the interaction of two spatially distributed population, predators and prey, with local diffusion. Our analysis, also illustrated by numerical simulations, reveals the impact of diffusion and hunting cooperation on the overall dynamics. We have investigated the stability of the coexistence equilibria and found the conditions guaranteeing the occurrence of Turing instability, lead11
10
10 0.7500
9
6.000
9 0.7000
8 7
0.6500
6
8
5.000
7 4.000
6 0.6000
5
5 3.000
4
0.5500
3
4 3
2.000
0.5000
2
2
1
0.4500
0
1.000
1 0
0
2
4
6
8
10
0
2
4
6
8
10
Figure 3: Plot of the state of the system for T =500 in presence of (unbalanced) diffusion: γ1 = 1.5; γ2 = 0.01 and parameters choice as in the previous Figure. Note that both populations (prey, in the left panel; predators in the right panel) show spatial patterns of great amplitude.
ing to the formation of isolated group. Specifically, we identified a scenario in which weak hunting cooperation allows predators and prey population to spatially redistribute to form isolated (spotted) groups, where predators density rises well above the level sustained by the same prey population in the absence of diffusion. All these results can have important applications for global environmental change, ecosystem adaptation and biological control using predators. Acknowlegements This paper has been performed under the auspices of the G.N.F.M. of INdAM. References [1] M. Alves, F. Hilker, Hunting cooperation and Allee effects in predators, J. Theor. Biol. 419 (2017) 13–22. [2] L. Berec, Impacts of foraging facilitation among predators on predatorprey dynamics, Bull. Math. Biol. 72 (2010) 94–121. [3] D. Boukal, M. Sabelis, L. Berec, How predator functional responses and Allee effects in prey affect the paradox of enrichment and population collapses, Theor. Popul. Biol. 72 (2007) 136–147. 12
[4] R. Cantrell, C. Cosner, Spatial ecology via reaction-diffusion equations, John Wiley and Sons, Ltd, Chichester, UK, 2003. [5] F. Capone, On the dynamics of predator-prey models with the Beddington-De Angelis functional response, under robin boundary conditions, Ricerche Mat. 57 (2008) 137–157. [6] F. Capone, M. Carfora, R. De Luca, I. Torcicollo, On the dynamics of an intraguild predator–prey model, Math. Comput. Simulation 149 (2018) 17–31. [7] F. Capone, V. De Cataldis, R. De Luca, On the stability of a SEIR reaction diffusion model for infections under Neumann boundary conditions, Acta Appl. Math. 132 (2014) 165–314. [8] F. Capone, R. De Luca, On the nonlinear dynamics of an ecoepidemic reaction-diffusion model, Int. J. Nonlinear Mech. 95 (2017) 307–314. [9] B. Datsko, V. Gafiychuk, Complex nonlinear dynamics in subdiffusive activator–inhibitor systems., Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 1673–1680. [10] B. Datsko, Y. Luchko, V. Gafychuk, Pattern formation in fractional reaction difusion systems with multiple homogeneous states, Int. J. Bifurcation Chaos. 22, (2012). 22 (2012) 1250087. [11] R. De Luca, On the long-time dynamics of nonautonomous predatorprey models with mutual interference, Ricerche Mat. 61 (2012) 275–290. [12] M. Du, Z. Wang, H. Hu, Measuring memory with the order of fractional derivative, Scientific Reports 3 (2013) 3431. [13] L. Dung, Dissipativity and global attractors for a class of quasilinear parabolic systems, Commun. Partial Diff. Eqns. 22 (1997) 213–433. [14] J. Flavin, S. Rionero, Qualitative Estimates for Partial Differential Equations: An Introduction, volume 8111, CRC Press Inc., Boca Raton, 1996. [15] D. Merkin, Introduction to the Theory of Stability, volume 24, Text in Applied Mathematic, Springer, 1997. 13
[16] J. Murray, Mathematical Biology II: Spatial Models and Biomedical Applications, 3rd ed., Springer-Verlag, Berlin Heidelberg, 1993. [17] M. Protter, H. Weinberger, Maximum Principles in Differential Equations, Springer, New York, 1999. [18] S. Rionero, A peculiar Lyapunov functional for ternary reaction-diffusion dynamical systems, Boll. U.M.I. 9 (2011) 393–407. [19] S. Rionero, Stability of ternary reaction-diffusion dynamical systems, Rend. Lincei Mat. Appl. 22 (2011) 245–268. [20] S. Rionero, I. Torcicollo, Stability of a continuous reaction-diffusion Cournot-Kopel duopoly game model, Acta Appl. Math. 132 (2014) 505– 513. [21] S. Rionero, I. Torcicollo, On the dynamics of a nonlinear reaction– diffusion duopoly model, Int. J. Non-Linear Mech. 99 (2018) 105–111. [22] X. Tang, Y. Song, T. Zhang, Turing-Hopf bifurcation analysis of a predator-prey model with herd behavior and cross-diffusion, Nonlinear Dyn. 86 (2016) 73–89. [23] A. Terry, Predator-prey models with component Allee effect for predator reproduction, J. Math. Biol. 71 (2015) 1325–1352. [24] I. Torcicollo, On the non-linear stability of a continuous duopoly model with constant conjectural variation, Int. J Non-Linear Mech. 81 (2016) 268–273. [25] A. Turing, The Chemical Basis for Morphogenesis, Philos. Trans. R. Soc. Lond. Ser. B, Biol. Sci. 237 (1952) 37–72. [26] H. Yin, X. Wen, Pattern formation through temporal fractional derivatives, Scientific Reports 8 (2018) 5070. [27] H. Yin, X. Xiao, X. Wen, Turing patterns in a predator-prey system with self-diffusion, Abstract and Applied Analysis 2013 (2013). [28] J. Zhang, W. Li, Y. Wang, Turing patterns of a strongly coupled predator-prey system with diffusion effects, Nonlinear Anal. 74 (2011) 847–858. 14