Hopf bifurcations in a three-species food chain system with multiple delays

Abstract This paper is concerned with a three-species Lotka-Volterra food chain system with multiple delays. By linearizing the system at the positive equilibrium and analyzing the associated characteristic equation, the stability of the positive equilibrium and existence of Hopf bifurcations are investigated. Furthermore, the direction of bifurcations and the stability of bifurcating periodic solutions are determined by the normal form theory and the center manifold theorem for functional differential equations. Finally, some numerical simulations are carried out for illustrating the theoretical results.


Introduction
The study on the dynamics of predator-prey system is one of the dominant subjects in ecology and mathematical ecology due to its universal existence and importance. As to our knowledge, delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and cause the populations to fluctuate. Thus, time delays of one type or another have been incorporated into mathematical models of population dynamics due to maturation time, capturing time or other reasons. In the last decades, many authors have explored the dynamics of systems with time delay and many interesting results have been obtained [1][2][3][4][5][6][7][8][9][10][11][12][13].
However, there may be more species in a habitat and they can construct a food chain. Therefore, it is more realistic to consider a multiple-species predator-prey system. Recently, Baek and Lee [14] proposed the following three-species food chain model with the Lotka-Volterra functional response 8 < : dx.t/ dt e 1 > 0 measure the hunting of the mid-level predator to the lowest-level prey and the top-level predator to the mid-level predator; d 1 > 0 and d 2 > 0 denote the death rates of the mid-level and top-level predator, respectively. In system .1/, the gestation periods and maturation time of some species are all omitted. However, some species need to take time to have the ability to reproduce and capture food. Based on this fact, by incorporating time delays into the maturation time, Cui and Yan considered the following delayed differential system in [15] 8 < : dx.t / dt D x.t /.a bx.t / cy.t //; where 1 0 denotes the time from birth to having the ability to predate for top-level predator and 2 0 represents the maturation time that the mid-level predator can be served as food for the top-level predator, respectively. Delayed models similar to .2/ have been investigated widely by many authors and lots of interesting results have been obtained. For more details see [16][17][18].
In system .2/, the author just considered the time for the top-level predator to have the ability to predate and the maturation time for the mid-level predator to be served as food for the top-level predator. But the time for the mid-level predator to have the ability to predate and the maturation time for the lowest-level prey to be served as food for the mid-level predator are omitted. Based on this fact, by incorporating delays into the above terms, we consider a more complicated system with multiple delay 8 < : where 1 0 and 2 0 denote the time from birth to having the ability to predate for the mid-level predator and the top-level predator. For a special case, we assume that 2 C 3 and 1 C 3 represent the maturation time for the lowest-level prey to be served as food for the mid-level predator and the maturation time for the mid-level predator to be served as food for the top-level predator, respectively. We address the question how the time delays that we incorporated affect the dynamical properties of the system .3/. So the aim of this paper is to study the dynamical behaviors of the system .3/, for which we investigate the stability and Hopf bifurcation of a three-species food chain system with multiple delays. We would like to mention that the bifurcation in a predator-prey system with a single or multiple delays had been investigated by many researchers [19][20][21][22]. However, to the best of our knowledge, few results for system .3/ have been obtained. Therefore, the research of this case is worth considering.

Stability of positive equilibrium and existence of local Hopf bifurcations
For convenience, we introduce new variables / and assume that D 1 C 2 C 3 , so that system .3/ can be written as the following system with a single delay: It is easy to see that the system .3/ has a unique positive equilibrium E W .x 0 ; y 0 ; z 0 / provided that the condition t / z 0 , and use relations a bx 0 cy 0 D 0, d 1 C c 1 x 0 e 1 z 0 D 0, d 2 C e 2 y 0 D 0, then system .4/ can be rewritten as the following equivalent system To study the stability of the equilibrium E , it is sufficient to study the stability of the origin for system .5/. The linearized system of system .5/ at origin is The characteristic equation of system .6/ is Next, we will investigate the distribution of roots of Eq:.7/. Obviously, D 0 is not a root of Eq:.7/. When D 0, the characteristic equation becomes It can be seen that B > 0, B.A C D/ C > 0, and C > 0. Therefore, if follows from the Routh-Hurwitz criteria that all roots of .8/ have negative real parts and thus the zero equilibrium of system .5/ is asymptotically stable when D 0. Now, we examine when the characteristic equation has pairs of purely imaginary roots. For > 0, if i !.! > 0/ is a root of Eq:.7/, then ! should satisfy the following equations Thus si n.! / D Squaring and adding both the equations of .9/, we have .ACD/ 2 , then Eq:.10/ becomes Let m D , D 0 D n 2 4 C m 3 27 . Similar to discussion in [24], assume that D 0 > 0, then from the Cardano's formula for the third-degree algebra equation we know that the equation f .z/ D 0 had only one real root z 1 . If D 0 D 0, then the equation f .z/ D 0 has three real roots z 1 ; z 2 and z 3 (where z 2 D z 3 ), and in this case we define z 2 by max fz 1 ; z 2 g. If D 0 < 0, we know that the equation f .z/ D 0 has three different real roots denoted by s 1 ; s 2 and s 3 . In this case ,we define z 3 D max fs 1 ; s 2 ; s 3 g. According to Lemma 2.2 in [24], without loss of generality, we can suppose that Eq:.10/ has four positive real roots, denoted by z 1 ; z 2 ; z 3 ; z 4 , respectively. Then Eq:.9/ should also have four positive real roots ! 1 D p Then . j k ; ! k / are solutions of Eq:.7/ and D˙i ! k are a pair of purely imaginary roots of Eq:.7/ with D j k .
Then 0 is the first value of such that Eq:.7/ has purely imaginary roots. In the following discussions, for the sake of convenience, we denoted j k by j .j D 0; 1; 2; / for fixed k 2 f1; 2; 3; 4g. Let . / D˛. /˙i !. / be the root of Eq:.7/ near D j satisfying Proof. Differentiate the two sides of Eq:.7/ with respect to . For the sake of simplicity, we denote j and ! 0 by , !, respectively, then This completes the proof of Lemma 2:1.
Since the multiplicity of roots with positive real roots of Eq:.7/ can change only if a root appears on or crosses the imaginary axis as time delay varies, similarly to Lemma 2:4 and Theorem 2:5 in [24], we have the following results

Direction of Hopf bifurcations and stability of bifurcating periodic solutions
In the previous section, we have already obtained that, under certain conditions, the system .4/ can undergo Hopf bifurcation at the positive equilibrium E W .x 0 ; y 0 ; z 0 / when takes some critical values D j ; .j D 0; 1; 2; /. In this section, by employing the normal form theory and center manifold theorem introduced by Hassard et al. [23], we shall present the formula determining the direction of the Hopf bifurcation and the stability of bifurcating periodic solutions of .4/.
Let u i .t/ D Q u i . t /, D j C , 2 R, then D 0 is the Hopf bifurcation value for system .4/, dropping the bars for simplification of notations, system .4/ becomes the following functional differential equation in C D C.OE 1; 0; R 3 /, where u.t / D .u 1 .t /; u 2 .t /; u 3 .t // T 2 R 3 , and L W C ! R 3 ; f W R C ! R 3 are given by In fact, we can choose where ı is the Dirac delta function. Â D 0: Then system .12/ is equivalent to where u t .Â / D u.t C Â / for Â 2 OE 1; 0. For q. 1/ D q.0/e i ! 0 j , then we obtain Similarly, we can obtain the eigenvector q .s/ D D.1; q 1 ; q 2 /e i! 0 j s of A corresponding to i ! 0 j , where In order to assure hq .s/; q.Â /i D 1, we need to determine the value of D. By .18/, we have Therefore, we can choose D as Next we will compute the coordinate to describe the center manifold C 0 at D 0. Let u t be the solution of .17/ when D 0. Define z.t / D hq ; u t i; On the center manifold C 0 , we have z and N z are local coordinates for center manifold C 0 in the direction of q and N q . Note that W is real if u t is real. We only consider real solutions for solutions u t 2 C 0 of .17/. Since D 0, we have where Following from .19/, .20/, .14/ and using the algorithms given in [23], we can get the coefficients which will be used to determine the important quantities g 20 D 2 j N D.a 11 C a 12 q 1 C a 21 N q 1 q 1 e i! 0 j C a 22 N q 1 q 1 q 2 C a 31 N q 2 q 1 q 2 e i! 0 j /; where and Similarly to the algorithms given in [25], we can obtain Similarly, we get ; .k D 0; 1; 2; /: Theorem 3.1.

Numerical simulations
In this section, we give some numerical simulations supporting our theoretical predictions. As an example, we consider the following system