# Pattern Dynamics in a Predator-Prey Model with Diffusion Network

Authors: Wenjie Yang; Qianqian Zheng; Jianwei Shen; Qing Hu
Journal: Complexity
 (2022)

---

## Abstract

Diffusion plays an essential role in the distribution of predator and prey. We mainly research the diffusion network’s effect on the predator-prey model through bifurcation. First, it is found that the link probability and diffusion parameter can cause Turing instability in the network-organized predator-prey model. Then, the Turing stability region is obtained according to the sufficient condition of Turing instability and the eigenvalues’ distribution. Finally, the biological mechanism is explained through our theoretical results, which are also illustrated by numerical simulation.


---

## Body

## 1. Introduction

Turing instability was first investigated in the reaction-diffusion system [1], and it is constantly being promoted to explain the dynamical mechanism [2, 3]. Asllani et al. pointed out that the directed network could induce Turing instability when an indirect network does not work [4, 5]. Meanwhile, tuning the topology structure of the system can create or destroy patterns in a reaction-diffusion system [6]. Mimar et al. proved that the pattern formation’s topological properties are determined by complex interaction [7]. Although spontaneous patterns [8] are associated with the dominance of eigenvectors and eigenvalues [9–11], the dynamical mechanism of the random network in pattern formation remains to be uncovered.Since the Lotka–Volterra model was proposed in the early twentieth century, some biological mechanisms were explained in the predator-prey system [12–18]. Chang et al. researched the dynamics in the predator-prey system on complex networks and found that Turing instability caused by delay can generate spiral waves [19]. Liu et al. showed the effect of network and diffusion on the ecological balance of the predator-prey system [20]. Upadhyay and Bhattacharya studied the differences between the aqueous and terrestrial environments in predator-prey networks and tried to explain their biological mechanism [21]. Astarloa et al. tried to use the joint species distribution modeling to reveal the coexistence problem of prey and predator in the Bay of Biscay [22]. Although previous work shows that diffusion and randomness influence the Turing pattern significantly, the random network’s effect on the distribution of predator and prey should be carefully evaluated.To reveal the natural mechanism of biological invasion, we intend to investigate the random diffusion network’s effect on the network-organized predator-prey model’s stability. First, we obtain the conditions under which Turing bifurcation arises. Second, we find an estimated region of all the eigenvalues of the Laplacian matrix, the sufficient stability conditions in the network-organized predator-prey model. Third, we explain the network-organized Turing instability by the mean-field approximation and comparison principle. Also, we estimate the Turing instability range about link probability and diffusion and try to explain the mechanism of biological invasion. Last, we illustrate our theoretical results through numerical simulation.


## 2. A Network-Organized Predator-Prey Model

For the convenience of subsequent research, we first give some necessary symbolic rules. The network Laplacian matrixL=Lij can be treated as Lij=Aij−kiδij, and all the eigenvalues of L are Λ=Λi. The eigenvectors ϕi∈Rn of Λi satisfy Lϕi=Λiϕi,i=1,…,n. L is real symmetric, and we choose an orthonormal basis for ϕiϕj=δij, where the degree of node i is ki, and δij is the Kronecker delta function.We consider the following predator-prey system:(1)dxdt=xr1−a11x−a12y,dydt=yr2+a21x−a22y,where x and y are prey and predator, respectively. r1 represents the intrinsic growth rate, a12 is the proportionality coefficient of predator and prey, and a21 is the growth rate of the predator; the predator increases exponentially with ratio r2 (actually r2<0 means the number of predators decreases exponentially). According to Samuelson’s assumptions [23], a11 and a22 represent increasing returns or decreasing returns, respectively. Among them, aii>0<0,i=1,2 correspond to increasing returns (decreasing returns). System (1) correspond to the mixed-income when a11a22≤0.System (1) always has three equilibria E0≜0,0,E1≜r1/a11,0, and E2≜0,r2/a22. System (1) also have fourth equilibrium point E3≜x∗,y∗≜r1a22−r2a12/a11a22+a12a21,r1a21+r2a11/a11a22+a12a21 when r1a22−r2a12/a11a22+a12a21>0,r1a21+r2a11/a11a22+a12a21>0. On the basis of Hopf bifurcation’s definition, we need to satisfy trJE3=0,detJE3>0. The critical value for bifurcation is a positive root of trJE3=0 and the bifurcation parameter is a12=a12∗ which satisfies detJE3>0,(2)a12=r1+r2a11a22+r1a21a22r2a11≜a12∗,r2r1a21+r2a11a11+a21<0.Then, we give the condition of Hopf bifurcation. Firstly, we verify the transversality condition [24].(3)dda12trJx∗,y∗|a12=a12∗=−a22a11+a21r1a21+r2a11a11a22+a12a212|a12=a12∗=−r2a112a22a11+a21r1a21+r2a11≠0.Then, we compute the first Lyapunov coefficient [25]. Let a12=a12∗, then x∗c,y∗c=−r2/a11+a21,r2a11/a22a11+a21. Making the following shift(4)X=x−x∗c,Y=y−y∗c.(1) becomes(5)X˙=α10X+α01Y+α20X2+α11XY,Y˙=β10X+β01Y+β11XY+β02Y2,where α10=−a11x∗c,α01=−a12x∗c,α20=−a11,α11=−a12,β10=a21y∗c,β01=−a22y∗c,β11=a21,β02=−a22.The first Lyapunov coefficient ata12=a12∗ can be computed by the formula(6)l1=−3π2α01Δ3/2α10α01β112−2α202+α20β11+α11β02+α10β10α112−2β022+α11β02+α01β10β11β02−α11α20−2α102β11β02+2α102α11α20=−3π2−a11+a21r2r1a21+r2a113/2r1a21+r1a11+r2a11a222a11+a21,where Δ=α10β01−α01β10>0,α10+β01=0, and if l1<0>0, the Hopf bifurcation is supercritical (resp. subcritical).Whenr1a22−r2a12/a11a22+a12a21>0,r1a21+r2a11/a11a22+a12a21>0, E3 changes its stability, and Hopf bifurcation occurs if a12=a12∗.Typically, the distribution of individuals is spatially heterogeneous. So, we research (1) with a reaction-diffusion and network as follows:(7)dxidt=xir1−a11xi−a12yi+d1∇2xi,dyidt=yir2+a21xi−a22yi+d2∇2yi,where d1,d2 are the diffusion constants.Generally, we can regard diffusion as an interplay between network nodes. In this article, we consider a distinctive interaction between nodes. So the network-organized system (7) is(8)dxidt=xir1−a11xi−a12yi+d1∑jLijxi,dyidt=yir2+a21xi−a22yi+d2∑jLijyi.The linearized network-organized system of the system (8) is(9)dxidt=α10xi+α01yi+d1∑jLijxi,dyidt=β10xi+β01yi+d2∑jLijyi.The general solution of the linear network-organized system can be expressed as [3].(10)xi=∑k=1Nckβkeλktϕik,yi=∑k=1Nckeλktϕik,where ∑jLijϕjk=Λiϕik.Substituting the general solution into system (9), the Jacobian matrix Bii=1,…,n is(11)Bi=α10+d1Λiα01β10β01+d2Λi,where all the eigenvalues of matrix L can be represented as Λi0=Λ1>Λ2>⋯>ΛN. Then, system (8) has the following characteristic function:(12)λ2+pΛiλ+qΛi=0,where pΛi=−α10+β01−d1+d2Λi,qΛi=d1d2Λi2+α10d2+β01d1Λi+α10β01+α01β10. Turing instability occurs when there is a Λi with Reλ>0. From (12),(13)4d1d2α10β01+α01β10+α10d2+β01d12=0.We can get the Turing instability’s critical value aboutΛ (Skim d1=0.1562d2) in the reaction-diffusion system Figure 1. Note k1c2,k2c2k1c2<k2c2 as two solutions of Reλk2=0, Reλk2>0 holds if k2∈B=λ|k1c2<λ<k2c2 Figure 2d2=2. Turing instability of the reaction-diffusion system is the Turing instability’s prerequisite in the network-organized system. Based on the Gershgorin circle theorem [9, 10], we have.(i)
Result 1:kmax=maxki,kmin=minki and Λi is the eigenvalue of the Laplacian matrix L, then, Λi∈C=Λi|−2kmax<Λi≤0.(ii)
Result 2: in a network-organized system, a system remains stable when no eigenvalue ofL stays at the instability range Λ∩B=Φ, Turing instability occurs when Λ∩B≠Φ, and B∩C≠Φ (C is the set of the eigenvalues of Laplacian matrix) induces the occurring of instability.Figure 1
The stability of the system (1) without diffusion when r1=1,r2=−1,a11=−1.158,a12=1,a21=2,a22=1.Figure 2
The stability of the system (1) without diffusion when r1=1,r2=−1,a11=−1.158,a12=1.7271,a21=2,a22=1.To study the mechanism of the network-organized system’s stability, we research (8) through the comparison principle: let xi=x∗+εx^i and yi=y∗+εy^i, where 0<ε≪1. Substituting xi,yi in (9), expanding in ε, we can get the linearized system:(14)dx^idt=α10x^i+α01y^i+d1Lx^i,dy^idt=β10x^i+β01y^i+d2Ly^i.We resolve the first-order perturbations intoϕi’s eigenfunction expansions, to consider the system’s stability. Let(15)x^i=X^iϕi,Y^i=Y^iϕi,for eachi=1,…,n.Substituting (15) into (20) and noting the properties of ϕi, we obtain(16)dX^idt=α10X^i+α01Y^i+d1ΛiX^i,dY^idt=β10X^i+β01Y^i+d2ΛiY^i.Using the comparison principle from [26], we note(17)d2Ydt2+PtdYdt+QtY=0,where Qt<0. We eliminate one of X^i or Y^i from (16), then(18)d2X^idt2+pΛidX^idt+qΛiX^i=0,d2Y^idt2+pΛidY^idt+qΛiX^i=0,where pΛi,qΛi are as in (12). Then, we obtain the generalized condition:(19)α10β01+α01β10+d1β01+d2α10Λi+d1d2Λi2<0.Assume that the instability condition (19) holds. Then, the homogeneous state xk,yk=x∗,y∗ is unstable under the ith Turing mode (15).We also consider the system (8) through the mean-field theory:(20)dxidt=fxi,yi+d1Hx−kixi,dyidt=gxi,yi+d2Hy−kiyi,where Hx=∑j=1nAijxj,Hy=∑j=1nAijyj and ki is the network’s degree. We let other nodes stay at equilibrium x∗,y∗, and rewrite the single-node system,(21)dxidt=xir1−a11xi−a12yi+d1x∗−kixi,dyidt=yir2+a21xi−a22yi+d2y∗−kiyi.Therefore, it is easy to obtain the characteristic equation:(22)λ2+pkiλ+qki=0,where pki=−α10+β01+d1ki+d2ki,qki=α10β01+α10d2+β01d1ki+d1d2ki2. Assume λ1 and λ2 are two eigenvalues of system (22). Then system (21) is stable when d1=d2=0, namely pki>0, so λ1+λ2=−pki<0. If qki>0qki<0. Then, the system (21) is stable(unstable). To sum up, we can use the comparison principle(the mean-field theory) to conclude that the system (8) is unstable when qΛi<0qki<0.


## 3. Numerical Simulation

We give some numerical analysis based on the earlier theoretical analysis. From Figure3, if we choose r1=1,r2=−1,a11=−1.158,a12=1,a21=2,a22=1, the equilibrium point x∗,y∗ is asymptotically stable, namely, ecological balance persists as long as there is no biological invasion. That means the predator and the prey can coexist. The equilibrium x∗c,y∗c=1.1326,1.2651 is unstable, and Hopf bifurcation occurs Figure 4, when a12 passes through the critical value a12∗=1.7271. Because l1=8.2123>0, the Hopf bifurcation is subcritical, the prey-predator system shows periodic changes; thus, this state is easy to destroy. The equilibrium point x∗,y∗ becomes unstable, when a12=1.8271>a12∗.Figure 3
Turing bifurcation aboutd2 and d1 in the continuous system.Figure 4
Linear stability analysis. The relationship betweenReλ about −Λik2 when r1=1,r2=−1,a11=−1.158,a12=1,a21=2,a22=1, and p=0.1,d2=2. The critical values are Λc=−3.35,d1c=0.3124.As diffusion is a vital factor in the distribution of predator and prey, we should not ignore migration. So we construct the random network and transform it into the Laplacian matrixLij. And we consider how the random network affects the equilibrium point’s stability when (7) is stable.Based on the above theoretical analysis, the Turing instability in the reaction-diffusion equation is a precondition for Turing instability in a network-organized system. Turing instability occurs in a reaction-diffusion equation when Turing bifurcation parametersd2 and d1 are in a certain range Figure 1. Namely, diffusion behavior is universal and allowed, but the relative diffusion rate needs to be within a specific range; otherwise, the ecosystem will be out of balance and destroyed.To observe Turing instability in the network-organized, we should guaranteed1<0.1562d2 while changing p. However, the Laplacian matrix eigenvalues are Λ1,Λ2,…,ΛN,N is the number of nodes, and k2 is continuous in a reaction-diffusion system. In other words, the distribution of −Λi determines the system stability. The critical point is d1c=0.3124 when d2=2. The range of eigenvalues Λi could be obtained by result 1. The eigenvalues of L are discrete and included in the continuous region Figure 2.The red region shows the estimated range of Turing instability aboutp,d1 in Figure 5. Besides, we obtain the estimated range of p, d2lnN/d1cN>p>1/N2, which determines the system dynamical behavior and network characteristics. Moreover, p=lnN/N is the critical value between the connected network and sparse network [27]. Anyway, the invasive rate of species and p in the predator-prey network play a vital role in balancing the native biological system.Figure 5
The region of instability aboutd1 and p.Let the value ofp,d1 lie in the blue region shown in Figure 5 (here, we set p=0.00006,d1=0.1124), then, there is no point Λi (k2) in the instability region Λ∩1.7,21.1=Φ (result 2) (Figure 6(a)) and the equilibrium is stable (Figure 6(a)). We verified the correctness of the previous analysis with bifurcation diagrams (Figure 6(b)). We can conclude that with a fixed predator invasion rate and a fixed prey invasion rate, the link probability between two different systems will negatively impact the entire predator-prey network stability. Even the predator-prey network’s diffusion induces the species extinction.Figure 6
(a) Instability region (left) of random network and pattern formation (right) in random network whenp=0.00006,d1=0.1124. (b) The bifurcation about d1 when p=0.00006.
(a)(b)Figure7 shows an example of p that belongs to the red region in Figure 5, such as p=0.08,d1=0.1124<d1c=0.3124, almost all of Λik2 in instability region Figure 7(a) left(up) and 7(a) right(up) shows that Turing instability occurs. If p=0.08 and d1=1.1124>d1c=0.3124 (i.e., p,d1 in the blue region in Figure 5), there is no point Λik2 in the instability region (Figure 7(a), left(down)), the system is stable (Figure 7(a), right(down)). We also verified the correctness of the previous analysis with a bifurcation diagram (Figure 7(b)).Figure 7
(a) The relationship betweenReλ about −Λik2 (left) and the pattern formation of system (3) (right). Top: when p=0.08,d1=0.1124, almost all points Λik2 lie in the instability region and the Turing instability occurs. Bottom: when p=0.08,d1=1.1124, no point Λik2 lies in the instability region and the system is stable. (b) The bifurcation occurs about d1 when p=0.08.
(a)(b)The ecological interpretation is that the system is less stable when the invasion rate of the prey is lower when the link probability of predator-prey systems is fixed. Conversely, if we want to keep the predator-prey system stable, the invasion rate of the prey needs to be decreased.


## 4. Conclusion

This paper theoretically derives Turing instability conditions [28–31] in a predator-prey network and carries out a detailed numerical study. We study the effects of diffusion and link probability on pattern formation in a random system. The smaller the invasion rate d1 of the prey undergoes, the easier the predator-prey network Turing instability. Meanwhile, Turing instability occurs when the link probability p falls in the approximate region d2lnN/d1cN>p>N2/1, which depends on the system characteristics.Regarding an explanatory biological mechanism, the link probability and the invasion rate significantly influence the entire predator-prey network stability. The lower the invasion rate of the prey is, the less stable it is. Conversely, the prey invasion rate can be increased if we want to keep the predator-prey system stable. That means we can find an equilibrium point between native and invasive species by adjusting the diffusion probability and species’ invasion rate. Finally, we obtain the estimated region ofp, d2lnN/d1cN>p>1/N2. Thus, invasion rate of species and the link probability in the predator-prey network play a key role in balancing the native biological system.



---

*Source: 9055480-2022-07-31.xml*