Bogdanov-Takens bifurcation in a predator prey model with age structure

The results obtained in this article aim at analyzing Bogdanov-Takens bifurcation in a predator prey model with an age structure for the predator. Firstly, we give the existence result of the Bogdanov-Takens singularity. Then we describe the bifurcation behavior of the parameterized predator prey model with Bogdanov-Takens singularity.


Introduction
The goal of this article is to study an example of Bogdanov-Takens normal form computation for a non standard class of partial differential equation. Bifurcation theory has a long history and for infinite dimensional dynamical systems an important reduction tool is the center manifold theory. Indeed, the existence of such a locally invariant manifold permits to identify an infinite dimensional system to a finite dimension system of ordinary differential equations. Nevertheless even if the reduction is theoretically possible its applicability is not always guaranteed because such an invariant manifold is only implicitly defined.
Normal form theory has been extended to various classes of partial differential equations. In the context of partial differential equations many classes of equations have been considered and we refer to [2,6,9,10,11,12,13,14,15,24,25,26,31,35,43,44] for a nice overview on this topic.
The center manifold that we are using here is extending the one obtained by Vanderbauwhede and Iooss [38]. For partial differential equations with homogeneous boundary condition we also refer to the book of Haragus and Iooss [17] for a very nice overview on this topic. Bogdanov-Takens normal form computation is also considered in [17, p.102] which corresponds to 0 2 -bifurcation. In the book of Magal and Ruan [30] the non densely defined case is considered. Roughly speaking the non-densely defined Cauchy problems correspond to the non-linear boundary conditions in some classes of partial differential equations.
Due to the nonlinear boundary condition for such a class of partial differential equation, such a system can only be rewritten as a non densely defined semilinear Cauchy problem. That is to say that v (t) = Lv(t) + F (v(t)), for t ≥ 0, and v(0) = x ∈ D(L), (1.1) where L = D(L) ⊂ X → X is a linear operator on a Banach space X, and F : D(L) → X is a sufficiently smooth map that is also Lipschitz on bounded sets with F (0 X ) = 0 X and DF (0 X ) = 0 L(X) .
The above system is called non-densely defined whenever D(L) = X. It is important to observe that in such a problem the solution belongs to a strictly smaller subspace D(L), namely v(t) ∈ D(L), ∀t ∈ [0, τ ].
Therefore the equation (1.1) generates a semiflow on the smaller subspace D(L). The linearized equation around the equilibrium 0 is w (t) = Lw(t) + DF (0)(w(t)), ∀t ≥ 0, and w(0) = x ∈ D(L). (1. 2) The major difficulty in the normal form computations is coming from the fact that the linearized equation generates a C 0 -semigroup on D(L) while the range of F is the larger space X = D(L). Therefore the normal form theory used in this article is non standard compared with the one by Haragus and Iooss [17]. Our motivation is coming from a predator prey system with an age structure for the predator. The system considered is the following for t ≥ 0 and with initial distribution is the number of prey at time t, a → v(t, a) is the population density of predator structured with respect to the chronological age a of the predator. That is to say that a2 a1 v(t, a)da is the number of predator with age between a 1 and a 2 . We assume that the predators become adults (i.e. can reproduce) when they reach the age τ > 0 and we set which is the number of adult predators.
Here we use the simplified Monod-Haldane (also called simplified Holling type IV) functional response, that was used by Andrew [1]. This functional response takes the following form It means that the rate of predation increases when the number of prey is small and then decreases when the number of prey is above some threshold value. This type of functional response has also been used by Xiao and Ruan [40,41], Pan and Wang [33] and many others.
By using (1.3), we obtain the following delay differential equations When τ = 0, the number of adult predators A(t) coincides with the total number of predators a)da, ∀t ≥ 0 and the system (1.3) is equivalent to the following ordinary differential equations (1.6) Global qualitative analysis and Bogdanov-Takens bifurcation were proved for (1.6) by Ruan and Xiao [32]. The goal of this article is to extend their results to the age structured model (1.3)-(1.4).
Consider the new time variable t = t τ and the new age variable a = a τ and set U t = U τ t and v t, a = τ v τ t, τ a . Then By using this rescaling (and dropping the hat notation) we obtain the following system for t ≥ 0 with In the rest of the paper we will consider the rescaled system (1.
Recently, codimension 2 bifurcations, such as Bogdanov-Takens bifurcation [18,4], zero-Hopf bifurcation [42,3] and others, have attracted the attention of numerous researchers in the context of population dynamics. We also refer to Haragus and Iooss [17] for more results going in that direction.
Bogdanov and Takens's bifurcation firstly revealed the bifurcation phenomena (Hopf and homoclinic bifurcations) caused by Bogdanov-Takens singularity in planar systems which is an important discovery in the bifurcation theory. Bogdanov-Takens bifurcation is a well-studied example of a bifurcation in twoparameter family of autonomous differential equations, meaning that two parameters must be varied for the bifurcation to occur, the two equilibria (a saddle and a nonsaddle) of the system collide via a saddle-node bifurcation and a limit cycle is generated by an Hopf bifurcation around the nonsaddle equilibrium which degenerates into an orbit homoclinic to the saddle and disappears via a homoclinic bifurcation. We refer to the books [5,16,23,34]. We also refer to [19,20] for some nice illustration in the context of population dynamics.
Age structured models have a long history and we refer to the books of [7,21,22,39]. Here we use the integrated semigroup to define a notion of solution. This was first used in this context by Thieme [36].
The first result of Bogdanov-Takens bifurcation for an age structured model was investigated in [29] by using the integrated semigroup theory and the normal form theory developped in [28] for non densely defined Cauchy problems. We refer to the book of Magal and Ruan [30] for more results on this subject. In [29], in order to derive Bogdanov-Takens bifurcation we were forced to extend the original age structured model to a larger class of partial differential equation.
Here we do not need to extend the system. Therefore all the properties obtained by using Bogdanov-Takens bifurcation for the planar system are inherited by the system (1.3)-(1.4).
The plan of the paper is the following. In section 2, we study the equilibrium solutions. We will rewrite the system (1.7) as an abstract non densely defined Cauchy problem in section 3. The linearized equation and the characteristic equation are considered in section 4 and section 5, respectively. In section 6, we study Bogdanov-Takens singularity. Section 7 is devoted to Bogdanov-Takens bifurcation. Finally, in the last section, we present some numerical simulations.

Equilibria
By using system (1.7) we obtain the following system for the equilibria From the second equation of (2.1), we have v(a) = e −τ Da v(0), ∀a ≥ 0, and by using (2.2) and the last equation (2.1), we obtain respectively. We can see that there are two cases. Case 1: Assume that v(0) = 0. Then Case 2: Assume that v(0) > 0. Then from equation (2.3) we deduce that Remark 2.1 As a consequence of the above formula we have So we should consider the zeros of the following second order function Therefore whenever p(x) has some real roots they must be all strictly positive. Moreover the discriminant of p(x) is Then we have two cases. Case 2)-(a): and in the case we have and we have two solutions Whenever U > 0 and by considering the first equation of (2.1), we deduce that To summarize, we have the following result.

Lemma 2.2 (Equilibria) The boundary equilibria
always exist. Moreover we have the following alternatives: then the interior equilibrium is unique and given by with v(0) = τ De τ D A and A is given by (2.6).
then we have the following alternatives: the system has exactly one positive interior equilibrium and this interior equilibrium is given by the system has exactly two positive interior equilibria and the interior equilibria are given by

Abstract Cauchy problem formulation
In this section we will rewrite the system (1.7) as an abstract non densely defined Cauchy problem. We adapt the idea of Thieme in [36]. Consider X = R 2 × L 1 (0, ∞) endowed with the usual product norm Consider the closed subspace of X (3.1) and Then we can observe that X 0 = D(L) and L is non densely defined. Let F : X 0 → X be the nonlinear operator defined by (3.5) The partial differential equation (1.7) can be rewritten as the following nondensely defined abstract Cauchy problem

Computation of the linearized equation
In the rest of the paper, we will make the following assumption.
where U and v are described in Lemma 2.2-(i). Then we have w ∈ D(L) and Lw + F (w) = 0.
In order to work around a 0-equilibrium we use the following change of variable and in the rest of the paper the solution x(t) we will consider is a mild solution of the abstract Cauchy problem The linearized equation of (4.2) around the equilibrium 0 is

Computation of the characteristic equation
Abstract Cauchy problem (4.2) can be rewritten as is a linear operator and satisfies H(0) = 0 and DH(0) = 0. Next we will study the spectral properties of the linearized equation of (5.1) in order to investigate the dynamical behavior for system (5.1). Denote By applying the results of Liu, Magal and Ruan [27], we obtain that if λ ∈ Ω, then λ ∈ ρ(L) and Moreover, L is a Hille-Yosida operator. Now it remains to precise the spectral properties of P : Then we obtain v = v and we should have As a consequence we obtain the following lemma.
(ii) If λ ∈ ρ(P ) ∩ Ω, we have the following formula for the resolvent where δ and γ are given by the following formula and ∆(λ) are defined in (5.3).
From the explicit formula for the resolvent we deduce the following lemma. Since L is a Hille-Yosida operator, and DF (w) is bounded, P is also a Hille-Yosida operator. Consequently P 0 (the part of P in X 0 ) generates a strongly continuous semigroup {T P0 (t)} on X 0 . In order to apply the center manifold theorem and normal form theory, we need to study the essential growth bound of P 0 . By using the perturbation result in Thieme [37] or Ducrot Liu and Magal [8] we obtain the following estimation.

Proposition 5.3
The essential growth rate of the strongly continuous semigroup generated by P 0 is strictly negative.
Since the matrix ∆(λ) in (5.3) is triangular we obtain Assumption 6.1 We assume that Remark 6.2 The condition (6.1) implies that µ 2D e −τ D < κ. Hence under Assumption 6.1, the system (1.7) has a unique interior equilibrium.
Now we compute the projectors on the generalized eigenspace associated to eigenvalue 0 of P. From the above discussion we already knew that 0 is a pole of (λI − P ) −1 of finite order 2. This means that 0 is isolated in σ (P ) ∩ Ω and the Laurent's expansion of the resolvent around 0 takes the following form The bounded linear operator B P −1,0 is the projector on the generalized eigenspace of P associated to 0. We remark that So we have the following approximation formula and by computation, we can obtain the following lemma. of order 2, and the projector on the generalized eigenspace of P associated to the eigenvalue 0 is given by From the above results, we obtain a state space decomposition with respect to the spectral properties of the linear operator P . More precisely, the projector on the linear center manifold is defined by we obtain the basis {χ 1 , χ 2 } of X c defined by The matrix of P c in the basis {χ 1 , χ 2 } of X c is given by From the above analysis, we obtain the following theorem.
Theorem 6.4 Suppose that Assumption 6.1 holds. Then the unique positive equilibrium U , v of (1.7) is a cusp of codimension 2, i.e., system (1.7) has a Bogdanov-Takens singularity.

Bogdanov-Takens bifurcation
It is more interesting to determine a versal unfolding for the original system (1.3) with a Bogdanov-Takens singularity, i.e., to determine which of the parameters can be chosen as bifurcation parameters such that system (1.3) exhibits Bogdanov-Takens bifurcation. For this, choose κ and D in system (1.3) as the bifurcation parameters, i.e., consider 1/κ + β 1 and D + β 2 , where β := (β 1 , β 2 ) vary in a sufficiently small neighborhood of (0, 0). Adding these perturbations to system (1.3), using the same procedure in section 1 and 3, we obtain which is corresponding to where L is defined by In the following, we assume Assumption 6.1 holds. By making the change of variables x(t) := w(t) − w, (7.1) becomes where H(β, x(t)) = F (β, x(t) + w) − F (0, w). To determine that system (7.1) is the versal unfolding of system (1.3) with Bogdanov-Takens singularity, we include the parameter β into the state variable. Consider the system Under Assumption 6.1, 0 is an equilibrium of system (7.3). Note that for w =  In order to rewrite (7.3) as an abstract Cauchy problem, we set X = R 2 × X and consider the linear operator A : D (A) ⊂ X → X defined by Since P is a Hille-Yosida operator, we can prove that A is also a Hille-Yosida operator. Consider F : D (A) → X the nonlinear map defined by Then we have Now we can reformulate system (7.3) as the following system Note that σ (P ) ∩ iR = {0} and ω 0,ess ((P ) 0 ) < 0.
For the proof of the following Lemma the reader is referred to Liu, Magal and Xiao [29, lemma 5.1].
Lemma 7.1 Let Assumption 6.1 be satisfied. Then and for each λ ∈ ρ (A) , Next we compute the projectors on the generalized eigenspace associated to eigenvalue 0 of A. Note that Then λ n B P n,0 , and lim λ→0 1 2 We first have the following lemma. The projector on the generalized eigenspace of A associated to 0 is given in the following lemma. Lemma 7.3 0 is a pole of order 3 of the resolvent of A and the projector on the generalized eigenspace of A associated to 0 is given by

Lemma 7.2 We have
In particular we have  From the above results, we obtain a state space decomposition with respect to the spectral properties of the linear operator A. More precisely, the projector on the linear center manifold is defined by Let G ∈ V 2 (X c , D(A) ∩ X h ) be a homogeneous polynomial of degree 2. Consider the following global change of variable . By using the results in Liu, Magal and Ruan [28, lemma 5.1], we obtain that the Taylor expansion of up to second order term on the local center manifold is given by the ordinary differential equation on X c : , ς c (t)) expressed in terms of the basis {e 1 , e 2 , e 3 , e 4 }. ς c can be expressed by By computation, we have Thus by choosing χ = 2 we obtain the Taylor expansion of up to second order term on the local center manifold as follows Following the procedure of deriving normal form in [5], we can transform the above system into where Theorem 7.4 Assume that Assumption 6.1 holds. Then system (7.2) can undergo Bogdanov-Takens bifurcation in a small neighborhood of the unique positive equilibrium as the bifurcating parameters (β 1 , β 2 ) vary in a small neighborhood of (0, 0). More precisely, there exist four bifurcation curves: 2 saddle-node bifurcation curves SN + and SN − , Hopf bifurcation curve H and homoclinic bifurcation curve HL, in the small neighborhood of (0, 0) of parameter plane (β 1 , β 2 ), such that system (7.2) has a unique stable limit cycle as (β 1 , β 2 ) lies between H and HL, and no limit cycle for system (7.2) outside this region. The local representations of these bifurcations curves are given by , and α 2 < 0 .
The corresponding bifurcation diagram is shown in Figure 1.

Numerical simulations
In order to illustrate the Bogdanov-Takens bifurcation for system (7.2), we present some numerical simulations.

Symbol
First  Table 1: The formulas used for κ and µ in the first column of the table serve to satisfy Assumption 6.1. We fix the parameter τ = 0.1 (respectively τ = 1) in the second column (respectively third column). The parameters in the second column are computed by using the formula for H. The parameters in the last line are computed by using the formula for HL.
The parameters values used for the simulation of system (7.2) are listed in Table 1. In Figure 2 and Figure 3 we run a simulation of the system (7.2) by using respectively the second and third column of Table 1. In the simulations we use the equilibrium solution of the system whenever β 1 = β 2 = 0. Namely we use the formula in Lemma 2.2-(i). Figure 2 corresponds to a simulation of the model (7.2) with the parameters fixed in the second column of Table 1. We observe that the solution is converging to a periodic orbit. Actually since τ = 0, A(t) coincides with the total number of predators V (t) and the Figure 2-(c) corresponds the convergence of the solution to a periodic orbit of system (7.2) in the phase plane (U, V ). Figure 3 corresponds to a simulation of the model (7.2) with the parameters fixed in the third column of Table 1. In this figure we fix in particular τ = 1. By comparing Figure 2-(c) and Figure 3-(c) we observe that the limit behavior is much more irregular when τ = 1 than when τ = 0.1. In Figure 3-(c) the limit cycle follows the U -axis and the A-axis. The picks in Figure 3-(a) and    Table 1. In figures (a) and (b) we plot the function t → U (t) and t → A(t) respectively. In figure (c) we plot the function t → (U (t), A(t)) in the phase plane (U, A). In figure (d) we plot the (t, a) → v(t, a) the age structured density of predators when the time t varies.