Bifurcation analysis in a diffusion mussel-algae interaction system with delays considering the half-saturation constant

In this paper, the kinetics of a class of delayed reaction-diffusion mussel-algae system under Neumann boundary conditions with the half-saturation constant is studied. The global existence and priori bounds as well as the existence conditions of positive equilibrium are obtained. The half-saturation constant affects the stability of the system and may result in Turing instability. When the half-saturation constant exceeds a certain critical value, the boundary equilibrium is globally asymptotically stable which means that the larger half-saturation constant forces the mussel population toward extinction. By analyzing the distribution of the roots of the characteristic equation with two delays, the stability conditions of the positive equilibrium in the parameter space are obtained. The stability of the positive equilibrium can be changed by steady-state bifurcation, Hopf bifurcation, Hopf-Hopf bifurcation or Hopf-steady state bifurcation, which can be verified by some numerical simulations. Among the parameters, the half-saturation constant and two delays drive the complexity of system dynamics.


Introduction
There has been great interest in dynamical characteristics (including stable, unstable, bifurcation and oscillatory behavior) of population models since Vito Volterra and James Lotka proposed the seminal models of predator-prey models in the mid-1920s. The musselalgal interaction model is one of them. Mussel beds are a typical system for studying pattern formation, in which patterns develop at two different scales, namely large scale banded patterns and small scale reticulated patterns. Recently, the emergence of nonlinear traveling waves in such systems has led to some researchers to propose this mechanism to better understand the periodic patterns observed in young beds of mussels feeding on algae. Many ecologists and mathematicians had used the mussel-algal model to study the development of mussel beds. Particularly, van de Koppel et al. [1] reported that these river beds occurred on the soft substrates of tidal pools at the edge of the Wadden Sea. Their model included a positive feedback to describe these promoting effects at higher mussel densities. From these observations, van de Koppel et al. [2] and Liu et al. [3] devised a field experiment containing a young mussel homogeneous substrate covered by a relatively stationary seawater layer including algae as a food source and observed various clumps or gapped type patterns. These laboratory experiments, consisting of young mussel beds, were initially distributed evenly over a homogeneous substrate covered with a static ocean layer containing algae, replicating patterns observed in the field. Van de Koppel et al. [1] proposed a coupled partial differential equation with two components containing the diffusion effect of mussel and advection effect of algae, but ignored the lateral diffusion effect of latter, which provided an explanation for the pattern formation of mussel. Based on the nonlinear numerical methods, Wang et al. [4] performed stability analyses for linear systems, compared numerical results for nonlinear systems, and focused on the sensitivity of model parameters to evaluate their effects on the velocity, amplitude, and wavelength of banded migration modes. Liu et al. [5] extended the simulation to two dimensions. Cangelosi et al. [6] modified the unidirectional advection formula of algal concentration into random Brownian dispersion, and obtained a Turing patterns similar to that of the original algal model. The system is given as follows: where A(x, t) and M(x, t) are the algae concentrations in the lower water layer overlying the mussel bed and the mussel biomass density on the sediment at time t and space x ∈ Ω. Ω ⊆ R n is a bounded domain with a smooth boundary ∂Ω. Δ is Laplace operator and ρ is the exchange rate between the upper and lower water layers. A u describes the uniform concentration of algae in the upper water body. The velocity V of the tidal current is assumed to be acting in the positive direction of x. c is the consumption constant and H is the height of the lower water layer. e is the biomass conversion constant for mussel to ingesting algae. The maximal mussel mortality rate per capita d M represents mussel loss due to predation or dislodgement. k M is the value of M(x, t) at which mortality is half maximal. D M is a measure of mussel motility caused by their pedal locomotion mechanism and its concomitant byssal thread deployment [4]. The diffusion of algae consists of two parts: an active lateral diffusion D A because their mech-anism of locomotion is flagella and a passive vertical diffusion ρ on account of the photo-gyro-gravitactic driven exchange between the upper and lower benthic water layers [7] (see Fig. 1). Note that although on the large scale, the algae was considered advection and tidal flows, on the small scale they actually dispersed in the liquid as Brownian particles. The experiments had shown that mussels could move actively within and between clusters, meaning that advection and tidal currents in small scales had little effect on mussel beds. Although advection and diffusion were two different ecological processes, in real mussel bed ecosystems, the two processes were usually coexisting and had the same activation-inhibition mechanism [8].
The advection and diffusion are equivalent for the emergence of spatial self-organizing patterns. Song et al. [9] and Jiang [10] studied Turing-Hopf bifurcations of system (1) under V = 0 and Neumann boundary conditions based on the normal form method, and obtained explicit dynamic classification at the critical point. Based on system (1), Shen and Wei [11,12] considered the effects of delay. Shen and Wei [13] considered that the mortality involves a positive feedback term resulting from the reduction of dislodgment and predation and a negative feedback term resulting from the intraspecific competition for mussel. Their results suggested that the regular patterning in mussel beds were caused by the high mobility of algae or the low diffusion of mussels. Zhong et al. [14] investigated the spatiotemporal dynamics for a diffusive musselalgae model near a Hopf bifurcation point. They found that the strip patterns were mainly induced by Turing instability in equilibrium and spot patterns were mainly induced by Turing instability in limit cycles by numerical simulations. In this paper, we use the suggestion in [15] to introduce saturated growth in the mussel population and revise system (1) as where β is the half-saturation constant. Fig. 1 Schematic for the mussel-algae interaction-diffusion system [6] By introducing the dimensionless variables: In 1986, Stepan carried out the first Hopf analysis in the presence of a single delay and the possibility of Hopf-Hopf bifurcation appeared [16]. It is well known that the delay may lead to periodic oscillations and the diffusion may cause Turing instability [17][18][19]. Inspired by the work of [10][11][12]15], we establish a delayed mussel-algae system with the Neumann boundary conditions in one-dimensional space as follows: where τ = max{τ 1 , τ 2 }, τ 1 represents the period of digestion of mussel, i.e. the consumption of algae in earlier times can make the mussel population in a later time in order to be more practical. τ 2 represents the mortality of mussel depending on the state whether they have be eaten in the past. A system with two time delays can produce complex dynamics [20,21]. This paper will try to explore the complex bifurcating phenomena of system (4). Compared with the work in [10][11][12]15], this article has two innovations: firstly, the half-saturation term is introduced into the mussel population growth function and the effect to the dynamics of system is investigated; secondly, two different delays are introduced and the bifurcation problem caused by delays is discussed.
Define the Sobolev space with real-value and its complexification X = X iX = {α + iβ| α, β ∈ X} with L 2 inner product defined as The structure of this paper is as follows. In next section, the existence conditions of the positive equilibrium are given. When the two delays do not exist, we get the global existence and a priori bound of solutions, at the same time, the conditions of Turing instability are given in Sect. 3. Furthermore, for the delayed system, when the half-saturation constant exceeds some critical value, the mussel population will die out. By analyzing the distribution of eigenvalues, we obtain the stability conditions for the positive equilibrium in the parameter space. The positive equilibrium can lose the stability by the steady-state bifurcation, Hopf bifurcation, Hopf-Hopf bifurcation or Hopf-steady state bifurcation. These results are shown in Sect. 4. At last, some numerical simulations verify that the dynamic behaviors of the system caused by the half-saturation constant are consistent with the theoretical findings and some new findings are foreseen.

Equilibrium
The equilibrium E(a, m) of system (4) satisfies the following equations: The boundary equilibrium E 0 (1, 0) always exists. The existence of positive equilibrium E * (a * , m * ) needs some analyses described as follows. From Eq. (5), it obtains that a must satisfy the following equation: where has a root a * satisfying 0 < a * < 1, then system (4) exists a positive equilibrium E * (a * , m * ), where The existence of positive root of Eq. (6) is given as follows.

Dynamics of system (3)
In this section, we focus on system (3) to investigate the property of solutions and Turing instability.

Existence and priori bound of solutions
In this part, the sufficient conditions for the existence of a positive solution of system (3) are given and a priori bound of solutions is derived.
The strong maximum principle implies that a(x, t) > 0 and m(x, t) > 0 when t > 0 for all x ∈ [0, lπ ]. This completes the proof of (i). It has known that a(x, t) ≤ a * (t) for t > 0, and a * (t) is the unique solution of to prove the boundedness of m(x, t), it only needs to prove that m * (t) is bounded since m * (t) is a uppersolution of m(x, t). According to a(x, t) . Since the randomicity . This completes the proof of (ii).

Stability of the nonspatial system
In this subsection, we discuss the dynamics of the following nonspatial system: The linearization of system (11) is For the equilibrium E 0 , it has b 11 and J 0 = α r (b+1−γ ) r (b+1) , which means that E 0 is locally asymptotically stable for b + 1 > γ and unstable for b + 1 < γ .
α and a ≥ 1, then E * is locally asymptotically stable. If b < 1 α and a < 1, then E * is locally asymptotically stable for 0 < a * < a and unstable for a * > a . Remark 1 2, the parameter b has directly effect to the stability of E * .
Taking r = 0.5, μ = 0.6, α = 0.01, we can draw the (b, γ ) plane plot in where the stability regions of E * are given (see Fig. 2A). In the region II, Hopf bifurcation line is given in Fig. 2B. We can choose γ = 10 and obtain b = 2.727 and 6.7151, the system exists stable periodic solutions and the phase plots are shown in Fig. 2C, D.
Obviously, r does not affect the existence of E * , but it does affect its stability. Solving T 0 = 0 for r , it has We have the following results on the stability and Hopf bifurcation of system (11) at E * .

Turing instability
In this subsection, we take the effect of diffusion into account on the dynamical behavior of system (3) and investigate the diffusion-driven instability for where where b i j (i, j = 1, 2) are given as (12). The characteristic equation of system (14) constrained by Neumann boundary conditions is where k is the wave number, If E * is stable in the absence of diffusion, but becomes unstable in the opposite case, we call this as Turing instability. Therefore, Turing instability only occurs when there is at least one positive integer k such that J k < 0. We get the following result.
Let k 2 min be the minimum value of p(k 2 ). Thus if μ < μ * , then If p(k 2 min ) < 0, the sufficient condition for the system to be unstable reduces to (3) under the following conditions:

Theorem 5 Turing instability occurs for system
At the same time, let the two roots of p(k 2 ) be k 2 − and k 2 + , then k 2 ∈ (k 2 − , k 2 + ) limits the range of Turing instability for a locally stable steady state.

Spatiotemporal dynamics of system (4)
Let E(a, m) be any equilibrium, define u = a − a and v = m − m, then system (4) becomes whose characteristic equation is where

Stability of boundary equilibrium
For equilibrium E 0 (1, 0) it has l 21 = 0, l 22 = 0 and let In the following, we prove the globally asymptotically stability of E 0 (1, 0) for system (4).
Since ε the sufficiently small and Corollary 2.1 of [26], we can obtain the conclusion. This completes the proof.

Stability of the positive equilibrium and bifurcation
From Theorem 6, it knows that E 0 (1, 0) is globally asymptotically stable when b + 1 > γ . Hence, E * is unstable when b+1 > γ which is contained in the conditions (A1) and (A3). Hence, we only need to investigate the stability of the positive equilibrium under the condition (A2). In this section, it uses E * (a * , m * ) instead of E * − (a * , m * ) as the positive equilibrium of system (4) and the parameters b, α, γ satisfy (A2).
At E * , system (15) becomes whose characteristic equation is where The stability/ instability of E * can be determined by the root distribution of the characteristic equation (20) for = +∞. If there exists some k satisfying q k < −1, then D k (0, τ 1 , τ 2 ) < 0. Hence, E * is unstable. In the following, we assume q k ≥ −1 to investigate the stability of E * . In addition, under the condition (H1), λ = 0 is not the root of Eq. (20).
In the absence of delay, we firstly seek for the conditions at which E * is stable for τ 1 = τ 2 = 0. Taking and if k 0 is just right an integer, then J 2k = 0 for k = k 0 , otherwise, Particularly, if b (b+1) 2 ≥ α holds, then all roots of Eq. (21) have negative real parts for any k ∈ N and any γ > 0. The following Lemma is from [27]. r, μ, γ, α, b, τ 1 and τ 2 change in the parameter space, the number (counting multiplicities) of eigenvalues of Re(λ) > 0 changes only when an eigenvalue passes through the imaginary axis in the complex plane.

Lemma 1 When
From lemma 1, if all roots of the characteristic equation (20) have negative real parts, then the trivial solution of system (19) will be asymptotically stable. The trivial solution of system (19) can only be unstable as the root of the characteristic equation changes through the imaginary axis. Therefore, the stability changes in parameter space at points where Eq. (20) has a root with zero real part.
Firstly, we consider the zero root of Eq. (20), λ = 0, which will occur when q k = −1 for some k < k 0 and λ = 0 is not the root for any k ≥ k 0 .
Furthermore, Hence Thus for some k < k 0 , the root λ = 0 is simple and passes through the zero with nonzero speed everywhere on the line q k = −1 except for one possibly isolated point given by Hence, for q k = −1, Eq. (20) has a simple zero root when J 1k + J 2k > 0, k < k 0 and τ 1 = τ * 1k . However, when q k = −1 and τ 1 = τ * 1k , which means that λ = 0 is the double zero root of Eq. (20). We summarize as follows.
Locating the pure imaginary roots is slightly more complicated. In the following, it assumes that (H1) hold. Firstly, let τ 2 = 0, then Eq. (20) becomes We investigate the existence of a pure imaginary root for Eq. (22). The analysis of Eq. (22) is based on the conclusion that the number of characteristic roots with positive real parts will change only when the characteristic roots enter the complex plane through the imaginary axis as the change of parameters. {τ 0 1k }. When it crosses a τ n 1k , the number of positive real parts root increases by 2 as the increase of n for fixed Q in (Q, τ 1 ) plane, where τ n 1k is shown in the equality (24) and Proof Without loss of generality, for a fixed k ∈ N, let iw be the root of the (k + 1)th equation of the equation (22) and separate the real and imaginary parts, it has and sin wτ 1 > 0. Furthermore, Squaring and adding the both sides of equations (23), it has which yields to From Eq. (25), if q 0 ≥ 1, then q k ≥ 1 and Eq. (25) has no any positive root for any k. On the contrary, if −1 < q k < 1, Eq. (25) has a uniquely positive root w k + . Thus, if k ∈ S 1 , then Eq. (22) exists purely imaginary roots as long as τ 1 takes the critical values determined by the equality (24).

Theorem 10 If
Remark 2 In fact, choosing τ 2 as bifurcating parameter for a fixed τ 1 , the stability switches may occur if F k (ω) = 0 has more than a positive root. This fact can be shown in Fig. 14 in Sect. 5.
In general, we can eliminate ω from Eq. (27) and get an equation for l, k, r, α, b, γ, μ, τ 1 and τ 2 , which defines a hypersurface in the parameter space. Since this cannot be done in practice, to describe this hypersurface we can represent Q and τ 1 as functions of ω. Firstly, for fixed k, we find Q by squaring and adding in both sides of Eq. (27): Taking the ratio of Eq. (27) and solving for τ 1 yields an expression involving the inverse tangent function.
Noting that the sign of cos ωτ 1 may be obtained from Eq. (27), we find that the appropriate expression for τ 1 is for j ∈ N and artan denotes the inverse tangent function which has range (−π/2, π/2). Clearly, the equality (31) represents an infinite curve family. We take notice of the following limits, This value τ * 1k is feasible only if k < k 0 .
In summary, the line q k = −1 represents the set of points in parameter space for which Eq. (20) has a zero eigenvalue. The curves (with multiple branches) defined by the equalities (30) Here τ * 1k > 0 holds because of k < k 0 . This completes the proof. Proof The λ = iω curves defined by the equalities (30) and (31) give the points in the (Q, τ 1 )-plane, where the real part of λ is zero. We can now prove our claim by considering the appropriate derivatives. Differentiating the equalities (31) with respect to ω, we get Taking the derivative of Eq. (20) with respect to Q and substituting λ = iω leads to the following expression: where dω > 0, then the real part of λ becomes negative, which completes the proof.
Remark 3 Lemma 1 indicates that with fixed k, α, b, μ, r, γ, τ 1 and τ 2 , a completely stable region can be obtained from these subsets by increasing until the λ = 0 line or λ = iω curve is reached. Theorems 11 and 12 describe how the real part of the root of the characteristic equation changes when the boundary of the stable region is crossed.
For fixed k, r, γ, α, b, μ and τ 2 , the equalities (30) and (31) describe the curves which lie in the right of the (Q, τ 1 ) plane, and are parameterized by ω. However, Q c k and Q s k is quite complicated functions of ω and the system parameters. Nevertheless, we can still make some assertions for the stability of E * . For fixed l, k, r, γ, α, b, μ and τ 2  Proof Consider first intersections of the curves defined by the equalities (30) and (31). Since these curves are defined according to ω parameter, the intersection of two different curves occurs when Q and τ 1 have the same values. If Q k H (ω) is monotone on ω then this is impossible. Now consider the curves defined by the equalities (30) and (31). Consideration lim

Theorem 13
is increasing about ω, then there are no intersection points of the line q k = −1 with these curves. This completes the proof.
Since the complexity of Q k H as the function of ω, we only can obtain the following monotonicity.
Hence there is no region for which all the roots of Eq. (20) have negative real parts for all positive values of τ 1 and τ 2 . This completes the proof.
In the following, we consider the behavior of system (4) near the λ = iω curve. We have proved that the characteristic equation has a pair of complex conjugate imaginary roots on this curve, so we can expect system (4) to exhibit a Hopf bifurcation on this curve. We firstly establish the usual nonresonant conditions.

Lemma 2
If λ c = iω c is a purely imaginary root of the characteristic equation (20), then λ c is simple and any root λ other than λ c and λ c satisfies λ = mλ c for any integer m.
Proof To prove that λ c = iω c is simple, we need to show that D k (λ c , τ 1 , τ 2 ) = 0, where Suppose that D k (λ c , τ 1 , τ 2 ) = 0. Substituting λ c = iω c , separating the real and imaginary parts, we get Consider the ratio of two equations in the equations (37), Furthermore, ω c must satisfy the ratio of two equations in equations (27) when Q = Q k H and λ = iω c . Furthermore, it can obtain we notice that H ω c is identical to the numerator of Eq. (35). Hence, by Theorem 9, H ω c > 0 for all τ 2 < τ (1) 2k . Moreover, zeros of H ω c occur when dQ k H (ω c )/dω = 0. Therefore, excluding the values of ω c satisfying dQ k H (ω c )/dω = 0, λ c = iω c is a simple root of Eq . (20). (ω) are defined by (30) and (31). Assume that ω c is neither a zero of dQ k H /dω nor of dτ j 1k /dω, then system (4) undergoes a Hopf bifurcation at Q = Q kc H . We conclude that if system (4) satisfies the appropriate non-resonance conditions, the curves (30) and (31) define the Hopf bifurcation lines on the (Q, τ 1 ) plane, and the vertical line q k = −1 defines the steadystate bifurcation. One way to determine the bifurcation properties of Hopf is to analyze the central manifold of the equation using the method in [28]. To determine the type of steady-state bifurcation in the system, we need to analyze the number and stability of non-trivial equilibrium points. To do this, we need to linearize the system about these equilibria, generate a new characteristic equation, and analyze its roots. This is a big task, so we omitted this process due to space constraints.
Codimension-2 bifurcation occurs when two of these different bifurcations happen simultaneously. Hopf-Hopf bifurcation point exists when the characteristic equation has two pairs of pure imaginary roots ±iω 1 and ±iω 2 . For the considered system, these points can occur for any set of the parameters l, k, r, α, b, γ, μ, τ 2 such that Q k H (ω 1 ) = Q k H (ω 2 ) and τ j 1k (ω 1 ) = τ jj 1k (ω 2 ), for some j, jj ∈ N. In general, these points cannot be solved in a closed form. However, they can be computed numerically (see Fig. 11A). Hopf-steady state interactions exist when the characteristic equations has both a zero root and a pure imaginary pair. For system (4), these points may be found by solving Q(ω) = −J 1k J 2k for ω and substituting in τ j 1k as given by the equality (31) where they appear as intersection points of the Hopf bifurcation curves with the line Q(ω) = −J 1k J 2k (see Fig. 12B). From Theorem 9, it is clear that neither type of codimension-2 point can occur in system (4) if τ 2 < τ (1) 2k . Codimension-2 points can be the source of more complicated dynamics such as multistability and quasiperiodicity.

Numerical simulation
Firstly, it chooses the same parameters as [11] except for l and b for system (4). i.e., It can obtain that the boundary equilibrium E 0 (1, 0) is globally asymptotically stable for system (4) (see Fig. 4).
Authors [11] considered the linear growth in mussel population, when τ = 3.6, they obtain that the positive equilibrium is unstable and there exists a stable periodic solution. Next, we choose τ = 3.6 and b = 0.5, it can find that E * is unstable and a stable periodic solution exists (see Fig. 6). Furthermore, it increases b to 0.8, the numerical simulation shows that E * becomes stable (see Fig. 7). In addition, in [11], when τ = 2, E * is stable. Here by numerical simulations, it finds that the stability of E * can not vary despite the change in b (see Fig. 8). These may explain that b can change the oscillation property to stability. That is, b contributes to stabilize the system.
In the following, it investigates the effect of b for the single delay system under the parameters (39). Let  Fig. 9A), but b = 0.8 can obtain that E * is stable (see Fig. 9B).
Furthermore, we choose an another parameters b = 0.5, l = 7 and the other parameters as (39). When τ 2 = 0, from Theorem 8, it can obtain the (Q-τ 0 1 ) plane plots for different k (see Fig. 10A). Under the curve k = 0, E * is stable, otherwise, unstable. Using the above parameters, it has Q = 0.1531. Choosing τ 1 = 1 and τ 1 = 3, respectively, it can obtain the stability of E * (see Figs. [10][11]. Furthermore, by choosing the different k, the system occurs the Hopf-Hopf bifurcation since there are two pairs of imaginary roots simultaneously (see Fig. 11B). This leads to more complicated dynamical behavior, for example, several cycles may coexist near the intersection point. Choosing the different τ 2 , it can get the different (Q, τ 1 ) plane plots where ensures the increasing (decreasing) number of the positive real part roots for the fixed k = 0 when the parameters cross these critical boundaries (see Fig. 12). It can choose "+" (increasing) or "-" (decreasing) depending on the transversal condition. From Fig. 12B, at P 0 the Turing-Hopf bifurcation occurs.
At last, we investigate the effects of two different delays. Choosing r = 0.5, μ= 1.0, α = 0.10, γ = 2, b = 0.5, l = 7. (40) It can obtain that τ 0 1k is the increasing function of k (see Fig. 13A). Hence, we only need compute τ 0 10 . Under the above parameters, w + = 0.3005 and τ 0 10 = 1.8398. Fixed τ 1 = 1.7, when τ 2 = 0, E * is stable for any k ∈ N. We can draw the (ω, F 0 (ω)) plot shown in Fig. 13B, which can see that there are two roots defined by ω 1 = 0.3129 and ω 2 = 0.3847. However, when k ≥ 1, F k (ω) = 0 has no any positive roots. That is, for τ 1 = 1.7, only two purely imaginary roots are feasible when k = 0, the other situations have no any purely imaginary roots. Furthermore, we get two series of τ 2 as follows τ  Choosing τ 2 = 5, 15, 22, 75, respectively, we can obtain the phase plots as Fig. 14. The stability switch phenomenon occurs. In addition, in [11], if two same delays were considered, the system can not occur the stability switches. This explains the importance of two different delays. These numerical simulations confirm the correctness of the theoretical analyses.

Conclusions
In this paper, we analyze in detail a mussel-algae system (4) which improves that in [11] by introducing the half-saturation constant b and two different delays τ 1 and τ 2 . The parameters can produce some complicate A B Fig. 13 A (k, τ 0 1k ) plot. B (ω,F 0 (ω)) plot effects, such as Turing instability, steady-state bifurcation, Hopf bifurcation, Hopf-Hopf bifurcation or Hopfsteady state bifurcation. Some main results are given as follows: -Under certain conditions (A1, A2 or A3), there may be one or two positive equilibrium points in system (4). The set of parameters for the existence of positive equilibrium is given. -For system (3), • The existence and uniqueness of solutions are proved by upper and lower solution method, under certain condition (b+1 > γ ), the optimal bounds of solutions are derived. • The stability of the positive equilibrium depends on the half-saturation parameter b without considering spatial effects. The parameter r does not affect the existence of the equilibrium, but it will affect the stability of the equilibrium and lead to the occurrence of Hopf bifurcation. • Turing instability can occur under certain conditions, the numerical simulations find that the range of Turing instability becomes larger as b increasing.
-For system (4), • When the b + 1 > γ , the boundary equilibrium E 0 is globally asymptotically stable. • By analyzing the two-parameter characteristic equation (20), it is found that Eq. (20) may have single or double zero roots under certain conditions. When the parameter τ 2 is not considered, the condition and critical value of τ 1 such that the roots of Eq. (20) have strictly negative real parts are obtained. • The change of positive real roots in the plane (Q, τ 1 ) is obtained. When (Q, τ 1 ) is fixed in the stable region, considering the influence of τ 2 , the conditions for E * absolute stability and the conditions for Hopf bifurcation occurrence are obtained. Through numerical simulation, it can find that τ 2 causes the system to produce stable switching phenomenon. • The possible intersection of steady-state branch line and Hopf branch line in (Q, τ 1 ) plane is discussed, codimension 2 bifurcation may occur.
Our results suggest that the positive equilibrium will lose its stability when τ 1 or τ 2 pass through some critical values, and there will be periodic oscillations in populations of species. Compared to the system in [11], by numerical simulations, it can be found that a large halfsaturation constant b can contribute to stabilize the system whether the system is single time delayed or multiple time delayed. Under certain conditions, the double delayed system can occur the stability switches phenomena, which is impossible when τ 1 = τ 2 . This shows that two single delays can make system more complicated. Beyond that, by numerical simulations, the system may occur more complicated dynamic behavior, such as codimension-2, it still needs furthermore to confirm from the theory. In addition, in this paper, it has investigated the effects of two different delays, but it does not give the stability changes of equilibrium in (τ 1 , τ 2 ) parameter plane, which needs the further consideration.