Delay-induced self-organization dynamics in a prey-predator network with diffusion

Considering that time delay (delay) is a common phenomenon in biological systems, reaction-diffusion equations with delay are widely used to study the dynamic mechanism of those systems, in which delay can induce the loss of stability and degradation of performance. In this paper, taking into account the inhomogeneous distribution of species in space and this can be considered as a random network, the pattern dynamics of a prey-predator network system with diffusion and delay are investigated. The effect of delay and diffusion on the network system is obtained by linear stability analysis, including the stability and Hopf bifurcation as well as Turing pattern. Our results show that the stability of the system changes with the value of delay. Moreover, we obtain Turing pattern related to the network connection probability and diffusion. Finally, the numerical simulation verifies our results.


Introduction
Colorful patterns emerge from reaction-diffusion equations and can explain a wide variety of natural phenom-Q. Hu · J. Shen [1], the Turing pattern has received increasing attention in biology, physics, chemicals, and other fields [2][3][4][5]. At present, it has been applied to more complex interdisciplinary subjects. The dynamics mechanics of the prey-predator system has always been a hot topic [6][7][8]. For the prey-predator model with Holling type II, Aziz-Alaoui and Daher [9] proposed a modified predator-prey model based on [10], and they studied the boundedness and global stability of the model.
Since the species in prey-predator systems are subjected to not only time evolution but also inhomogeneous space evolution, and Turing [1] first used the reaction-diffusion system to successfully discovered that the introduction of diffusion could cause the stable state to become unstable under certain conditions and obtained Turing pattern phenomena, so many scholars began to study the pattern formations of prey-predator models by using the reaction-diffusion equation [11][12][13][14][15], where the diffusion mainly refer to the movement of species within a community. Camara and Aziz-Alaoui [16] investigated a modified Holling-Tanner model and found Turing and Hopf bifurcations, they also derived the bifurcation conditions. Then, for a diffusive predator-prey model with the Holling-Tanner functional response, Shi et al. [17] investigated the stability and the existence of a positive solution.
However, the inhomogeneous diffusion (such as the global spread of the virus, the interaction and migration of populations, and the global economic crisis) exists not only in one region but in multiple regions, and these regions can be abstracted as a complex network in the real world, so we can study the pattern dynamics of a prey-predator diffusion system in network. For the network diffusion system, Othmer and Scriven [18] proposed a general method for analyzing the instability of reaction-diffusion model in small networks, and they obtained the effect of network topology changes on the stability and pattern. Then, Aly and Farkas [19] studied the Turing bifurcation of a diffusion prey-predator system in two patches, where the migration of populations among patches causes diffusion. For the multi-patch systems, Jansen and Lloyd [20] used the eigenvalues of matrix to study the local stability of spatially homogeneous solutions. Subsequently, Turing patterns in complex organized networks were studied in [21][22][23]. For the FitzHugh-Nagumo model in a network-organized system, Zheng and Shen [24] recently studied Turing instability, where the network structure and diffusion lead to Turing instability. Liu et al. [25] placed the preypredator model into complex networks and found that the network topology and diffusion rate can affect the Turing pattern formation.
It is well known that delay is a universal phenomenon in nature, such as the food we digest, the gestation or maturation period of population reproduction in the biological system and the latent period of some infectious disease. And with the further study, researchers found that the prey-predator system with delay can be closer to the reality [26][27][28][29][30]. For the delay system, [31] and [32] provided two methods to investigate the influence of delay on the stability and bifurcations. By constructing a Lyapunov function for a delay predator-prey model, Nindjin et al. [33] investigated the global stability of the positive equilibrium, and obtained the sufficient condition. For a Lotka-Volterra competition system, Song et al. [34] studied the stability and the local Hopf bifurcation as well as the Hopf bifurcation direction of the positive equilibrium. Then, Yang and Zhang [35] considered a diffusive preypredator system with delay, the conditions of Turing instability and Hopf bifurcation are obtained by using linear analysis.
If we study the prey-predator system with diffusion and delay in network, the dynamic results obtained will reflect real population development trend more accurately. However, although many researchers investigated the network diffusion systems, there are few studies on the network systems with delay. Recently, Chang et al. [36] investigated a predator-prey diffusive model in complex networks and found that delay can induce pattern dynamics. Wang and Tian [37,38] found that the network structure hardly has an impact on the stability of this model with small delay, and a large time delay can induce the thick-tailed phenomenon of evolution patterns. Therefore, considering the random migration and inhomogeneous spatial distribution of prey and predators among different communities, we mainly investigate the effect of delay and diffusion network on dynamical behaviors of a prey-predator diffusion network with Holling type II and modified Leslie-Gower functional response.
The rest of this paper is organized as follows. For the prey-predator diffusion network system with delay, we derive the stability condition of the positive equilibrium without delay by linear stability analysis in Sect. 2. Then, in Sect. 3, we use two methods to study the effect of delay on the stability, and obtain the conditions of Hopf and Turing bifurcations caused by delay and network diffusion. In Sect. 4, consistent results are obtained by numerical simulation. Finally, the above results are summarized in Sect. 5.

Stability analysis of the prey-predator diffusion network without delay
The modified Leslie-Gower prey-predator model with a Holling type II functional response can be written as follows: where U and V represent the density of prey and predator at time T , respectively. K , r 1 , r 2 , D, A, γ , W are model parameters with positive values. K denotes the carrying capacity and r 1 represents the intrinsic growth rate of prey. DU U +A is known as the Holling type II functional response [6], D is the maximum capture rate per predator, A is the half-saturation constant which satisfies DU U +A | U =A = D 2 . r 2 is the intrinsic growth rate of predator, γ V U +W is known as the modified Leslie-Gower term [9], γ is the number of prey required per predator at equilibrium point when V is equal to U +W γ , and W measures the degree of which the environment protects predators.
For simplicity, make the following dimensionless transformation: Then, the modified reaction-diffusion Leslie-Gower prey-predator model is where d 1 , d 2 are diffusion coefficients. Considering the inhomogeneous distribution of species and the effect of delay τ (τ > 0) caused by predator gestation, model (1) with delay in the random network can be written as: and A i j is the adjacency matrix's element. N nodes make up the random network, and the nodes are connected with the connection probability p [24]. If τ = 0, the system (2) becomes: Clearly, the equilibrium points of system (3) are derived by the following equations: When a − δ β b > 0, we obtain only one positive equilibrium point E (u * , v * ): We assume that system (3) has only a positive equilibrium point E(u * , v * ) in the following, namely a − δ β b > 0 is always hold.

Stability analysis of the model without delay and diffusion
By linear stability analysis, the Jacobian matrix of (3) Then, the corresponding characteristic equation is where The roots of Eq. (5) are According to the Routh-Hurwitz Criterion, Eq. (5) has roots with negative real parts when tr J < 0 and det J > 0, namely, Reλ 1, And from Fig. 1, we can get that system (3) without diffusion is finally stabilized to the equilibrium point when condition (7) holds. (7) holds, for sys-

Lemma 1 Supposing that condition
In the following, we assume that tr J < 0 and det J > 0. Namely, condition (7) is satisfied.

Stability analysis of the network diffusion model without delay
For system (3), the eigenvalues Λ α and their cor- .., N , and Λ α ≤ 0 because L is a negative semidefinite real matrix. We assume that |Λ α | is in decreasing order, namely, Suppose Eq. (3) has the general solution of the following form: Substituting Eq. (8) into Eq. (3) give the following Thus, the characteristic equation of system (3) is where The roots of Eq. (10) are λ α = p 1 ± p 2 1 −4q 1 2 . Then, for Eq. (10), the real part of the roots is negative when p 1 < 0 and q 1 > 0. (7) holds, then for system (3), E(u * , v * ) is asymptotically stable if the diffusion coefficients d 1 and d 2 satisfy

Stability analysis of the predator-prey diffusion network with delay
For system (2), the equilibrium point is clearly E(u * , v * ). Here, we use two methods [31,32] to study the bifurcation dynamics property of system (2).

Method I
In this section, according to the method in [31] and regarding the delay τ as the bifurcation parameter, we Therefore, the corresponding characteristic equation is where When τ = 0, Eq. (12) becomes (10). (7) and (11) hold and regardless of what the network connection probability p is, then Eq. (12) has at least a simple pair of complex roots with zero real part ±iω α+ (ω α+ > 0) when τ = τ α j , where

Method II
In this section, according to the method in [32] and assuming the delay τ to be small, we consider the Turing instability of system (2).
Because τ is small enough, we replace v i (t − τ ) = v i (t) − τ dv i dt , and then, system (2) can be written as where Expanding (16) in the Taylor series and neglecting the higher order nonlinearities, we have According to the research method in Sect. 2.2, we obtain the corresponding characteristic equation of the system at E(u * , v * ) as where Because τ is small, 1 1+τ a 22 > 0, namely τ < τ 1 = − 1 a 22 . And let τ 2 = − tr J det J , we can calculate that τ 1 > τ 2 under condition (7). Supposing τ < τ 1 and condition (7) always hold in the following.

Numerical results
In this section, based on the above theoretical results, we present some numerical simulations. Here, we construct a random network with N = 100 nodes and consider the parameters as: For this case, the positive equilibrium point is E(0.1640, 0.2040), trJ = −0.201, det J = 0.1645, s 1 = 6.8386, and s 2 = 1.5078. These parameters satisfy condition (7). Therefore, according to Lemma 1, the system (3) without diffusion is asymptotically stable (Fig. 1). If we choose d 2 = 0.1, condition (11) holds, then the equilibrium E(u * , v * ) of system (3) is asymptotically stable, which is shown in Fig. 2.
For system (2), taking d 2 = 0.1 and p as randomly selected, we obtain that Λ α = −28.7563 and τ * = 0.8657. Since Λ α ≤ 0 is always true, and this satisfies Λ α ∈ (Λ α , 0], so the stability of equilibrium point E(u * , v * ) is only related to the delay τ no matter what p is. Therefore, considering delay τ as the bifurcation parameter and from Theorem 2, E(u * , v * ) is asymptotically stable when τ ∈ [0, τ * ), which is shown in Fig. 3a, b [Here we choose p = 0.1 or 0.6]. While the system becomes unstable when τ > τ * [Fig. 3c, d], where prey and predator will coexist and periodic oscillations will occur. And when τ reaches the critical value τ * , the system undergoes Hopf bifurcation when τ = τ * (Fig. 3e, f). We obtain the bifurcation diagram of u i about τ when p = 0.1 or 0.6, which is shown in Fig. 4, where the vertical axis represents the maximum and minimum values of u i . From those figures, we can get that system (2) would undergoes a transition from stable state to unstable state as the delay τ changes no matter what p is, and Hopf bifurcation occurs when τ = τ * .
When τ is small, we obtain τ 1 = 1.6667 and τ 2 = 1.2285. Considering d 2 and τ as the bifurcation parameters and from Theorem 3, E(u * , v * ) is asymptotically stable when τ < τ 2 and s 2 d 1 < d 2 < s 1 d 1 , which is shown in Fig. 5a when τ = 0.5 and d 2 = 0.1. In contrast, the system is unstable regardless of what p is when τ 2 < τ < τ 1 and s 2 d 1 < d 2 < s 1 d 1 , which is shown in Fig. 5b, c when τ = 1.5 and d 2 = 0.1 [Here we only choose p = 0.1 or 0.6]. And the bifurcation diagram of u i about τ is shown in Fig. 5d when p = 0.1 or 0.6, which indicates that system (2) is stable when τ is small, while the instability occurs when τ is beyond the critical bifurcation value. Due to the higher order nonlinearities are ignored for Taylor series expansion in the theoretical analysis, τ 2 is not equal to τ * , while τ 2 is equal to τ * in numerical simulation Figs. 4 and 5d.
Turing instability occurs when τ < τ 2 and d 2 > s 1 d 1 for existing Λ α ∈ [Λ α− , Λ α+ ]. When τ = 0.5 < τ 2 and d 2 = 0.5 > s 1 d 1 , Fig. 6a and b show the distribution of network eigenvalues Λ α , there exist some Λ α satisfy Λ α ∈ [Λ α− , Λ α+ ] for p = 0.1, while Λ α ∈ [Λ α− , Λ α+ ] is violated for p = 0.6. And Fig. 6c shows that Turing instability occurs when p = 0.1, while Fig. 6d shows that the system (2) is always stable when p = 0.6. For p = 0.1, Fig. 7a shows the system is stable when d 2 is small, and Turing instability occurs when d 2 is beyond the critical bifurcation value s 1 d 1 . While for p = 0.6, Fig. 7b shows the system is always stable, which means that no bifurcation occurs in system (2). Then, we get the approximate Turing bifurcation region of system (2) with respect to d 2 and p when τ = 0.5, which is shown in Fig. 8, where p approximately satisfies 1 N 2 ≤ p ≤ s 1 ln N N ( Due to the influence of network topology, there is no specific theoretical calculation methods about p at present, so we only give the approximate range of p based on the simulation results).

Conclusion
In this paper, a prey-predator diffusion network with delay is studied, and the stability of the positive equilibrium point E(u * , v * ) is investigated. From the linear stability analysis, we find that delay and diffusion can change the stability and produce Hopf and Turing bifurcations. Then, the numerical simulations are given.
First, for the network model without delay, we give the conditions of asymptotic stability, and the corresponding figures are shown in Figs. 1 and 2. Then, for the network system with delay, we use two methods to analyze the influences of delay and diffusion on the stability, Hopf and Turing bifurcations are obtained (Theorems 2 and 3). When condition (11) holds and regardless of the network connection probability p, the network system is always stable if τ < τ * or τ < τ 2 , while the system is unstable if τ > τ * or τ 2 < τ < τ 1 and undergoes Hopf bifurcation if τ = τ * (Figs. 3, 5ac), the corresponding bifurcation diagrams about τ are shown in Figs. 4 and 5d. Namely, regardless of how the network nodes are connected, prey and predators will coexist, and the system will oscillate periodically.
Then, if condition (11) is not held (namely when d 2 > s 1 d 1 ), Turing instability occurs when τ < τ 2 and p approximately satisfies 1 N 2 ≤ p ≤ s 1 ln N N (Λ α ∈ [Λ α− , Λ α+ ]), the numerical simulations about distribution of network eigenvalues and the time series diagram as well as the bifurcation diagrams are shown in Figs. 6, 7, and 8. Therefor, we can get that Turing instability is affected by the network structure. From those conclusions and figures, it is not difficult to find that Turing and Hopf bifurcation exist in real delay ecological networks. And through the theoretical analysis methods in this paper, we can qualitatively control some conditions in reality to achieve the results we want, such as protecting endangered species, eliminating pests and maintaining ecological balance. However, since the specific relationship between p and the eigenvalues of L is not obtained, and the region of Turing instability is approximate, there are some problems that we have not studied, such as the application of real data, the stability and period as well as the direction of the bifurcating periodic solution. Thus, there is a lot of work to do.