Bifurcation and temporal periodic patterns in a plant–pollinator model with diffusion and time delay effects

ABSTRACT This paper deals with a plant–pollinator model with diffusion and time delay effects. By considering the distribution of eigenvalues of the corresponding linearized equation, we first study stability of the positive constant steady-state and existence of spatially homogeneous and spatially inhomogeneous periodic solutions are investigated. We then derive an explicit formula for determining the direction and stability of the Hopf bifurcation by applying the normal form theory and the centre manifold reduction for partial functional differential equations. Finally, we present an example and numerical simulations to illustrate the obtained theoretical results.


Introduction
It is believed that the explosive diversification and present-day abundance of flowering plants is due to their co-evolution with animal pollinators, especially insects [13]. The interactions between flowering plants and their insect pollinators remain an important ecological relationship crucial to the maintenance of both natural and agricultural ecosystems [15]. Mathematical modeling plays a useful role in pollination research and various mathematical models have been proposed to study plant-pollinator population dynamics, see Soberon and Del Rio [24], Lundberg and Ingvarsson [19], Jang [14], Neuhauser and Fargione [20], Fishman and Hadany [8], Wang et al. [29], Wang [26], and the references cited therein.
Consumer-resource systems model some biological phenomena and relationships between consumer and resource in the real world. A resource is considered to be a biotic population that helps to maintain the population growth of its consumers, whereas a consumer exploits a resource and then reduces its growth rate. Consumer-resource systems have been extensively studied by many researchers (see Chamberlain and Holland [3], Holland and DeAngelis [11], Li et al. [17], Neuhauser and Fargione [20], Wang and DeAngelis [27], Wang et al. [28]). Bi-directional consumer-resource interactions occur when each species acts as both a consumer and a resource of the other. Uni-directional consumer-resource interactions occur when one acts as a consumer and the other as a material and/or energy resource, but neither acts as both.
Recently, Wang, DeAngelis and Holland [29] derived a plant-pollinator model based on unidirectional interactions between plants and pollinators [11]. Pollinators travel from their nest to a foraging patch, collecting food, flying back to their nests, and unloading food. Interacting with flowers individually, the pollinators remove nectar, contact pollen, and provide pollination service. Therefore, the plants provide food, seeds, nectar, and other resources for the pollinators, while the pollinators have both positive and negative effects on the plants. Let N 1 and N 2 represent the population densities of plants and pollinators, respectively. The plant-pollinator model takes the following form: (1) where a, b, r 1 , β 1 , d 1 , d 2 , α 12 , and α 21 are positive constants. The parameter r 1 is the intrinsic growth rate of the plants and d 1 the self-incompatible degree. Following Fishman and Hadany [8], the positive effect of pollinators on plants is described by the Beddington-DeAngelis functional response aN 1 N 2 /(1 + aN 1 + bN 2 ), where the parameter a is the effective equilibrium constant for (undepleted) plant-pollinator interaction, which combines traveling and unloading times spent in central place pollinator foraging, with individual-level plant-pollinator interaction. b denotes the intensity of exploitation competition among pollinators (Pianka [21]). Since a is fixed, the parameter α 12 is regarded as the plants efficiency in translating plant-pollinator interactions into fitness (Beddington [2], DeAngelis et al. [6]) and α 21 is the corresponding value for the pollinators. β 1 denotes the per-capita negative effect of pollinators on plants. d 2 is the per-capita mortality rate of pollinators. Wang et al. [29] studied the globally asymptotically stability of the positive equilibria and demonstrated mechanisms by which interaction outcomes of this system vary with different conditions. In particular, it was shown in [29] that system (1) has no periodic orbits or cycle chains in the positive quadrant. In order to reflect the dynamical behaviours of models depending on the history, it is necessary to incorporate time delay into the models. Following Adams et al. [1], we assume that there is a time delay τ > 0 in the process when the pollinators translate plant-pollinator interactions into the fitness. Also, as pollinators travel between their nests and foraging patches, we further introduce the spatial diffusion with zero-flux boundary conditions. Thus, the plant-pollinator model with diffusion and time delay effects is described by the following delayed reaction-diffusion system: x ∈ , t > 0, where D 2 > 0 denotes the diffusion coefficient of pollinators. is a bounded open domain in R n (n ≥ 1) and its boundary ∂ is smooth, = ∂ 2 /∂x 2 1 + · · · + ∂ 2 /∂x 2 n is the Laplacian operator in R n , ν is the outer normal direction on ∂ , and the homogeneous Neumann boundary conditions reflect the situation where the population cannot move across the boundary of the domain.
Throughout this paper, without of loss of generality, we consider the domain = (0, π). Thus, = ∂ 2 /∂x 2 . We also assume that (φ, ψ) ∈ C := C([−τ , 0], X) and X is a suitable Hilbert space. For example, we can take with the inner product ·, · . The rest of the paper is organized as follows. In Section 2, we consider the corresponding characteristic equation of system (2) and give conditions on the stability of the positive constant steady-state and the existence of Hopf bifurcation. In Section 3, by applying the normal form theory and centre manifold reduction of partial functional differential equations (PFDEs) (Wu [30], Faria [7]), an explicit algorithm for determining the direction and stability of the Hopf bifurcation is given. Finally, some numerical simulations are included to support our theoretical predictions in Section 4 and a brief discussion is given in Section 5.

Stability and Hopf bifurcation
In this section, we consider the local stability of the positive constant steady-state and the Hopf bifurcation of system (2) by regarding the time delay τ as the bifurcation parameter. We assume that 21 , We can prove that, if (A1) or (A2) hold, then system (2) has two boundary equilibria E 0 (0, 0), E 1 (r 1 /d 1 , 0), and a unique positive constant steady-state E * (N * 1 , N * 2 ), where Then system (2) can be rewritten as The positive equilibrium E * (N * 1 , N * 2 ) of system (2) is transformed into the zero equilibrium of system (3). Let By the definition of the above functions, for i, j, l ∈ N 0 = {0, 1, 2 . . .}, define f (1) ij (i + j ≥ 1) and f (2) ijl (i + j + l ≥ 1) as follow: in particularly Obviously, we have α 1 + γ 2 < 0. By Taylor expansion, Equation (3) becomes Let u 1 (t) = u(t, ·), u 2 (t) = v(t, ·) and U = (u 1 , u 2 ) T . Then system (4) can be rewritten as an abstract differential equation in the phase space C := C([−τ , 0], X), where The linearized system of system (5) at (0, 0) has the form: and its characteristic equation is where y ∈ dom( )\{0} and dom( ) ⊂ X. It is well known that the Laplacian operator on X has eigenvalues −k 2 , k = 0, 1, 2, . . . , with corresponding eigenfunctions form a basis of X. Thus, any y ∈ X can be expanded as Fourier series in the following form: Hence, we conclude that the characteristic equation (7) is equivalent to the following sequence of characteristic equations: (9) Notice that (8) with τ = 0 is the characteristic equation of the linearization of (2) without delay at the positive equilibrium. Because D 2 k 2 − α 1 − γ 2 > 0, so the characteristic equation (8) with τ = 0 does not have a pair of purely imaginary roots for any k ∈ N 0 with N 0 := {0, 1, 2, . . .}. According to the Hopf bifurcation theorem, we obtain the following result.
Let λ = iω(ω > 0) be a purely imaginary root of Equation (8) Separating the real and imaginary parts in the above equation, we obtain which imply that i.e.
Set z = ω 2 , (12) is transformed into If (13) has only one positive root which is denoted by z k . Hence Equation (12) has only one positive root w k + = √ z k . From Equation (10), we know that Equation (8) with k ∈ N 0 has a pair of purely imaginary roots ±iw k with Therefore, λ = iw k + is a simple root of (8) for k ∈ N 0 .
Then λ(τ ) satisfies the following transversality condition: Proof: Differentiating both sides of Equation (8) with respect to τ yields From Equation (10), we have By inserting the expression of (w k + ) 2 into the last expression, we obtain that sign Re dλ dτ The proof is complete.

Theorem 2.10: Assume that
The following statements are true: is a Hopf bifurcation value of system (2) and these Hopf bifurcations are all spatially inhomogeneous.

Properties of Hopf bifurcations
In this section, we shall study the direction, stability and the period of bifurcating periodic solution by applying the normal form theory and the centre manifold theorem presented in [7,10,30]. Let τ k j ∈ F ∪ {τ 0 j , j ∈ N 0 }. Normalizing the delay τ in system (4) by the timescaling t → t/τ , Equation (4) is transformed into Let τ = τ k j + α, α ∈ R, u 1 (t) = u(t, ·), u 2 (t) = v(t, ·), and U = (u 1 , u 2 ) T . Then system (16) can be rewritten in the abstract form in the space C := C([−1, 0], X) as where L(α)(·) : C → X and F : C × R → X are defined by , f (ϕ, α), Consider the linear equation According to results in Section 3, we know that the origin (0, 0) is an equilibrium of Equation (16), and under some conditions, the characteristic equation of (19) has a pair of simple purely imaginary eigenvalues k = {iw k + τ k j , −iw k + τ k j }. We now consider the ordinary functional differential equation: By the Riesz representation theorem, there exists a 2 × 2 matrix function η(θ, τ k j ), θ ∈ [−1, 0], whose entries are of bounded variation such that for φ ∈ C([−1, 0], R 2 ). In fact, we can choose Let A(τ k j ) denote the infinitesimal generator of the semigroup induced by the solutions of system (20) and A * be the formal adjoint of A(τ k j ) under the bilinear pairing for ψ ∈ C([0, 1], R 2 ), φ ∈ C([−1, 0], R 2 ). From the previous section, we know that A(τ k j ) has a pair of simple purely imaginary eigenvalues ±iw k + τ k j . Because A(τ k j ) and A * are a pair of adjoint operators (see Hale [9]), ±iw k + τ k j are also eigenvalues of A * . Let P and P * be the centre subspace, that is, the generalized eigenspace of A(τ k j ) and A * associated with k respectively. Then P * is the adjoint space of P and dim P = dim P * = 2. Direct computations yield the following results.

Lemma 3.1: Let
Then form a basis of P with k and form a basis of P * with k .
and construct a new basis for P * by = ( 1 , 2 ) T = ( * , ) −1 * . Then ( , ) = I 2 , where I 2 is the identity matrix. In addition, f k := (β 1 k , β 2 k ), where Then the centre subspace of linear equation (19) is given by P CN C, where and we can decompose C([−1, 0], X) as C = P CN C ⊕ P S C, in which P S C denotes the complement subspace of P CN C in C.
Let A τ k j be the infinitesimal generator induced by the linear system (19), and Equation (17) can be rewritten as the following abstract form: where By the decomposition of C, the solution of Equation (17) can be written as where and h(x 1 , x 2 , α) ∈ P S C, h(0, 0, 0) = 0, Dh(0, 0, 0) = 0. In particular, the solution of (17) on the centre manifold is given by Then we obtain Hence, Equation (28) can be transformed into where From Wu [30], z satisfies where Let g(z,z) = g 20 z 2 2 + g 11 zz + g 02z Notice that . Then by computation, we obtain the following quantities: g 02 = g 20 , cos kx, cos kx where with (2) 101 +ξ f (2) 011 ) + f (2) 200 +ξξf and Then system (2) becomes By computation, we have E * (N * 1 , N * 2 ) = (0.427839, 1.380211), w + 0 = 0.123963, τ 0 0 = 12.518011. First we choose τ = 10 < τ 0 0 and plot the solutions N 1 (t, x) and N 2 (t, x) by using the software Matlab in Figure 1. From the numerical simulations we can see that the solutions of system (38) with τ = 10 tend asymptotically to the positive equilibrium E * (N * 1 , N * 2 ) = (0.427839, 1.380211). Under the same initial values, now we choose τ = 20 > τ 0 0 and plot the graphs of N 1 (t, x) and N 2 (t, x) in Figure 2. From Figure 2, we see that there exists a family temporal periodic solutions, which implies that Hopf bifurcation occurs for system (38) at τ 0 0 .
Most of these models are described by ordinary differential equations. Since pollinators travel between their nests and foraging patches, we believe that reaction-diffusion equations are more suitable to model the interactions between the plants and pollinators. We also assumed that there is a time delay in the process when the pollinators translate plant-pollinator interactions into the fitness and considered a plant-pollinator model with diffusion and time delay effects. As far as we know, there are no results for system (2) with diffusion and time delay. Firstly, by considering the distribution of eigenvalues of the corresponding linearized equation, stability of the positive constant steady-state and existence of spatially homogeneous and spatially inhomogeneous periodic solutions were studied. Secondly, by applying the normal form theory and the centre manifold reduction for partial functional differential equations, an explicit formula for determining the direction and stability of the Hopf bifurcation was given. Finally, to explain the obtained results, numerical simulations were presented.
Our results showed that if α 21 > ad 2 and either (A1) a 1 < 0, a 2 1 − 4a 0 a 2 = 0 or (A2) 4a 0 a 2 < 0 holds, where 21 , then system (2) has a unique positive constant steady-state E * (N * 1 , N * 2 ), in which The first inequality α 21 > ad 2 ensures the existence of a 0 , a 1 , a 2 , and N * 1 . Recall that α 21 is regarded as the pollinators efficiency in translating plant-pollinator interactions into fitness, a is the effective constant for plant-pollinator interaction, and d 2 is the percapita mortality rate of pollinators. This inequality means that the efficiency in translating plant-pollinator interactions into fitness of the pollinators must be greater than their mortality rate; otherwise the pollinators even cannot survive.
The inequality a 1 < 0 in (A1) is equivalent to which indicates that the ratio of the efficiencies in translating plant-pollinator interactions into fitness of the plants and pollinators is greater than a certain value. In this case, an additional condition a 2 1 − 4a 0 a 2 = 0 is needed to ensure the existence of E * (N * 1 , N * 2 ). Under the assumption (A2), it requires that 4a 0 a 2 < 0. Note that now a 0 > 0, so the condition is equivalent to a 2 < 0, which, in turn, is equivalent to The last inequality means that the intrinsic growth rate r 1 of the plants must be large enough compared to the death rates of the plants and pollinators.
We were interested in not only the effect of diffusion but also the effect of delay [4,12,31]. We found that system (2) without delay cannot undergo Hopf bifurcations at the positive constant steady-state. But, under certain conditions, system (2) undergoes Hopf bifurcations at the positive constant steady-state under the effect of delay. Recall that then the positive equilibrium E * is locally asymptotically stable if the time delay is less than a critical value τ < τ 0 , unstable when τ > τ 0 , and a family of periodic solutions bifurcates from E * when τ passes through τ 0 via Hopf bifurcation. Moreover, the direction, stability and period of the bifurcating periodic solutions can be determined analytically. Notice that Wang et al. [29] showed that the ODE model (2) does not have periodic solutions and Wang et al. [25] proved that the unique positive steady-state solution of a reaction-diffusion plant-pollinator model is a global attractor. Our results thus indicate that the time delay causes bifurcations and induces temporal periodic patterns in the diffusive plant-pollinator model. Such properties have been observed in many delay differential equation models [5,16]. This is similar to the observation in our other work [18] that oscillations occur in age-structured resource-consumer (plant-pollinator) models. Wang et al. [29] and Wang [26] indeed investigated three species plant-pollinator-robber models. Since the movement of the nectar robbers plays an important role in their invasibility and coexistence of all species, it will be very interesting to study the population dynamics of the three species diffusive plant-pollinator-robber models. We leave this for future consideration.