# Dynamics of a Three-Patch Prey-Predator System with the Impact of Dispersal Speed Incorporating Strong Allee Effect on Double Prey

Authors: S. Biswas; D. Pal; P. K. Santra; Ebenezer Bonyah; G. S. Mahapatra
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

This paper presents the dispersal effect on a prey-predator model with three-patch states incorporating a strong Allee effect on the first two prey populations. The prey species are considered to be mutable to exhibit a balanced dispersal between the patches. The dispersal among patches is directed through lower fit patches to higher fit patches. This paper derives a new approach for dispersal and Allee’s effect with the specific condition on the stability on a three-patch of a three-species prey-predator anticipated system. The persistence of the system is observed because of the dispersal effect on the three-patch prey-predator system. The numerical simulations of analytical findings are presented using hypothetical parameter values to relate real-world prey-predator situations for balancing ecology.


---

## Body

## 1. Introduction

In mathematical ecology, most of the researchers considered the prey-predator system in a homogeneous environment. But in a real life situation, the environment is heterogeneous, which contains different patches connected through migration. Many ecologists and researchers ([1–5]) have studied the impact of predator species migration on prey-predator interactions. In their work, the prey density was significant and the predator species was considered to remain in the specified patch. Many researchers comprehensively studied the dispersal model in a multipatch environment ([6–20]). Amarasekare [13] studied two-patch models of single species with local density-dependent dispersal and spatial heterogeneity. Padrón and Trevisan [18] considered a single species’ logistic growth model composing several habitats connected by linear migration. Stephen Cantrell et al. [19] examined evolutionary stability strategies for dispersal in heterogeneous patchy environments. In most cases, the randomness of the dispersal rate between different patches was assumed to be fixed. Kang et al. [20] formulated a Rosenzweig-MacArthur prey-predator model in two patch environments. The rate of dispersal has a great influence on stable prototypes and the persistence of the species shown in various research works ([21–24]). Ruxton [22] investigated the stability behaviour of a population model by adding density-dependent migration between nearest-neighbour populations. Rohani and Ruxton [23] presented a two-species interaction to establish the non-stabilization of the system for density-independent movement between populations. Fretwell and Lucas [25] proposed ideal free distribution, i.e. individuals in different patches possess similar fitness. Several researchers have already discussed this concept after that ([10, 26, 27]).In balance dispersion [10], the dispersion rate remains constant when the population is in the equilibrium state, providing ideal free distribution. But constant rates of dispersal for the passive animal can lead the prey-predator system to a stable situation ([10, 28]). Furthermore, passive dispersal may create negative density dependence staffing rates of the population to stabilize the prey-predator system at an equilibrium state. Therefore, due to heterogeneity in patch or dispersal rates, an unstable population may behave like a stable population. Dispersal rates cannot be considered very high because they may harmonize local stability behaviour between patches. Therefore, to show the effect of dispersal, researchers have developed their models with the help of fixed dispersal rates. The population density alters very fast because of the dispersal process in most of their models. Krivan and Sirot [27] observed that, the stable population dynamics for two competitive species becomes unstable due to rapid adaptive animal dispersal. Abrams et al. ([29, 30]) also observed the profound impact of the dispersal process on the population dynamics.Modeling the prey-predator system incorporating the dispersal concept is an active field of research. The dispersal rate of prey species steadily increases with the enlargement of predators in a patch ([31, 32]). Another exciting matter in this field is the Allee effect, a mechanism where the prey population shows a negative growth rate at low density ([33, 34]). There are several reasons to consider the Allee effect: mate restrictions, stumpy probability of victorious meeting [35], food mistreatment, supportive defence, inbreeding dejection [36], etc. Predator evasion [37] for evolutionary alterations is another important cause to consider the Allee effect. The strong Allee effect ([38–43]) has been considered here, whereas the concept of the weak Allee effect ([44, 45]) also exists in this regard. This classification depends on the per capita enlargement rate of the population at low density. The mathematical formulation of the growth equation due to the Allee effect takes the form Ψy=ry1−y/ly/l0−1 where, y represents the population size at time t, r is the per capita growth rate, environmental carrying capacity is presented by l and l0 is the Allee edge. Here, 0<l0⋘. The population change rate continuously diminishes ([46, 47]) and finally departs to extermination ([42, 43]) if the relation 0<y<l0 holds. On the contrary, at low population density, the population intensification rate decreases but stays positive in the weak Allee effect. Researchers ([48–51]) observed the impact of the Allee effect on the dynamical behaviour of the prey-predator system. Celik and Duman [52] enormously studied the stability behaviour of a prey-predator model with and devoid of the Allee effect. The impact of the Allee effect on a prey-predator harvesting system was investigated by Javidi and Nyamorady [53]. Wang et al. [54] discussed the prey-predator system that exhibits intricate population dynamics due to the Allee effect. Zhou et al. [55] studied the population dynamics of a prey-predator system with the Allee effect. Kent et al. [56] used the concept of Allee effect as prey outfluxes. Therefore, numerous works ([57–60]) on prey-predator interactions have already been presented in this field of ecological modeling. A small number of researchers studied the dispersal dynamics of the prey-predator system in a patchy environment with a strong Allee effect.Recently, Pal and Samanta [61] presented the dispersal dynamics of the Allee effect in one prey of a prey-predator system in two patch environment. Saha and Samanta [62] studied the dispersal dynamics of a prey-predator system in two patch environments with two prey species considered the Allee effect. Consequently, we are introducing a prey-predator system with dispersal and strong Allee effect in three-patch environments, namely Patch 1, Patch 2 and Patch 3. Every patch consists of a pair of prey-predator species and a strong Allee effect in the prey population escalation in the first two patches. Our main objective is to study the impact of dispersal speed and the Allee effect on the population dynamics of our desired system. Ecology always demands a balance among its inhabitants, so stability of the population is always expected in this regard. We have studied the effect of dispersal and the Allee effect on the persistence of species. So, we hope that this study may help ecologists in their development of the environment.The next section reflects the mathematical structure of the prey-predator interaction in three-patch environment. Section3 is equipped with the positivity and boundedness characteristics for our proposed prey-predator system. Equilibrium points and their stability conditions have been discussed in Section 4. Bifurcation behaviour of the proposed model has been analysed in Section 5. Numerical verifications have been done in Section 6 and the concluding remarks have been given in Section 7.


## 2. Formation of Three-Patch Prey-Predator System

We opt for a six species (three each prey and predators) prey-predator system in three-patchi=1,2,3 environments based on the following notations and assumptions.Notations:xi: prey population density in patch i.yi: predator population density in patch i.ri: growth rate of prey species in patch i.ki: environmental carrying capacity of the prey in ith patch in nonattendance of predation and dispersal.k¯i: threshold value of the Allee effect in patch i without predation as well as dispersal 0<k¯i≪ki.λi: predation rate in patch i.δ: dispersal speed among patches.dij: probability of dispersion from the ith patch to the jth patch j=1,2,3;i≠j.ei: conversion rate of prey biomass to predator biomass in patch i.mi: mortality rate of the predator species in patch i.Assumptions:(i)
Prey population growth rate is affected by the Allee effect in the first two patches.(ii)
All three predator species are free from the Allee effect.(iii)
Prey species are movable to higher fitness patches.(iv)
Conversion rate of prey biomass to predator biomass less than the predation rateei≤λi∀i.(v)
We consider balanced dispersal ([11, 63]), mathematically d12/d21=k2/k1, d13/d31=k3/k1, d23/d32=k3/k2 i.e, there is no mesh progress among the patches.The graphical view of our proposed system of three-patch environment system is presented in Figure1.Figure 1
Flow diagram of the three-patch prey-predator interaction.Based on the above notations, assumptions, and flow diagram, the population-dispersal dynamics can be presented by the following set of nonlinear differential equations:(1)dx1dt=r1x11−x1k1x1k¯1−1−λ1x1y1+δd21x2+d31x3−d12x1−d13x1,dx2dt=r2x21−x2k2x2k¯2−1−λ2x2y2+δd12x1+d32x3−d21x2−d23x2,dx3dt=r3x31−x3k3−λ3x3y3+δd13x1+d23x2−d31x3−d32x3.dy1dt=e1x1y1−m1y1,dy2dt=e2x2y2−m2y2,dy3dt=e3x3y3−m3y3.with primary information xi0>0 and yi0>0i=1,2,3.Since we considered the balanced dispersal, and hence mathematically, we can writed12/d21=k2/k1, d13/d31=k3/k1, d23/d32=k3/k2. Therefore, the carrying capacities, which are identical to the population abundances in all patches, always imply no mesh progress among the patches. Consequently, this situation is precisely the same as without dispersal. Therefore, without dispersal speed, the proposed model equation (1) converts to the following system of equations:(2)dx1dt=r1x11−x1k1x1k¯1−1−λ1x1y1≡x1f1,dx2dt=r2x21−x2k2x2k¯2−1−λ2x2y2≡x2f2,dx3dt=r3x31−x3k3−λ3x3y3≡x3f3,dy1dt=e1x1y1−m1y1≡y1f4,dy2dt=e2x2y2−m2y2≡y2f5,dy3dt=e3x3y3−m3y3≡y3f6.with initial densities xi0>0 and yi0>0 for all i=1,2,3.Whenδ=0, the population in each patch progresses separately. In that case, the non-trivial equilibrium can be found by solving the system of equations fi=0i=1,2,3,4,5,6. Our primary target is to find the interior equilibrium point x1,x2,x3,y1,y2,y3 of the proposed system equation (1) with respect to k1, k2, k3, m1, m2, m3 and we shall also observe the impact of the dispersal rate δ on the said equilibrium point.


## 3. Positivity and Boundedness of the Proposed System

Let us consider the following theorems to ensure that the anticipated model equation (1) is well-posed.Theorem 1.
Each solution of the anticipated system equation (1) starting from ℝ+6 stay positive for all time.Proof.
The right-hand side of the model equation (1) is locally Lipschitzian in the space of continuous functions.
Therefore, the solutionx1t,x2t,x3t,y1t,y2t,y3t along with the positive initial solution x10,x20,x30,y10,y20,y30 exists and is unique in the interval 0,θ where 0<θ<+∞.
Firstly we prove thatx1t>0 for all 0≤t<θ.
If it is false, then there existT1∈0,θ so that x1T1=0, dx1/dt<0 as well as x1t>0 for all t∈0,T1.
Now, the first equation of the system equation (1) provides,(3)dx1T1dt=r1x1T11−x1T1k1x1T1k¯1−1−λ1x1T1y1T1+δd21x2T1+d31x3T1−d12x1T1−d13x1T1=δd21x2T1+d31x3T1>0.which contradicts dx1T1/dt<0 and hence x1t>0 for all t∈0,θ.
Next, we claim thatx2t>0 for all 0≤t<θ. If it is false, then there exist T2∈0,θ so that x2T2=0, dx2/dt≤0 as well as x2t>0 for all 0≤t<T2.
Now, the second equation of the system equation (1) gives(4)dx2T2dt=r2x2T21−x2T2k2x2T2k¯2−1−λ2x2T2y2T2+δd12x1T2+d32x3T2−d21x2T2−d23x2T2=δd12x1T2+d32x3T2>0.which again contradicts dx2/dt≤0 and consequently, x2t>0 for all t∈0,θ.
We further claim thatx3t>0 for all 0≤t<θ. If our claim is false, then their exist T3∈0,θ so that x3T3=0, dx3/dt≤0 as well as x3t>0 for all 0≤t<T3.
Now, from the third equation of equation (1), we get(5)dx3T3dt=r3x3T31−x3T3k3−λ3x3T3y3T3+δd13x1T3+d23x2T3−d31x3T3−d32x3T3=δd13x1T3+d23x2T3>0.this also contradicts dx3/dt≤0. Hence, x3t>0 for all t∈0,θ.
Again, from the last three equations of the system equation (1) we can obtain,(6)yit=yi0exp∫0teixis−mids>0,i=1,2,3.therefore, from the above discussion, we can conclude xit>0 as well as yit>0i=1,2,3 for all t≥0. Hence, the theorem is proved.Theorem 2.
For balanced dispersal, every solution of the anticipated system equation (1) which starts from ℝ+6 is uniformly bounded there.Proof.
To prove this theorem, we have the following cases.Case 1.
x10≤k1, x20≤k2, x30≤k3.
Forx10≤k1, let us claim that x1t≤k1. Assume that the claim is false. Therefore, there exist t1 and t2 in such a way that x1t1=k1 as well as x1t>k1 for all t∈t1,t2 where x2t<k2 and x3t<k3 for all t∈t1,t2 as k2, k3 are the respective carrying capacities of patch 2 and patch 3.
Now, for allt∈t1,t2 we have,(7)x1t=x10exp∫0tϕx1s,x2s,x3s,y1s,y2s,y3sds.where,(8)ϕx1s,x2s,x3s,y1s,y2s,y3s=r11−x1sk1x1sk¯1−1−λ1y1s+δd21x2sx1s+d31x3sx1s−d12−d13.
Now,(9)ϕ=r11−x1k1x1k¯1−1−λ1y1+δd12d21d12x2x1−1+δd13d31d13x3x1−1=r11−x1k1x1k¯1−1−λ1y1+δd12k1k2x2x1−1+δd13k1k3x3x1−1asd21d12=k1k2,d31d13=k1k3<r11−x1k1x1k¯1−1−λ1y1+δd12x2k2−1+δd13x3k3−1as x1t>k1.
Therefore,(10)x1t=x10exp∫0t1ϕx1s,x2s,x3s,y1s,y2s,y3sds×exp∫t1tϕx1s,x2s,x3s,y1s,y2s,y3sds=x1t1exp∫t1tϕx1s,x2s,x3s,y1s,y2s,y3sds=k1exp∫t1tϕx1s,x2s,x3s,y1s,y2s,y3sdsasx1t1=k1.
But for allt∈t1,t2, ϕx1,x2,x3,y1,y2,y3<0 as x1t>k1, 0<k¯1<k1 and x2<k2, x3<k3.
Thus,x1t<k1 which is a contradiction and proves that the claim is valid. Hence x1t≤k1.
Similarly, forx20≤k2 we can prove that x2t≤k2 taking d12/d21=k2/k1 and d32/d23=k2/k3.
Forx30≤k3, we claim x3t≤k3. Again, if we assume that the claim is not true, then there exist t′, t″ such that x3t′=k3 and x3t>k3 for all t∈t′,t″ where x1t<k1 and x2t<k2 as k1, k2 are the respective carrying capacities of patch 1 and patch 2. Now, for all t∈t′,t″ we have,(11)x3t=x30exp∫0tΨx1u,x2u,x3u,y1u,y2u,y3udu.where,(12)Ψx1u,x2u,x3u,y1u,y2u,y3u=r31−x3uk3−λ3y3u+δd13x1ux3u+d23x2ux3u−d31−d32.
Now,(13)Ψ=r31−x3k3−λ3y3+δd31d13d31x1x3−1+δd32d23d32x2x3−1=r31−x3k3−λ3y3+δd31k3k1x1x3−1+δd32k3k2x2x3−1asd13d31=k3k1,d23d32=k3k2<r31−x3k3−λ3y3+δd31x1k1−1+δd32x2k2−1asx3t>k3.
Then,(14)x3t=x30exp∫0t′Ψx1u,x2u,x3u,y1u,y2u,y3udu×exp∫t′tΨx1u,x2u,x3u,y1u,y2u,y3udu=x3t′exp∫t′tΨx1u,x2u,x3u,y1u,y2u,y3udu=k3exp∫t1tΨx1u,x2u,x3u,y1u,y2u,y3uduasx3t′=k3.
Butx3>k3 implies Ψx1u,x2u,x3u,y1u,y2u,y3u<0 as x1<k1 and x2<k2. Hence x3t<k3 which contradicts our assumption. Consequently, our claim is true and x3t≤k3.
Therefore, we getx1t≤k1, x2t≤k2, x3t≤k3.Case 2.
x10≤k1, x20≤k2, x30>k3.
Forx10≤k1 and x20≤k2 we can easily prove that x1t≤k1 and x2t≤k2 (as in Case 1).
Forx30>k3, we claim x3t≤k3. If we assume that the claim is not true, then there exist t′,t″ such that x3t′=k3 and x3t>k3 for all t∈t′,t″ where x1t<k1 and x2t<k2 as k1, k2 are the respective carrying capacities of patch 1 and patch 2. Now, for all t∈t′,t″ we have,(15)x3t=x30exp∫0tΦx1u,x2u,x3u,y1u,y2u,y3udu.where,(16)Φx1u,x2u,x3u,y1u,y2u,y3u=r31−x3uk3−λ3y3u+δd13x1ux3u+d23x2ux3u−d31−d32.
Now,(17)Φ=r31−x3k3−λ3y3+δd31d13d31x1x3−1+δd32d23d32x2x3−1=r31−x3k3−λ3y3+δd31k3k1x1x3−1+δd32k3k2x2x3−1asd13d31=k3k1,d23d32=k3k2<r31−x3k3−λ3y3+δd31x1k1−1+δd32x2k2−1asx3t>k3.
Then,(18)x3t=x30exp∫0tΦΦx1u,x2u,x3u,y1u,y2u,y3udu×exp∫t1tΦx1u,x2u,x3u,y1u,y2u,y3udu=x3t′exp∫t′tΦx1u,x2u,x3u,y1u,y2u,y3udu=k3exp∫t1tΦx1u,x2u,x3u,y1u,y2u,y3uduasx3t′=k3.
Butx3>k3 implies Φx1u,x2u,x3u,y1u,y2u,y3u<0 as x1<k1 and x2<k2. Hence x3t<k3 which contradicts our assumption. Consequently, our claim is true and x3t≤k3.
Therefore, we getx1t≤k1, x2t≤k2, x3t≤k3.Case 3.
x10≤k1, x20>k2, x30≤k3.
Forx10≤k1 and x30≤k3 we have x1t≤k1 and x2t≤k2 (as proved in Case 1).
Forx20>k2 let us claim that x2t≤k2. Assume that the claim is false. Therefore, there exist t1 and t2 in such a way that x2t1=k2 as well as x2t>k2 for all t∈t1,t2 where x1t<k1 and x3t<k3 for all t∈t1,t2 as k1, k3 are the respective carrying capacities of patch 1 and patch 3.
Now, for allt∈t1,t2 we have,(19)x2t=x20exp∫0tϕx1s,x2s,x3s,y1s,y2s,y3sds,where,(20)ϕx1s,x2s,x3s,y1s,y2s,y3s=r21−x2sk2x2sk¯2−1−λ2y2s+δd12x1sx2s+d32x3sx2s−d21−d23.
Now,(21)ϕ=r21−x2k2x2k¯2−1−λ2y2+δd21d12d21x1x2−1+δd23d32d23x3x2−1=r21−x2k2x2k¯2−1−λ2y2+δd21k2k1x1x2−1+δd23k2k3x3x2−1asd12d21=k2k1,d32d23=k2k3<r21−x2k2x2k¯2−1−λ2y2+δd21x1k1−1+δd23x3k3−1asx2t>k2.
Therefore,(22)x2t=x20exp∫0t1ϕx1s,x2s,x3s,y1s,y2s,y3sds×exp∫t1tϕx1s,x2s,x3s,y1s,y2s,y3sds=x2t1exp∫t1tϕx1s,x2s,x3s,y1s,y2s,y3sds=k2exp∫t1tϕx1s,x2s,x3s,y1s,y2s,y3sdsasx2t1=k2.
But for allt∈t1,t2,ϕx1,x2,x3,y1,y2,y3<0[ as x2t>k2,0<k¯2<k2 and x1<k1, ]x3<k3.
Thus,x2t<k2 which is a contradiction and proves that the claim is valid. Hence x2t≤k2.
Therefore, in this case we also getx1t≤k1, x2t≤k2, x3t≤k3.Case 4.
x10>k1, x20≤k2, x30≤k3.
Forx10>k1 we can prove that x1t≤k1 using the method of Case 3. Again for x20≤k2 and x30≤k3 we get x2t≤k2 and x3t≤k3 (from Case 1).
Therefore, we getx1t≤k1,x2t≤k2,x3t≤k3.Case 5.
x10≤k1, x20>k2, x30>k3.
Forx10≤k1 we have x1t≤k1 (from Case 1) and we get x2t≤k2 for x20≤k2 (from Case 3). Again x3t≤k3 gives x30≤k3 (from Case 2).
Therefore, we getx1t≤k1, x2t≤k2, x3t≤k3.Case 6.
x10>k1, x20≤k2, x30>k3.
Forx10>k1 we can prove that x1t≤k1 (as in Case 3). In a similar way, considering Case 1, we get x2t≤k2 for x20≤k2. Lastly, for x30≤k3 we can prove x3t≤k3 (from Case 2).
Consequently, we getx1t≤k1, x2t≤k2,x3t≤k3.Case 7.
x10>k1, x20>k2, x30≤k3.
Forx10>k1 and x20>k2, we can prove x1t≤k1 and x2t≤k2 (using Case 3). Again, for x30≤k3, we get x3t≤k3 (using Case 1).
Therefore, we getx1t≤k1,x2t≤k2,x3t≤k3.Case 8.
x10>k1, x20>k2, x30>k3.
Here we claimlimsupt⟶∞x1t≤k1,limsupt⟶∞x2t≤k2,limsupt⟶∞x3t≤k3.
Let us assume thatlimsupt⟶∞x1t≰k1. Then x1t>k1 for all t>0.
Now taking−λ1y1t+δd21x2t/x1t+d31x3t/x1t−d12−d13.
ϕx1t,x2t,x3t,y1t,y2t,y3t=r11−x1t/k1x1t/k¯1−1 we have(23)ϕ=r11−x1k1x1k¯1−1−λ1y1+δd12d21d12x2x1−1+δd13d31d13x3x1−1=r11−x1k1x1k¯1−1−λ1y1+δd12k1k2x2x1−1+δd13k1k3x3x1−1asd21d12=k1k2,d31d13=k1k3<r11−x1k1x1k¯1−1−λ1y1+δd12x2k2−1+δd13x3k3−1asx1t>k1.
Now,k2, k3 are the respective carrying capacities of patch 2 and patch 3. So, x2t<k2 and x3t<k3 (carrying capacity means the maximum population that can be achieved by the environment). Thus, ϕx1t,x2t,x3t,y1t,y2t,y3t<0 as x1t>k1≫k¯1.
Again, the first equation of (1) gives,(24)x1t=x10exp∫0tϕx1s,x2s,x3s,y1s,y2s,y3sds<x10.
Also, we have,dx1/dt<x1′/k¯1−1r1x11−x1/k1 where x1′=limsupt⟶∞x1t.
This implies thatlimsupt⟶∞x1t≤k1 which contradicts the assumption and so limsupt⟶∞x1t≤k1.
Similarly, we can prove that,limsupt⟶∞x2t≤k2.
Lastly let us assume thatlimsupt⟶∞x3t≰k3. Then x3t>k3 for all t>0.
Now takingΨx1t,x2t,x3t,y1t,y2t,y3t=, r31−x3t/k3−λ3y3t+δd13x1t/x3t+d23x2t/x3t−d31−d32 we have,(25)Ψ=r31−x3k3−λ3y3+δd31d13d31x1x3−1+δd32d23d32x2x3−1=r31−x3k3−λ3y3+δd31k3k1x1x3−1+δd32k3k2x2x3−1asd13d31=k3k1,d23d32=k3k2<r31−x3k3−λ3y3+δd31x1k1−1+δd32x2k2−1asx3t>k3.
Now,k1, k2 are the respective carrying capacities of patch 1 and patch 2. So, x1t<k1 and x2t<k2 (carrying capacity means the maximum population that can be achieved by the environment). Thus Ψx1t,x2t,x3t,y1t,y2t,y3t<0 as x3t>k3.
Again, the third equation of (1) gives,(26)x3t=x30exp∫0tΨx1u,x2u,x3u,y1u,y2u,y3udu<x30.
We also have,dx3/dt<r3x31−x3/k3 from the third equation of (1). So, limsupt⟶∞x3t≤k3. This again contradicts the assumption and we get limsupt⟶∞x3t≤k3.
Therefore, we concludelimsupt⟶∞x1t≤k1,limsupt⟶∞x2t≤k2,limsupt⟶∞x3t≤k3.
LetU=∑i=13Axit+Biyit where Bi=Aλi/ei for i=1,2,3. Also let γi=riki+k¯i+ri2ki2+k¯i2−kik¯i+3rikik¯i/3ri for i=1,2; γ3=k3/2r31+r3 and A=1/γ3−r3/k3γ3+r3+1+∑i=12γi−ri/kik¯iγi2+ri1/ki+1/k¯iγi−ri+1.
Then,(27)dUdt=∑i=13Adxidt+Bidyidt=A∑i=12rixi1−xikixik¯i−1+xi+Ar3x31−x3k3+x3−Aλ1m1y1e1−Aλ2m2y2e2−Aλ3m3y3e3−A∑i=13xi=Ar3x31−x3k3+x3+∑i=12rixi1−xikixik¯i−1+xi−∑i=13Axi+Bimiyi.
Now, if we takeH=r3x31−x3/k3+x3+∑i=12rixi1−xi/kixi/k¯i−1+xi, then the value of H is maximum at γ1,γ2,γ3 and Hmax=1/A.
Therefore,(28)dUdt≤AHmax−μUwhereμ=min1,m1,m2,m3
This implies,(29)dUdt≤1−μU
Now applying the theory of differential inequality, we get,(30)Ut≤1μ1−e−μt+Ux10,x20,x30,y10,y20,y30e−μt.
Takingt⟶∞,(31)Ut≤1μ
Then all the solutions of (1) that initiate in ℝ+6 are confined in the region:(32)Γ=x1,x2,x3,y1,y2,y3:xi∈0,ki,i=1,2,3;0≤Ut≤1μ+ε,ε>0.
This proves the theorem.


## 4. Equilibria and Stability of the Proposed Model

For the non-trivial equilibrium point of the proposed model (1) in the absence and presence of dispersal, we have δ=0 and δ≠0 for the respective cases.Case 9.
δ=0.
Clearly, the non-trivial equilibriumE1x10∗,x20∗,x30∗,y10∗,y20∗,y30∗ is a positive solution of the following system of equations:(33)r11−x1k1x1k¯1−1−λ1y1=0,(34)r21−x2k2x2k¯2−1−λ2y2=0(35)r31−x3k3−λ3y3=0(36)e1x1−m1=0(37)e2x2−m2=0(38)e3x3−m3=0
From equations (37), and (39) we have xi0∗=mi/ei>0i=1,2,3. Also, from (34) we obtain,(39)yi0∗=riλi−xi2kik¯i+1ki+1k¯ixi−1xi=miei,i=1,2=riλikik¯i−mi2ei2+ki+k¯imiei−kik¯i,i=1,2.
Thus,yi0∗>0i=1,2 if(40)ki+k¯imiei−kik¯i>mi2ei2,i=1,2.
Now, let us take that,(41)ki+k¯i>2miei,i=1,2.
Also, letk1=m1/e1+η1, k¯1=m1/e1−η2, k2=m2/e2+η3 and k¯2=m2/e2−η4 where η1>η2; η2<m1/e1,η3>η4; η4<m2/e2. Then (41) holds for both values.(42)Now,k1+k¯1m1e1−k1k¯1−m12e12=2m12e12+η1−η2m1e1−m1e1+η1m1e1−η2−m12e12=η1η2>0,k2+k¯2m2e2−k2k¯2−m22e22=2m22e22+η3−η4m2e2−m2e2+η3m2e2−η4−m22e22=η3η4>0.
Therefore, the condition (40) is verified and we get yi0∗>0, i=1,2.
Again from (36) we get,(43)y30∗=r3λ31−x3k3x3=m3e3=r3λ3k3k3−m3e3.
Therefore,y30∗>0 if(44)k3−m3e3>0⇒k3>m3e3.
Thus, we conclude(45)xi0∗=mieii=1,2,3,yi0∗=riλikik¯i−mi2ei2+ki+k¯imiei−kik¯iforki+k¯i>2mieii=1,2,y30∗=r3λ3k3k3−m3e3fork3>m3e3.Case 10.
δ≠0.
In this case, we have to consider the following set of equations to evaluate the interior equilibriumE2x1δ∗,x2δ∗,x3δ∗,y1δ∗,y2δ∗,y3δ∗(46)r1x11−x1k1x1k¯1−1−λ1x1y1+δd21x2+d31x3−d12x1−d13x1=0,(47)r2x21−x2k2x2k¯2−1−λ2x2y2+δd12x1+d32x3−d21x2−d23x2=0,(48)r3x31−x3k3−λ3x3y3+δd13x1+d23x2−d31x3−d32x3=0,(49)e1x1−m1=0,(50)e2x2−m2=0,(51)e3x3−m3=0,
The positive solution of equations (46)–(51) is the interior equilibrium E2x1δ∗,x2δ∗,x3δ∗,y1δ∗,y2δ∗,y3δ∗ of the proposed model (1) with dispersal (δ≠0). Now solving (46)–(51), we get the interior equilibrium E2x1δ∗,x2δ∗,x3δ∗,y1δ∗,y2δ∗,y3δ∗ as(52)xiδ∗=mieii=1,2,3,y1δ∗=r1λ11−m1e1k1m1e1k¯1−1+δλ1e1d21m2e2m1+e1d31m3e3m1−d12−d13=r1m1e2e3e1k1−m1m1−e1k¯1+δk1k¯1e12m2e1e3d21+m3e1e2d31−m1e2e3d12+d13λ1m1k1k¯1e12e2e3,y2δ∗=r2λ21−m2e2k2m2e2k¯2−1+δλ2e2d12m1e1m2+e2d32m3e3m2−d21−d23=r2m2e1e3e2k2−m2m2−e2k¯2+δk2k¯2e22m1e2e3d12+m3e1e2d32−m2e1e3d21+d23λ2m2k2k¯2e1e22e3,y3δ∗=r3λ31−m3e3k3+δλ3e3d13m1e1m3+e3d23m2e2m3−d31−d32=r3m3e1e2e3k3−m3+δk3e3m1e2e3d13+m2e1e3d23−e1e2m3d31+d32λ3m3k3e1e2e3,where x1δ∗, x2δ∗, x3δ∗ , y1δ∗, y2δ∗ and y3δ∗ are all considered positive.


### 4.1. Stability Condition of the Proposed System without Dispersal

At this pointE1x10∗,x20∗,x30∗,y10∗,y20∗,y30∗, the variation matrix is provided by,(53)VE1=−r1x10∗k1k¯12x10∗−k1−k¯100−λ1x10∗000−r2x20∗k2k¯22x20∗−k2−k¯200−λ2x20∗000−r3x30∗k300−λ3x30∗e1y10∗000000e2y20∗000000e3y30∗000,Theorem 3.
The system (1) is unstable without dispersal δ=0 under the condition ki+k¯i>2mi/ei.Proof.
The characteristic equation of the variation matrixVE1 is given by,(54)Δ=−r1k1k¯1x10∗2x10∗−k1−k¯1−λ00−λ1x10∗000−r2x20∗k2k¯22x20∗−k2−k¯2−λ00−λ2x20∗000−r3x30∗k3−λ00−λ3x30∗e1y10∗00−λ000e2y20∗00−λ000e3y30∗00−λ=0.
Evaluating, we get (see Appendix A),(55)cλ2+r3x30∗k3λ+e3λ3x30∗y30∗λ2+r2x20∗k2k¯22x20∗−k2−k¯2λ+e2λ2x20∗y20∗,λ2+r1x10∗k1k¯12x10∗−k1−k¯1λ+e1λ1x10∗y10∗=0.
Therefore, we have,(56)λ2+r3x30∗k3λ+e3λ3x30∗y30∗=0,(57)λ2+r2x20∗k2k¯22x20∗−k2−k¯2λ+e2λ2x20∗y20∗=0,(58)λ2+r1x10∗k1k¯12x10∗−k1−k¯1λ+e1λ1x10∗y10∗=0,
Now, we observe thatr3x30∗/k3>0 and e3λ3x30∗y30∗>0. Therefore, the roots of (56) are either all negative or have negative real parts.
Again2xi0∗−ki−k¯i<0⇔ki+k¯i>2mi/ei, i=1,2.
Then all the roots of (57) and (58) have positive real parts. Thus, the system (1) is unstable without dispersal δ=0 with the assumption ki+k¯i>2mi/ei, i=1,2.


### 4.2. Persistence of the Proposed System in the Presence of Dispersal

The biological meaning of persistence is the survival of all populations in future time. For the Kolgomorov type equations, persistence was discussed by Freedman and Waltman [64]. Mathematical definition for persistence of a system can be written as follows:Definition 1.
(Persistence): ifx0>0 and limt⟶∞infxt>0, xt persists. The system is said to persist if each component of the system persists.Theorem 4.
The system (1) persists.Proof.
Consider the matrix:(59)Af=D+diagf10,f20,f30,f40,f50,f60,where, D=dij is a matrix whose elements are given by,(60)dij=Dijfor i≠j,−∑k=1nDikαikfor i=j.
HereDij denotes the nonnegative diffusion coefficient of the species from the j-th patch to the i-th patch for i≠j and Dii=0, i=1,2,3,4,5,6. For the continuous diffusion case, the parameter αij≥0 corresponds to the boundary conditions. Now we have,(61)f1=r11−x1k1x1k¯1−1−λ1y1,f2=r21−x2k2x2k¯2−1−λ2y2,f3=r31−x3k3−λ3y3,f4=e1x1−m1,f5=e2x2−m2,f6=e3x3−m3.
Here we can recall the theorem which states that,sAf>0 implies liminft≥0xit≥δ>0 and liminft≥0yit≥δ>0 for i=1,2,3 where, sAf=maxReλ:λ∈the set of eigenvalues ofAf. Interested readers can follow reference [65] about the proof of this theorem. The characteristic equation for the matrix Af is,(62)detAf−λI=−r1−λ000000−r2−λ000000r3−λ000000−δd12+d13−m1−λδd12δd13000δd21−δd21+d23−m2−λδd23000δd31δd32−δd31+d32−m3−λ=0.
Which implies, (see Appendix B)(63)lλ+r1λ+r2λ−r3λ3+P+Q+Rλ2+PQ+QR+RP−δ2d12d21+d13d31+d23d32λ+PQR−δ2Pd23d32+Qd13d31+Rd12d21−δ3d12d23d31+d13d21d32=0.where, P=δd12+d13+m1, Q=δd21+d23+m2 and R=δd31+d32+m3. Clearly, one of the roots of (62) is positive. Therefore sAf>0. Consequently, the system (1) persists.


### 4.3. Stability in the Presence of Dispersal

Let us consider the understated theorem for the stability of the equilibriumE2x1δ∗,x2δ∗,x3δ∗,y1δ∗,y2δ∗,y3δ∗.Theorem 5.
The equilibriumE2x1δ∗,x2δ∗,x3δ∗,y1δ∗,y2δ∗,y3δ∗ is stable under some conditions.Proof.
The characteristic equation forVE2 can be written as,(64)r1k1k¯1x1δ∗2k1+2k¯1−3r1x1δ∗2k1k¯1−r1−λ1y1δ∗−δd12+d13−λδd21δd31−λ1x1δ∗00δd12r2k2k¯2x2δ∗2k2+2k¯2−3r2x2δ∗2k2k¯2−r2−λ2y2δ∗−δd21+d23−λδd320−λ2x2δ∗0δd13δd23r3−2r3k3x3δ∗−λ3y3δ∗−δd31+d32−λ00−λ3x3δ∗e1y1δ∗00e1x1δ∗−m1−λ000e2y2δ∗00e2x2δ∗−m2−λ000e3y3δ∗00e3x3δ∗−m3−λ=0.
Taking(65)M1=r1k1k¯1x1δ∗2k1+2k¯1−3r1x1δ∗2k1k¯1−r1−λ1y1δ∗−δd12+d13,M2=r2k2k¯2x2δ∗2k2+2k¯2−3r2x2δ∗2k2k¯2−r2−λ2y2δ∗−δd21+d23,M3=r3−2r3k3x3δ∗−λ3y3δ∗−δd31+d32,M4=e1x1δ∗−m1,M5=e2x2δ∗−m2,M6=e3x3δ∗−m3.
And expanding the determinant we obtain,M6−λM5−λδ3d12d23d31+δ3d13d21d32+λδ2d12d21+λδ2d13d31+λδ2d23d32−λ3−δ2d23d32M1+λ2M1−δ2d13d31M2+λ2M2−λM1M2−δ2d12d21M3+λ2M3−λM1M3−λM2M3+M1M2M3M4−λ−e1δ2d23d32λ1x1δ∗−λ2λ1x1δ∗+λM2λ1x1δ∗+λM3λ1x1δ∗−λ1M2M3x1δ∗y1δ∗+e2λ2x2δ∗−δ2d13d31+λ2−λM1−λM3+M1M3M4−λ+e1−λλ1x1δ∗+λ1M3x1δ∗y1δ∗y2δ∗+e3λ3x3δ∗M5−λ−δ2d12d21+λ2−λM1−λM2+M1M2M4−λ+e1−λλ1x1δ∗+λ1M2x1δ∗y1δ∗+e2λ2λ2x2δ∗−λλ2M1x2δ∗−λλ2M4x2δ∗+λ2M1M4x2δ∗+e1λ1λ2x1δ∗x2δ∗y1δ∗y2δ∗y3δ∗=0
Simplifying we get,(66)λ6+N1λ5+N2λ4+N3λ3+N4λ2+N5λ+N6=0.where, N1,N2,N3,N4,N5,N6 are constant (see Appendix C).
Now, the Hurwitz matrix(67)L=N110000N3N2N1100N5N4N3N2N110N6N5N4N3N2000N6N5N400000N6,∇1=N1,∇2=N11N3N2,∇3=N110N3N2N1N5N4N3,∇4=N1100N3N2N11N5N4N3N20N6N5N4,∇5=N11000N3N2N110N5N4N3N2N10N6N5N4N3000N6N5.
The roots of the auxiliary equation have negative real parts if and only if all the principal diagonal minors of the Hurwitz matrix are positive, provided that∇1>0,∇2>0,∇3>0,∇4>0,∇5>0 and N6>0. Therefore, the equilibrium E2x1δ∗,x2δ∗,x3δ∗,y1δ∗,y2δ∗,y3δ∗ is stable iff ∇1>0,∇2>0,∇3>0,∇4>0,∇5>0 and N6>0.


## 4.1. Stability Condition of the Proposed System without Dispersal

At this pointE1x10∗,x20∗,x30∗,y10∗,y20∗,y30∗, the variation matrix is provided by,(53)VE1=−r1x10∗k1k¯12x10∗−k1−k¯100−λ1x10∗000−r2x20∗k2k¯22x20∗−k2−k¯200−λ2x20∗000−r3x30∗k300−λ3x30∗e1y10∗000000e2y20∗000000e3y30∗000,Theorem 3.
The system (1) is unstable without dispersal δ=0 under the condition ki+k¯i>2mi/ei.Proof.
The characteristic equation of the variation matrixVE1 is given by,(54)Δ=−r1k1k¯1x10∗2x10∗−k1−k¯1−λ00−λ1x10∗000−r2x20∗k2k¯22x20∗−k2−k¯2−λ00−λ2x20∗000−r3x30∗k3−λ00−λ3x30∗e1y10∗00−λ000e2y20∗00−λ000e3y30∗00−λ=0.
Evaluating, we get (see Appendix A),(55)cλ2+r3x30∗k3λ+e3λ3x30∗y30∗λ2+r2x20∗k2k¯22x20∗−k2−k¯2λ+e2λ2x20∗y20∗,λ2+r1x10∗k1k¯12x10∗−k1−k¯1λ+e1λ1x10∗y10∗=0.
Therefore, we have,(56)λ2+r3x30∗k3λ+e3λ3x30∗y30∗=0,(57)λ2+r2x20∗k2k¯22x20∗−k2−k¯2λ+e2λ2x20∗y20∗=0,(58)λ2+r1x10∗k1k¯12x10∗−k1−k¯1λ+e1λ1x10∗y10∗=0,
Now, we observe thatr3x30∗/k3>0 and e3λ3x30∗y30∗>0. Therefore, the roots of (56) are either all negative or have negative real parts.
Again2xi0∗−ki−k¯i<0⇔ki+k¯i>2mi/ei, i=1,2.
Then all the roots of (57) and (58) have positive real parts. Thus, the system (1) is unstable without dispersal δ=0 with the assumption ki+k¯i>2mi/ei, i=1,2.


## 4.2. Persistence of the Proposed System in the Presence of Dispersal

The biological meaning of persistence is the survival of all populations in future time. For the Kolgomorov type equations, persistence was discussed by Freedman and Waltman [64]. Mathematical definition for persistence of a system can be written as follows:Definition 1.
(Persistence): ifx0>0 and limt⟶∞infxt>0, xt persists. The system is said to persist if each component of the system persists.Theorem 4.
The system (1) persists.Proof.
Consider the matrix:(59)Af=D+diagf10,f20,f30,f40,f50,f60,where, D=dij is a matrix whose elements are given by,(60)dij=Dijfor i≠j,−∑k=1nDikαikfor i=j.
HereDij denotes the nonnegative diffusion coefficient of the species from the j-th patch to the i-th patch for i≠j and Dii=0, i=1,2,3,4,5,6. For the continuous diffusion case, the parameter αij≥0 corresponds to the boundary conditions. Now we have,(61)f1=r11−x1k1x1k¯1−1−λ1y1,f2=r21−x2k2x2k¯2−1−λ2y2,f3=r31−x3k3−λ3y3,f4=e1x1−m1,f5=e2x2−m2,f6=e3x3−m3.
Here we can recall the theorem which states that,sAf>0 implies liminft≥0xit≥δ>0 and liminft≥0yit≥δ>0 for i=1,2,3 where, sAf=maxReλ:λ∈the set of eigenvalues ofAf. Interested readers can follow reference [65] about the proof of this theorem. The characteristic equation for the matrix Af is,(62)detAf−λI=−r1−λ000000−r2−λ000000r3−λ000000−δd12+d13−m1−λδd12δd13000δd21−δd21+d23−m2−λδd23000δd31δd32−δd31+d32−m3−λ=0.
Which implies, (see Appendix B)(63)lλ+r1λ+r2λ−r3λ3+P+Q+Rλ2+PQ+QR+RP−δ2d12d21+d13d31+d23d32λ+PQR−δ2Pd23d32+Qd13d31+Rd12d21−δ3d12d23d31+d13d21d32=0.where, P=δd12+d13+m1, Q=δd21+d23+m2 and R=δd31+d32+m3. Clearly, one of the roots of (62) is positive. Therefore sAf>0. Consequently, the system (1) persists.


## 4.3. Stability in the Presence of Dispersal

Let us consider the understated theorem for the stability of the equilibriumE2x1δ∗,x2δ∗,x3δ∗,y1δ∗,y2δ∗,y3δ∗.Theorem 5.
The equilibriumE2x1δ∗,x2δ∗,x3δ∗,y1δ∗,y2δ∗,y3δ∗ is stable under some conditions.Proof.
The characteristic equation forVE2 can be written as,(64)r1k1k¯1x1δ∗2k1+2k¯1−3r1x1δ∗2k1k¯1−r1−λ1y1δ∗−δd12+d13−λδd21δd31−λ1x1δ∗00δd12r2k2k¯2x2δ∗2k2+2k¯2−3r2x2δ∗2k2k¯2−r2−λ2y2δ∗−δd21+d23−λδd320−λ2x2δ∗0δd13δd23r3−2r3k3x3δ∗−λ3y3δ∗−δd31+d32−λ00−λ3x3δ∗e1y1δ∗00e1x1δ∗−m1−λ000e2y2δ∗00e2x2δ∗−m2−λ000e3y3δ∗00e3x3δ∗−m3−λ=0.
Taking(65)M1=r1k1k¯1x1δ∗2k1+2k¯1−3r1x1δ∗2k1k¯1−r1−λ1y1δ∗−δd12+d13,M2=r2k2k¯2x2δ∗2k2+2k¯2−3r2x2δ∗2k2k¯2−r2−λ2y2δ∗−δd21+d23,M3=r3−2r3k3x3δ∗−λ3y3δ∗−δd31+d32,M4=e1x1δ∗−m1,M5=e2x2δ∗−m2,M6=e3x3δ∗−m3.
And expanding the determinant we obtain,M6−λM5−λδ3d12d23d31+δ3d13d21d32+λδ2d12d21+λδ2d13d31+λδ2d23d32−λ3−δ2d23d32M1+λ2M1−δ2d13d31M2+λ2M2−λM1M2−δ2d12d21M3+λ2M3−λM1M3−λM2M3+M1M2M3M4−λ−e1δ2d23d32λ1x1δ∗−λ2λ1x1δ∗+λM2λ1x1δ∗+λM3λ1x1δ∗−λ1M2M3x1δ∗y1δ∗+e2λ2x2δ∗−δ2d13d31+λ2−λM1−λM3+M1M3M4−λ+e1−λλ1x1δ∗+λ1M3x1δ∗y1δ∗y2δ∗+e3λ3x3δ∗M5−λ−δ2d12d21+λ2−λM1−λM2+M1M2M4−λ+e1−λλ1x1δ∗+λ1M2x1δ∗y1δ∗+e2λ2λ2x2δ∗−λλ2M1x2δ∗−λλ2M4x2δ∗+λ2M1M4x2δ∗+e1λ1λ2x1δ∗x2δ∗y1δ∗y2δ∗y3δ∗=0
Simplifying we get,(66)λ6+N1λ5+N2λ4+N3λ3+N4λ2+N5λ+N6=0.where, N1,N2,N3,N4,N5,N6 are constant (see Appendix C).
Now, the Hurwitz matrix(67)L=N110000N3N2N1100N5N4N3N2N110N6N5N4N3N2000N6N5N400000N6,∇1=N1,∇2=N11N3N2,∇3=N110N3N2N1N5N4N3,∇4=N1100N3N2N11N5N4N3N20N6N5N4,∇5=N11000N3N2N110N5N4N3N2N10N6N5N4N3000N6N5.
The roots of the auxiliary equation have negative real parts if and only if all the principal diagonal minors of the Hurwitz matrix are positive, provided that∇1>0,∇2>0,∇3>0,∇4>0,∇5>0 and N6>0. Therefore, the equilibrium E2x1δ∗,x2δ∗,x3δ∗,y1δ∗,y2δ∗,y3δ∗ is stable iff ∇1>0,∇2>0,∇3>0,∇4>0,∇5>0 and N6>0.


## 5. Hopf Bifurcation

When the parameter of a dynamical system changes due to a sudden qualitative change in its behavior, then a bifurcation of the system occurs. Hopf bifurcation occurs in the case of nonhyperbolic nonlinear equations for two or more dimensions. It is typically happened in a differential equation for switching the eigenvalues to become purely imaginary. The fixed point switches from a stable focus to an unstable one when the real part of the eigenvalue changes from negative to positive; such a bifurcation is known as supercritical. On the contrary, the bifurcation is subcritical when the fixed point changes from an unstable focus to a stable one for switching the real part of the eigenvalue from positive to negative.We observed a special type of Hopf bifurcation for the system proposed in this study. We find a pair of complex conjugate eigenvalues that pass through the imaginary axis of the Jacobian matrix, and all other eigenvalues consist of negative real parts. Bifurcated limit cycles can be observed in a supercritical Hopf bifurcation both physically and numerically. To distinguish this type of bifurcation with the same for the eigenvalues on the right half-plane, let us call them a simple Hopf bifurcation. In this study, we are interested in exploring the possible occurrence of a simple Hopf bifurcation aroundE2x1δ∗,x2δ∗,x3δ∗,y1δ∗,y2δ∗,y3δ∗ with δ as a bifurcating parameter.The following theorem provides the necessary and sufficient condition for the occurrence of a simple Hopf bifurcation of the system (1):Theorem 6.
The necessary and sufficient conditions for the system (1) to undergo Hopf bifurcation at δ=δ0 (point of bifurcation) around the interior equilibrium point E2x1δ∗,x2δ∗,x3δ∗,y1δ∗,y2δ∗,y3δ∗ are that N6δ0>0, ∇1δ0>0, ∇2δ0>0, ∇3δ0>0, ∇4δ0>0, ∇5δ0=0 and d∇5δ0/dδ≠0.Proof.
The characteristic equation forVE2 is(68)λ6+N1δλ5+N2δλ4+N3δλ3+N4δλ2+N5δλ+N6δ=0.where, Niδi=1,2,3,4,5,6 is a smooth function of δ.
The Hurwitz matrix(69)Lδ=N1δ10000N3δN2δN1δ100N5δN4δN3δN2δN1δ10N6δN5δN4δN3δN2δ000N6δN5δN4δ00000N6δ,∇1δ=N1δ,∇2δ=N1δ1N3δN2δ,∇3δ=N1δ10N3δN2δN1δN5δN4δN3δ,∇4δ=N1δ100N3δN2δN1δ1N5δN4δN3δN2δ0N6δN5δN4δ,∇5δ=N1δ1000N3δN2δN1δ10N5δN4δN3δN2δN1δ0N6δN5δN4δN3δ000N6δN5δ.
Using the condition for simple Hopf bifurcation [66] at δ=δ0 around E2x1δ∗,x2δ∗,x3δ∗,y1δ∗,y2δ∗,y3δ∗, we get N6δ0>0,∇1δ0>0,∇2δ0>0,∇3δ0>0,∇4δ0>0,∇5δ0=0 and d∇5δ0/dδ≠0. Hence the theorem is proved.


## 6. Numerical Simulation

For the numerical verification of analytical findings of the previous sections, computer simulations have been performed to obtain a different graphical presentation of the proposed model. We used Matlab and Wolfram Mathematica for this purpose, and the figures obtained are extremely imperative from a realistic viewpoint. We have used a set of hypothetical data to examine the effect of dispersal speed on the patches. So, we decided to do the simulations in two different cases. Case11 attributes to the prey-predator system when the prey population budges independently in its patches, i.e., dispersal does not happen δ=0. For Case 12, we consider the shift in the prey population between patches δ≠0. To simulate the proposed model, we take the value of the parameters as provided in Table 1.Table 1
Numerical values of the model parameters.
ParameterValueParameterValuer12.6e20.381k14.29e30.825k1¯0.43m10.75λ10.69m20.8r20.5m31.6k24.19d120.419k2¯0.47d210.429λ20.6d130.36r31.43d310.429k33.6d230.36λ30.71d320.419e10.32δ7.7Case 11.
Simulations without dispersalδ=0.
Based on the model parameters as prearranged in Table1, we achieve the inner equilibrium point E12.34375,2.09974,1.93939,7.60821,1.44154,0.929056 for the system (1) without dispersal. Now we observe that, 4.72=k1+k1¯>2m1/e1=4.6875 and 4.66=k2+k2¯>2m2/e2=4.1995. Therefore ki+ki¯>2mi/ei for i=1,2 and the inner equilibrium point E12.34375,2.09974,1.93939,7.60821,1.44154,0.929056 is unstable in the absence of dispersal (Theorem 3). Time series graph of the species with initial condition x10,x20,x30,y10,y20,y30=3.1,3.5,3.7,0.3,0.4,0.5 depicted in Figure 2.
Clearly, the system becomes unstable with time progress in the absence of dispersal (Figure2). Without dispersal speed δ=0, trophic cascades of the prey and predator species with respect to ki in all patches are presented through Figures 3 to 6.
From Figures3 to 5, we observe that with the increasing values of ki, the i th i=1,2,3 predator yi gradually increases and remaining species remain constant throughout. Following the above cases, we found that, the changing carrying capacity only affects a certain patch. Now in Figure 6, we consider that, all the environmental carrying capacities are equal, that is, k1=k2=k3=k (say) and observe that, all predator species gradually increase, while all prey species remain constant throughout as k increases.
Now, the impact of the mortality rate of the predator speciesmii=1,2,3 on the equilibrium components of the prey-predator populations is represented through Figures 7 to 10 when there is no dispersal.
From Figures7 and 8 we observed that, as mi increases in the i th i=1,2 patch, then the xii=1,2 species gradually increases and the yii=1,2 species follow the parabolic path and all other species remain constant throughout. It is clear from Figure 9 that, when m3 increases, x3 species gradually increases and the corresponding predator species y3 gradually decreases and all other species of the system remain constant throughout. In Figure 10 we consider the death rate of all predator species as equal i.e., m1=m2=m3=m (say) and observe that the three prey species xi gradually increase, predator species yii=1,2 follow the parabolic path and the remaining predator species y3 gradually decrease as m progresses.
We have studied the influence of environmental carrying capacity and death rate of predators on species without dispersal. We found that parameters such ask1, k2, k3, m1, m2 and m3 play a vital role in the population dynamics of the system (1).Figure 2
Unstable behavior of the proposed prey-predator system as time progress.Figure 3
Behaviour of all species with respect tok1.Figure 4
Behaviour of all species with respect tok2.Figure 5
Behaviour of all species with respect tok3.Figure 6
Behaviour of all species with respect tok.Figure 7
Effect ofm1 on all the species.Figure 8
Effect ofm2 on all the species.Figure 9
Effect ofm3 on all the species.Figure 10
Effect ofm on all the species.Case 12.
Simulations with dispersalδ≠0.
Here we assume the positive dispersal of the system (1) and taking δ=7.7 and other parameters of the system are the same as given in Table 1. This dispersal speed gives the interior equilibrium point E22.34375,2.09974,1.93939,7.16543,2.28463,0.677693. Now, the main objective is to study the effect of dispersal on the equilibrium point E2 as well as its stable nature.
The time series graph of the system (1) forδ=7.7 with the initial condition x10,x20,x30,y10,y20,y30=3.1,3.5,3.7,0.3,0.4,0.5 is presented in Figure 11.
Figure11 reflects the stable behavior of the proposed prey-predator system in the presence of dispersal. Here we found that each of the species of system (1) gradually converges to the inner equilibrium point E2 as time progresses. This observation concludes that the proposed system (1) persists in the presence of dispersal (Theorem 4). Therefore, the stability of the proposed model is completely dependent on the dispersal speed δ.
The effect ofk1,k2, k3, and kk1=k2=k3=k on the system (1) is presented in Figures 12 to 15.
From Figure12, we observe that the predator species of the first patch gradually increases and all other species remain constant with the progression of k1. In Figures 13 and 14, the predator species of the ith patch gradually increase as kii=2,3 progress. In Figure 15, we consider k1=k2=k3=k (say) and observe that, all predator species gradually increase as k progresses and all prey species maintain a constant level throughout. Therefore, we conclude that carrying capacity can affect only the predator levels with the dispersal of prey species.
Again, to observe the effect of predator mortality rate on species, we consider Figures16 to 19.
From Figure16, we observe that, the predator species in patch one gradually decreases. Consequently, prey species in patch one gradually increases. Also, predator species in patches two and three increase progressively, and the remaining two prey species maintain constant levels. The same types of observation are seen at in Figures 17 and 18, respectively. In Figure 19, we consider m1=m2=m3=m (say) and observe that the predator species y1 follow a parabolic type path, the species y2 and y3 increases first, then decreases and the three prey species gradually increase as m increases. Therefore, we conclude that the death rate of predators affects all the population levels.
Finally, we observe the effect of dispersal speedδ on the equilibrium components of the system (1). For different values of dispersal, the equilibrium points are provided in Table 2.
Table2 reflects that the population level of the prey species xii=1,2,3 at the equilibrium level is independent of δ as we obtained xi=mi/ei for i=1,2,3 in the theoretical section. Therefore, equilibrium levels of the prey population are not affected by the dispersal speed in any patch. On the contrary, the population level of predator species is changing gradually.
The pictorial presentation of the trophic cascade of the system (1) with increasing dispersal speed is described in Figure 20.
Here we found that, the predator population is very sensitive to changing values ofδ and therefore the dispersal speed is important for the persistence of the three patch-based systems.
Now, if we gradually increase the value ofδ, then the solutions of the proposed model system near E2 shift from its unstable periodic limit cycle oscillations to the stable steady-state through Hopf bifurcation. This dynamical property via Hopf bifurcation has been shown in Figure 21. Here, we have taken all the parameter values as previously and δ as the bifurcation parameter.
All figures show a subcritical Hopf bifurcation localized atδ=0.35. The proposed system undergoes a limit cycle for δ<0.35. The equilibrium point of co-existence is stable locally and asymptotically for δ>0.35 as shown in Figure 21.Figure 11
Time series plot of the system (1) with initial condition 3.1,3.5,3.7,0.3,0.4,0.5.Figure 12
Variation of the system (1) with respect to k1.Figure 13
Variation of the system (1) with respect to k2.Figure 14
Variation of the system (1) with respect to k3.Figure 15
Variation of the system (1) when k=k1=k2=k3.Figure 16
Effect ofm1 on the system (1).Figure 17
Effect ofm2 on the system (1).Figure 18
Effect ofm3 on the system (1).Figure 19
Effect ofm=m1=m2=m3 on the system (1).Table 2
Equilibrium points for differentδ.
δEquilibrium point (E2)02.34375,2.09974,1.93939,7.60821,1.44154,0.92905622.34375,2.09974,1.93939,7.4932,1.66052,0.86376742.34375,2.09974,1.93939,7.37819,1.87951,0.79847862.34375,2.09974,1.93939,7.26319,2.0985,0.73318982.34375,2.09974,1.93939,7.14818,2.31748,0.667899102.34375,2.09974,1.93939,7.03317,2.53647,0.60261122.34375,2.09974,1.93939,6.91817,2.75546,0.537321142.34375,2.09974,1.93939,6.80316,2.97444,0.472032162.34375,2.09974,1.93939,6.68815,3.19343,0.406743182.34375,2.09974,1.93939,6.57315,3.41242,0.341454202.34375,2.09974,1.93939,6.45814,3.6314,0.276165Figure 20
Trophic cascade of all species in our proposed system with changing dispersal rate.Figure 21
Bifurcation diagrams with respect toδ.


## 7. Discussion and Conclusion

In this study, a three-patch based prey-predator system is presented with a strong Allee effect in the first two patches. This is the first attempt to study the dynamics of the prey-predator system in three-patch environments with the Allee effect in the first two patches. The impact of the dispersal speed on the stability and persistence of the proposed system is elaborately studied in this paper. For a better understanding of the proposed model, we provide the concept of balanced dispersal. We observed that the equilibrium level of the prey species depends on the mortality rate of the predator species as well as the conversion rate of the prey biomass to the predator biomass. These population levels are not affected by the speed of dispersal. But the speed of dispersal has a definite effect on the equilibrium level of predator species. In addition, the stable nature of the equilibrium point depends on the mortality rate of the predator in each patch. The stability criterion of the system in the absence of dispersal (Theorem3) can help to maintain the ecological balance. But the system becomes stable in the presence of positive dispersal speed (Theorem 5). Therefore, we conclude that the mortality rate of the predator species and the dispersal speed play a vital role in the dynamics of the supposed ecological system.The numerical verification of the analytical findings is vital to the biological relevance of this type of model. Moreover, graphical views always support the findings and contribute to the vast future development of the ecosystem. Therefore, the proposed system may contribute to the dynamics of the prey-predator in the assumed situations.Finally, we conclude that the proposed system in three-patch environments is very interesting and shows complex dynamics. But the model can also be modified by taking the Allee effect in all three patches. Predator species can also be assumed to move between patches, making the model more interesting and realistic.



---

*Source: 7919952-2022-06-20.xml*