Stability switches in a ring-structured predator–prey metapopulation model with dispersal delay

In this paper, we consider a predator–prey metapopulation model with a ring-structured configuration of an arbitrary and finite number of patches. The prey are assumed to disperse between the connected patches with a constant dispersal delay. We show that the dispersal delay can induce stability switches exhibiting both stabilizing and destabilizing roles in the stability of the symmetric coexistence equilibrium. Numerical simulations are presented to further illustrate the effects of the dispersal delay, the dispersal rate, the fraction of dispersal due to predation avoidance and the network topology on the number of stability switches.


Introduction
The predator-prey relationship is one of the fundamental relationships in ecological systems and has been extensively studied in the literature [1]. Recently, there has been a growing interest in integrating the spatial heterogeneity into modeling metapopulation dynamics. A metacommunity is a set of local communities linked by dispersal of species. It has been shown that metacommunity dynamics can be greatly influenced by dispersal between discrete patches [2,3].
Many factors such as lack of food, competition, sex, age, climate, seasons and overpopulation in a patch contribute to the dispersal of either prey, or predator or both between patches [4]. For instance, prey may choose to move to other patches on the basis of resource availability and predation risk, while predator tend to move to patches with more prey. Many mathematical models have been proposed to study the impacts of dispersal over patches on metacommunity dynamics. For example, the dispersal of a single species is considered in [5][6][7][8][9], while the case with the dispersal of both prey and predator is studied in [10][11][12][13][14][15]. It has been shown that the dispersal plays an important role on the persistence and stability of the resulting ecosystems.
Since it always takes time for species to complete the dispersal from one patch to another, it is natural to incorporate the travel time into modeling. Along this line, dispersal delay is considered in [16][17][18][19] and it is shown that the dispersal delay may induce very complex dynamics. For instance, under the framework of a two-patch predator-prey model with delayed dispersal of predator, it is found that if the dispersal rates are in an intermediate range and the mean time that the predator spent in one patch is much shorter than the timescale of reproduction of the prey and is larger than the double mean time of capture of prey, the dispersal delay can induce stability switches such that an otherwise unstable coexistence equilibrium can be stabilized over a finite number of stability intervals [19].
Most of the aforementioned work considered the case that dispersal occurs only between two patches. If the environment involves more than two patches, then there could be many different configurations of network topology. How the impacts of dispersal differ under different configurations of network topology has not been well addressed in the literature [20]. Recently, Mai et al. considered a predator-prey metapopulation model over an arbitrary and finite number of identical patches [18], where the dispersal is assumed to occur among all patches in the environment with the same rates. Such figuration is referred to as the fully connected configuration. In this paper, we shall consider the ringstructured configuration: the patches under consideration form a ring structure, and dispersal occurs only between two neighbouring patches. We allow the number of patches to be arbitrary and finite. In order to have a comparison with the results with the fully connected configuration, we also assume that all patches are identical and the dispersal rates of the prey are the same.
The rest of the paper is organized as follows. In Sect. 2 we present the metapopulation model. The stability analysis of the symmetric coexistence equilibrium is given in Sect. 3. A specific example is presented in Sect. 4. Some numerical simulations based on our theoretical analysis are reported in Sect. 5. We summarize our results in Sect. 6.

The model
We use a general predator-prey model given in [21] to describe the interaction between prey and predator in an isolated patch, (2.1) Here x represents the density of prey, and y denotes the density of predator. In the absence of predator, prey have the per capita growth rate g(x). p(x) is the functional response, and ε is the conversion ratio of loss of prey due to predation. q(x) denotes the death rate (which may depend on the prey density) of predator. We assume g(x), p(x) and q(x) satisfy the following assumptions: (A1) g(0) > 0 and there exists K > 0 such that g(K) = 0 and (x - Typical functions satisfying the above assumptions include which are widely used in the literature [21,22]. It then follows from the assumptions (A2) and (A3) that there exists a unique positive number, which we denote by x * , such that p(x * ) = q(x * ). If further, x * < K , then (A1) implies that y * = x * g(x * )/εp(x * ) > 0 and (x * , y * ) is the unique coexistence equilibrium of system (2.1).
For simplicity, we assume that all patches under the ring-structured configuration are identical and the dynamics in each patch is governed by (2.1). In addition, only prey are assumed to disperse from its habitat patch to its two neighboring patches with a dispersal rate d, while the predator does not move between patches. As in [17], we assume that the dispersal of prey is due to two factors: a random effect and predation avoidance, and we use α ∈ [0, 1] to denote the fraction of prey dispersal due to predation avoidance. Let τ denote the finite travel time during the process of dispersal. Then we can use the following metapopulation model to describe the interactions of predator and prey under a ringstructured configuration: for i = 1, 2, . . . , n, n ≥ 3, where γ is a scaling factor proportional to ε. To close the ring, we define x 0 = x n , x n+1 = x 1 , y 0 = y n and y n+1 = y 1 .

Stability of the symmetric coexistence equilibrium
Clearly, under our assumptions, system (2.2) always admits a symmetric coexistence equilibrium E * n = (x * , y * , . . . , x * , y * ) ∈ R 2n provided that x * ∈ (0, K). In this section, we are concerned with the linear stability analysis of the coexistence equilibrium E * n .

Instantaneous dispersal: τ = 0
We first consider the special case where the dispersal of prey is instantaneous, i.e., τ = 0. Linearizing system (2.2) at the equilibrium E * n , we obtain the associated characteristic equation det J n = 0, where the 2n × 2n matrix J n is given by It follows from (A2)-(A3) that B > 0, D > 0 and μ > 0. Therefore, J n is a block-circulant matrix with 2 × 2 blocks. Using the formula in [23], we find that Thus, the eigenvalues are determined by the equations det J 1 + 2J 2 · cos 2kπ n , k = 0, 1, . . . , n -1. In the sequel, we always assume that k ∈ K. When τ = 0, the characteristic equations reduce to By the Routh-Hurwitz criterion [24], it is easy to obtain the following result. (1) If A < 0, the coexistence equilibrium E * n is locally asymptotically stable.
Proof If A < 0, then A -2db(1 -cos 2kπ n ) < 0. Therefore, all characteristic roots of (3.3) have negative real parts and hence E * n is locally asymptotically stable. If A > 0, then (3.3) with k = 0 admits two roots with positive real parts. Thus E * n is unstable.
Remark 3.2 Note that the stability condition of the coexistence equilibrium of the ringstructured system is the same as that of the coexistence equilibrium of the single-patch system. This indicates that the instantaneous dispersal of prey does not affect the stability of the coexistence equilibrium.

Random dispersal with delay: τ > 0 and α = 0
Linearizing system (2.2) at the equilibrium E * n , we obtain the characteristic equations of system (2.2): where k ∈ K. When α = 0, the characteristic equations (3.4) become It is well known that the coexistence equilibrium E * n of system (2.2) is stable if and only if all characteristic roots of (3.5) have negative real parts [25]. To explore how the dispersal delay affects the dynamics of the coexistence equilibrium, next we use the dispersal delay τ as the bifurcation parameter. As τ increases, the stability of the equilibrium changes only when some characteristic roots cross the imaginary axis transversely [25][26][27]. It is clear that zero is not a characteristic root of (3.5), we thus look for purely imaginary roots. Since complex roots appear in pairs, we assume that λ = iω, with ω > 0. Substituting λ = iω into (3.5), and separating the real and imaginary parts, we obtain 2d -A = 2d cos ωτ · cos 2kπ n , -ω 2 + εμB = 2ωd sin ωτ · cos 2kπ n . (3.6) If C k > 0, then (3.7) admits the following two positive roots:
To analyze the stability of the system (2.2) at E * n , next we consider two cases: A < 0 and A > 0, respectively.
A < 0 implies that C k < 0 for all k ∈ K. Thus (3.7) does not admit any solutions. That is, no characteristic roots of (3.4) would cross the imaginary axis. By Lemma 3.1, E * n is stable with τ = 0, thus, E * n remains stable for τ > 0. Case 2: A > 0. In this case, the coexistence equilibrium E * n is unstable when τ = 0 by Lemma 3.1. Next we examine if the unstable coexistence equilibrium E * n still remains unstable as τ increases, or can it be stabilized for some values of τ . To this end, we consider three cases: In this case, it follows from (3.7) that C k < 0 for all k ∈ K. Therefore, there are no such purely imaginary roots. Consequently, the coexistence equilibrium E * n remains unstable as τ increases. ( For such k, (3.7) admits two positive solutions: ω 1k and ω 2k , given in (3.8). Therefore, we distinguish three different conditions: We denote by n(M) the number of elements in the set M. Let Under condition (1), there are n(K 1 ) characteristic equations (3.5). For these characteristic equations, together with (3.6), we have -1 < cos ωτ < 0. Then Thus, sin ω 1k τ > 0 and sin ω 2k τ < 0.
We define Then θ (k) 1 ∈ ( π 2 , π) and θ (k) 2 ∈ (π, 3π 2 ). This shows that purely imaginary roots of the characteristic equations (3.5) under condition (1) are obtained at the following sequences of positive values of τ : 3) that, as τ increases, the characteristic roots cross the imaginary axis through ±iω 1k at τ = τ (k) l,1 from right to left and the number of characteristic roots with positive real parts is reduced by 2. The characteristic roots cross the imaginary axis through ±iω 2k at τ = τ (k) l,2 from left to right and the number of characteristic roots with positive real parts is increased by 2.
Similarly, under condition (2), there are n(K 2 ) characteristic equations (3.5) satisfying cos 2kπ n < 2d-A 2d < 0. For these characteristic equations, in a similar manner, we can show that there are purely imaginary roots at the following sequences of (Hopf bifurcation) values of τ : where k ∈ K 2 , l = 0, 1, 2, . . . . By Lemma 3.3, as τ increases, the characteristic roots cross the imaginary axis through ±iω 1k at τ = τ (k) l,3 from right to left and the number of characteristic roots with positive real parts is decreased by 2, while the crossing through ±iω 2k occurs at τ = τ (k) l,4 from left to right, and the number of characteristic roots with positive real parts is increased by 2.
From the above analysis, we see that, for case (ii), n(K 3 ) characteristic equations have no purely imaginary roots; while for n(K 1 ) + n(K 2 ) characteristic equations, there are sequences of critical values of τ at which the characteristic roots cross the imaginary axis (either from left to right or from right to left). For case (iii), n(K 6 ) characteristic equations have no purely imaginary roots; while for n(K 4 ) + n(K 5 ) characteristic equations, there are sequences of critical values of τ at which the characteristic roots cross the imaginary axis.
To further investigate the possibility of undergoing stability switches, we introduce some notation. For the above obtained sequences of critical values of τ , we define Γ j = τ (k) 0,j , τ (k) 1,j , . . . , τ (k) k,j , . . . , j = 1, 2, 3, 4, 5, 6, 7, 8, where k ∈ K 1 if j = 1, 2; k ∈ K 2 if j = 3, 4; k ∈ K 4 if j = 5, 6 and k ∈ K 5 if j = 7, 8. For each case (case (ii) and case (iii)), we sort the obtained sequences of τ as a single increasing sequence with τ j < τ j+1 , j = 1, 2, . . . . For instance, if n(K 1 ) = 1, n(K 2 ) = 1 and τ (2) 0,4 < τ (0) 0,1 < τ (0) 0,2 < τ (2) 0,3 < · · · , then the new sequence reads We associate each element τ j of the sequence with a number σ j defined by At τ = τ j , if σ j = +1, then it means a pair of purely imaginary roots cross the imaginary axis at τ = τ j from left to right, and the number of characteristic roots with positive real parts increases by 2 and we say the delay τ j has a positive jump, otherwise, if σ j = -1, then a pair of purely imaginary roots cross the imaginary axis at τ = τ j from right to left, and the number of characteristic roots with positive real parts is reduced by 2 and we say the delay τ j has a negative jump. Define a sequence {s j } ∞ j=1 as Note that characteristic equations (3.5) with τ = 0 reduce to Here, we denote by n(Re(λ) > 0) the number of characteristic roots with positive real parts of characteristic equation (3.9).
We denote by N the number of elements in the set S j , j = 2 or j = 3, i.e., N = n(S j ), j = 2, 3. (3.12) Then N determines the number of stability intervals.

Predator-induced dispersal with delay: τ > 0 and α > 0
In this section, we study the stability of the coexistence equilibrium E * n of system (2.2) with α > 0 and τ > 0.
For any fixed α > 0, we use the same method as for the case with α = 0 to find sequences of critical values of τ at which the characteristic equations (3.4) admit purely imaginary characteristic roots, and to compute directions of the characteristic roots crossing the imaginary axis and then to determine if there are stability intervals.
Following the results in Sect. 5 of Cooke and Grossman [26], the characteristic equation (3.4) admits purely imaginary roots if and only if the following inequalities hold: (3.13) can never hold, thus, all characteristic equations (3.4) have no purely imaginary roots. If d > A 2 4(ã+Ab) andC k > 0, then it is possible for characteristic equations (3.4) with some k ∈ K to admit purely imaginary roots and thus the interplay of α and τ may affect the stability of coexistence equilibrium. Substituting λ = iω (withω > 0) into (3.4) and separating the real and imaginary parts, we obtain , sinωτ =ω˜b (εμB-ω 2 )+ãA 2d cos 2kπ n (b 2ω2 +ã 2 ) .

An example
In this section, we give an example to illustrate our theoretical analysis in detail. Take n = 4. Then the characteristic equations (3.4) of system (2.2) reduce to where k = 0, 1, 2.

Numerical simulations
In this section, we present some numerical simulations to demonstrate our theoretical results. We consider the following ring-structured patchy model:  and we have an N = 1 stability interval. The coexistence equilibrium E * n is locally asymptotically stable for τ ∈ (τ 1 , τ 2 ) = (3.50, 4.73). A bifurcation diagram is presented in Fig. 1.
Next we fix α = 0.5 and vary the value of the dispersal rate d. We find that the number of stability intervals changes from 0 to 1 near d = 0.0421, then it changes to 2 when d = 0.0721, and back to 1 at about d = 0.2651. This shows that the dispersal rate can also induce finite number of stability switches. This is illustrated in Fig. 2 . Now we fix d = 0.13 and vary the value of α. We find that the number of stability intervals changes from 1 to 2 and back to 1 near α = 0.023 and α = 0.939, respectively (see Fig. 3). This indicates that the fraction of prey dispersing due to predation avoidance can also induce stability switches.
For the case with α > 0 and A < 0, we take parameter values κ = 1.8, μ = 0.3, ε = 1.1, n = 4, α = 0.6, d = 0.6, γ = 0.5. Regarding τ as the bifurcation parameter, we find that there are N = 30 stability intervals. We show some of them in Fig. 4. the right panel: the distribution of stability intervals of system (5.1). Other parameter values are the same as in Fig. 2 Next, using Figs. 5 and 6, we illustrate that the configuration of network topology also has an impact on the stability switches of the coexistence equilibrium by comparing our ringstructured patchy model with the fully connected model studied in [18]. The differences are significant when n is odd, and when α is small and d is large. The parameter values are κ = 2, μ = 0.2, ε = 1, γ = 0.5

Summary
In this paper, we have considered a metapopulation predator-prey model with delayed dispersal within a ring-structured configuration of network topology. Prey in one patch are assumed to move between two neighboring patches on the ring. We have shown that the dispersal delay can induce stability switches for the case with random dispersal only and also for the case with random dispersal and density-dependent dispersal. We have also shown that the dispersal rate (d) and the fraction of dispersal due to predation avoidance (α) can induce stability switches. In addition, by comparing our ring-structured patchy model with the fully connected model studied in [18], we find that the configuration of network topology also has an impact on the stability switches of the coexistence equilib- Figure 6 The number of stability intervals N against α, d and n for the fully connected patch model in [18].
The parameter values are κ = 2, μ = 0.2, ε = 1, γ = 0.5 rium, especially when the number of patches n is odd, the fraction α of dispersal due to predation avoidance is small and the dispersal rate d is large.