Bautin bifurcation in delayed reaction-diffusion systems with application to the Segel-Jackson model

In this paper, we present an algorithm for deriving the normal forms of Bautin bifurcations in reaction-diffusion systems with time delays and Neumann boundary conditions. On the center manifold near a Bautin bifurcation, the first and second Lyapunov coefficients are calculated explicitly, which completely determine the dynamical behavior near the bifurcation point. As an example, the Segel-Jackson predator-prey model is studied. Near the Bautin bifurcation we find the existence of fold bifurcation of periodic orbits, as well as subcritical and supercritical Hopf bifurcations. Both theoretical and numerical results indicate that solutions with small (large) initial conditions converge to stable periodic orbits (diverge to infinity).


INTRODUCTION
In recent years, Hopf bifurcation in reaction-diffusion systems with time delays has been widely studied [1,2]. The results are used to interpret the oscillation behaviors in many subjects, such as, biology [3,4], chemistry [5], and epidemics [6], etc. The authors usually investigated the change of stability of equilibria and the existence of Hopf bifurcations. Furthermore, they calculated the normal forms near the Hopf bifurcations to determine the stability of bifurcating periodic solutions and direction of Hopf bifurcations. The framework of deriving normal forms is due to the center manifold reduction technique introduced in [4,[7][8][9].
As is known to all, supercritical Hopf bifurcations with negative first Lyapunov coefficient lead a system to stable periodic orbits at the side of bifurcation point that the equilibrium becomes unstable; whereas subcritical Hopf bifurcations with positive first Lyapunov coefficient induce unstable periodic orbits at the side of bifurcation point that the equilibrium becomes stable. General theory about Hopf bifurcation can be found in [10][11][12] for ordinary differential equations, in [7] for functional differential equations, and in [8,13] for partial differential equations or partial functional differential equations. If we introduce another parameter and with the change of that parameter, Hopf bifurcation points forms a curve in the two-parameter plane. Thus, at the intersection of the supercritical Hopf bifurcation and subcritical Hopf bifurcation curves, there will be a point with zero first Lyapunov coefficient that makes possible Bautin bifurcation happen which is of codimension two [10][11][12]14].
Near a Bautin bifurcation, the dynamical behavior of a system is quite complicated: the bifurcation point separates branches of subcritical and supercritical Hopf bifurcations in the parameter space. For parameter values near the singularity, the system has two limit cycles which collide and disappear via a fold bifurcation of periodic orbits.
So far, there are few papers on codimension two bifurcations in delay equation, most of which focus in retarded functional differential equations [15,16], and neutral functional differential equations [17,18]. The study on codimension two bifurcations of delayed reactiondiffusion equations just begins in recent few years, and most of which are about Turing-Hopf bifurcation [19], Hopf-zero bifurcation [20] or Double Hopf bifurcation [21]. Near these bifurcations the spatially non-homogeneous steady-state solutions have been found as well as periodic solutions and quasi-periodic solutions. Bautin bifurcations have been found and analyzed in ordinary differential equations [22][23][24][25] or in delay differential equations [26][27][28].
However, to our best knowledge, there are not any research papers found in literature about Bautin bifurcations in reaction-diffusion systems with time delays.
Motivated by such a consideration, in this paper we will give universal derivations of the first and second Lyapunov coefficients for Bautin bifurcations in reaction-diffusion equations with time delays, then the normal form of Bautin bifurcations is given and has two different kinds of bifurcation diagrams. The results are then applied to a Segel-Jackson prey-predator model [29,30] In such a system with time delay and homogeneous boundary conditions, we find that Bautin bifurcation appears in (k, τ ) plane. Both subcritical and supercritical Hopf bifurcations are found theoretically near the bifurcation point.
The rest part of the paper will be organized as follows. In section 2, the fifth-order normal form near a pair of imaginary roots in a general reaction-diffusion system of N equations with time delays are derived, from which detailed and explicit formulae to determine the first and second Lyapunov coefficients for Bautin bifurcation are given. In section 3, we apply the results to a Segel-Jackson model, in which we find a Bautin bifurcation, and both subcritical and supercritical Hopf bifurcations are theoretically analyzed. To support the results some simulations are performed to illustrate the dynamical behaviors near the Bautin bifurcation. Finally a conclusion section completes this paper.

TIME DELAYS
We consider the following general form of a reaction-diffusion system with delays, Here the state variable U maps [−1, +∞) to X, and the state space is Notice that here we use homogeneous Neumann boundary conditions, and a one dimensional space interval [0, lπ]. Using the general notation of delayed equations [8], . . , F N ) T and satisfies F (µ, 0) = 0 and D ϕ F (µ, 0) = 0, which means U t = 0 is a steady states of system (2).
The linearized equation of system (2) at U t = 0 is Since we are about to calculate the fifth order normal form at a Bautin bifurcation, we first assume that when µ ∈ V , with V representing some neighbourhood of the origin in R 2 , the characteristic equation of (3) has a pair of roots α(µ) ± iω(µ) for some n = n 0 where α(µ) and ω(µ) are both continuously differentiable, and α(0) = 0, ω(0) = ω 0 > 0. We further assume that all the other roots of (4) have negative real parts for µ ∈ V . Thus the center manifold will be locally attractive and the local dynamical behaviors of (3) will be governed by the normal form on the center manifold [9].
The solution operator of Eq.(2) generates a C 0 -semigroup with the infinitesimal generator A µ defined by The domain of A µ is chosen by We rewrite Eq.(2) as an abstract ordinary differential equation with the nonlinearity Also denote by β n = (b n , b n , . . . , b n ) T , then for any φ = (φ (1) , φ (2) , . . . , φ (N ) ) T ∈ C , we can Now we have the action of A µ on φ n b n , denoted by A µ,n as with the function of bounded variation η n satisfying Consider the adjoint operator of A 0 , A * given by for some ϕ ∈ C * = C ([0, 1], X) and ϕ n =< ϕ, b n >. According to the formal adjoint theory given in [8,31], we define the bilinear form by In fact, Thus, we only need to define (ϕ n , φ n ) c from [7,8], as Let q(θ)b n 0 and q * (s)b n 0 be the eigenfunctions of A 0 and A * corresponding to iω 0 and −iω 0 respectively, with (q * , q) c = 1, (q * ,q) c = 0.
Write the centerspace P = {zqb n 0 +zqb n 0 |z ∈ C} and its orthogonal complement space then C = P Q and the state variable U t of Eq.(6) could be decomposed by where W (t, ·) ∈ Q.
Obviously z(t) = (q * b n 0 , U t ) and W (t, θ) = U t (θ) − 2Re{z(t)q(θ)b n 0 }. Thus, we havė where the nonlinear term is given by Using the center manifold theorem [8,9], W (t, ·) can be expressed by a series of z andz, starting at quadratic terms where z andz are approximately the local coordinates for center manifold C 0 in two directions qb n 0 andqb n 0 , respectively. For solution U t ∈ C 0 , we have U t = z(t)q(·)b n 0 +z(t)q(·) + W (z(t),z(t), θ) and denote by We write the Taylor formulaF Now the equation (9) on the center manifold is expressed bẏ In order to calculate the normal form, we expand According to the theory about normal form at a pair of imaginary roots [12]. We have the normal form till fifth order given bẏ and Notice that some g ij 's may depend on the form of W ij . By using (6) and (8), we havė Again, expending H(z,z, θ) as we have Recalling the definition of A 0 in (5), we know (18) is actually a group of ordinary differential equations about W ij (θ) with some given boundary conditions, solving which the explicit expressions of W ij and g ij are finally obtained.
Following the results given in [10], we know that near the Bautin bifurcation, when l 2 (0) < 0 two limit cycles (region III) collide on the curve standing for fold bifurcation of periodic orbits, then disappear, and leave the system a stable equilibrium (region I), which undergoes a supercritical Hopf bifurcation and lead the system to stable periodic oscillations (region II). This is shown in Figure 1 (a). When l 2 (0) > 0, the situation is shown in Figure   1 (b), where the fold bifurcation of periodic orbits eliminates two limit cycles (region III), and leave the system an unstable equilibrium (region I), which undergoes a subcritical Hopf bifurcation and unstable periodic orbits appear (region II).
In the coming section, an example of Segel-Jackson model will be analyzed, where we find a Bautin bifurcation with l 2 (0) > 0 appears.

EXAMPLE: SEGEL-JACKSON MODEL
In this section, we shall investigate the Bautin bifurcation for the Segel-Jackson model.
u(x, t) stands for the density of the prey and v(x, t) the predator. d is the diffusion rate of the prey whereas the diffusion rate of predator is normalized by 1. 1 + ku(x, t) stands for the reproduction rate per capita of prey. a is the ability of the predator hunting the prey.
Equipping (19) with homogeneous Neumann boundary condition ∂u ∂x x=0,lπ = ∂v ∂x x=0,lπ = 0, there are two constant equilibria of system (19): one is the trivial equilibrium E 0 = (0, 0) and the other is the unique positive equilibrium The linear equation of system (19) Then we have the characteristic equation of system (20) is In the following, we will discuss the stability of E * , regarding time delay τ and k as bifurcation parameters.
When τ = 0, Eq.(21) becomes In fact, this is just the case investigated in [29], and we state the main results here: when holds true, then the positive equilibrium E * is locally asymptotically stable.
In the following part, we will investigate the existence of Hopf bifurcations destabilizing E * as τ increases, where we always assume (H 1 ) holds true.
When τ > 0, following the method given in [32], we plug λ = iω into Eq.(21) and separate the real and imaginary parts, then obtain This is a linear equation about unknowns cos ωτ and sin ωτ . Once ω and n are fixed, we can solve from (23) that cos ωτ :=C n (ω) and sin ωτ :=S n (ω).
Adding the square of both sides of Eq.(23), we have an algebraic equation about ω 2 , denoted by (24) has two positive roots, given by Else if M 2 < 0 holds true, Eq.(24) has only one positive root, denoted by Plugging ω + n or ω ± n into Eq.(23), at most two sequences of possible Hopf bifurcation values are obtained, which are denoted by We know when τ = τ ± n,j , n, j = 0, 1, 2, · · · , Eq.(21) has a pair of imaginary roots. By re- the tildes for convenience, we have system (19) becomes where u = u(x, t), v = v(x, t). For any (φ 1 , φ 2 ) ∈ C := C([−1, 0], X) and X := {(u, v) ∈ H 2 (0, lπ) × H 2 (0, lπ) : du dx = dv dx = 0 at x = 0, lπ}, two nonlinear functionals are defined by and We will calculate the first and second Lyapunov coefficients of the normal forms when τ = τ * ∈ {τ ± n,j , n, j = 0, 1, 2, · · · } and for any fixed k, while Eq.(21) has a pair of imaginary roots λ = ±iω ± n for n = n 0 . For simplification, we denote ω ± n 0 by ω * . Let τ = τ * + ǫ. We can rewrite system (26) in an abstract form in the phase space and and Thus, the linearized equation of system (26) at the equilibrium (0, 0) has the forṁ We have that the solution operator of system (26) forms a C 0 semigroup with the infinitesimal generator A ǫ , The domain of A ǫ is chosen by In order to use the results we obtained in the precious section, in the space C we rewrite Eq.(29) as the abstract formU where the nonlinear term is given by In the state space X, we know, cos nx l (0, 1) T , cos nx l (1, 0) T , n = 0, 1, 2, · · · , are eigenfunctions of ∆ with the no-flux boundary conditions. Moreover, they form a basis of X which are normalized by b n = cos( nx l ) cos nx l L 2 . For convenience we use β n = {b n (1, 0) T , b n (0, 1) T } to be the basis of X.
Using the notations introduced in Section 2, for any φ = (φ (1) , φ (2) ) T ∈ C , we denote as the coordinates of φ on β n in X.
Obviously, we have, from the linear Eq.(32) Similar as the previous section we can define the adjoint operator A * of A 0 as As we stated at the beginning of this section, we can calculate the eigenfunctions of A 0 and A * at the eigenvalue iω * τ * and −iω * τ * , which are, respectively, qb n 0 and q * b n 0 .
By direct calculations, we have q(θ) = (1, q 1 ) T e iω * τ * θ and q * (s) = M(q 2 , 1) T e iω * τ * s . Here such that (q * , q) c = 1. By the formal adjoint theory, we have (q * , q) c = (q * ,q) = 0. By using the standard theory of phase space decomposition, we have the eigenvalues Λ = {±iω * τ * } can be used to decompose C by C = P Q. The center space P is given by P = {zqb n 0 +zqb n 0 |z ∈ C} and its orthogonal complement space is Thus, the state of system (34) could be decomposed, correspondingly, by where W (t, ·) ∈ Q.
Notice the relation between q * and q, we calculate z(t) by Thus, W (t, θ) = U t (θ) − 2Re{z(t)q(θ)b n 0 }. Then we obtain the reduced equation onto the center spaceż where the nonlinear term is given by with F 1 and F 2 already defined in (30) and (31).
Using the center manifold theorem [9], there exists a center manifold C 0 and on C 0 we have W (t, ·) can be expressed by series For solution U t ∈ C 0 , we have U t = z(t)q(·)b n 0 +z(t)q(·) + W (z(t),z(t), θ) and denote by We write the Taylor expressioñ In fact, from Eqs. (27)(28) and (30-31), we havẽ and all the rest higher-order derivatives are zero.

Now rewrite the equation on the center manifold (36) aṡ
In order to calculate the normal form, we expand g(z,z) = g 20 z 2 2 + g 11 zz + g 02z By comparing the same order of terms, we obtain the expressions of g ij 's, which, together with the detailed derivations, are given in the Appendix.
According to the results given in Section 2, the first and second Lyapunov coefficients can be expressed by and The transversality condition can be verified, numerically, by Thus the map (k, τ ) → ( α(k,τ ) ω(kτ ) , l 1 (k, τ )) is regular at (k * , τ * ), and the system undergoes a Bautin bifurcation at (k * , τ * ) with the bifurcation diagram topologically equivalent to Figure 1 (b). When k < k * = 0.3075, the Hopf bifurcation is subcritical, as shown in Figure   3 (a). The Bautin bifurcation diagrams is shown in Figure 2 (b). When k > k * , the Hopf bifurcation is supercritical. As shown in Figure 3 (b), the fold bifurcation of periodic orbits eliminates the stable and unstable periodic orbits.
In fact, we find that when k = 0.8 and τ ∈ (0.9566, 1.1550), solutions with small initial values converge to a stable periodic orbits (Figure 4), and those with large initial values diverge to infinity as time goes to infinity ( Figure 5).

CONCLUSION
In this paper, a universal and explicit method to calculate the first and second Lyapunov coefficients at a pair of imaginary eigenvalues in reaction-diffusion system with time delays is given, which can be used to determine the dynamics near a Bautin bifurcation point. As an example, the method is applied to the Segel-Jackson model. Near the Bautin bifurcation we theoretically proved that solutions with small (large) initial values are convergent to stable periodic oscillations (diverge to infinity).

ACKNOWLEDGMENTS
The authors wish to express their gratitude to the editors and the reviewers for the helpful comments. This research is supported by National Natural Science Foundation of China (11701120,11771109).

APPENDIX
In this Appendix, we will give detailed formulae to calculating g ij 's, then the first and second Lyapunov coefficients are determined.