Transient spatio-temporal dynamics of a diffusive plant–herbivore system with Neumann boundary conditions

ABSTRACT In many existing predator–prey or plant–herbivore models, the numerical response is assumed to be proportional to the functional response. In this paper, without such an assumption, we consider a diffusive plant–herbivore system with Neumann boundary conditions. Besides stability of spatially homogeneous steady states, we also derive conditions for the occurrence of Hopf bifurcation and steady-state bifurcation and provide geometrical methods to locate the bifurcation values. We numerically explore the complex transient spatio-temporal behaviours induced by these bifurcations. A large variety of different types of transient behaviours including oscillations in one or both of space and time are observed.


Introduction
Plant-herbivore interactions have been extensively studied by many researchers. Generally speaking, plant-herbivore interactions fall into the category of predator-prey systems. A general introduction to the theory of population dynamics and predator-prey interactions can be found in the book by Turchin [25]. Mathematically, a predator-prey system consists of two coupled differential equations describing the dynamics of prey abundance, N, and predator abundance, P: Here, μ is the mortality rate of predators in the absence of prey, the function h is dynamics of plants in the absence of predator. The function f is called the predators' functional response [14], and the function g is referred as the numerical response function [22]. Typical functional response functions include the Holling type I, II, III [14], and the Beddington-DeAngelis type [2,6]. Rich dynamics are possible for predator-prey models with non-monotonic response functions (also called Holling type IV) [21,31], age-structure [5], delays [20,28] and diffusions [7].
Many existing predator-prey models assume that the numerical response is proportional to the functional response, i.e. g = f , where is the conversion coefficient denoting the energetic efficiency in converting consumption into reproduction. In general, the numerical response is not necessarily proportional to the functional response [19]. For instance at a very high level of consumption, additional consumption or just killing prey without eating them may not lead to proportionally more reproduction of predators [10,11]. This is particularly true if the chemical defense of plants is considered as the negative effect of plant toxin on herbivore can lead to a decrease in the growth rate [9]. Recently, a complex toxin-determined functional response function is proposed in [4,8,17] describing the reduction in the consumption of plants by the herbivore due to chemical defenses.
In this paper, we study the dynamics of plant-herbivore interactions with plant defense by assuming the numerical response function is not necessarily proportional to the functional response. More precisely, we consider the following diffusive plant-herbivore system: where u(x, t) and v(x, t) represent the population densities of plant and herbivore at location x ∈ and time t > 0, respectively. For simplicity, we take = (0, π). The two constants d 1 and d 2 are the diffusion coefficients; r is the intrinsic growth rate of the plant in the absence of herbivore; κ is the death rate of the herbivore; the positive constants β and σ are the scale parameters of functional response (characterized by the function u/(1 + βu)) and numerical response (characterized by the function u/(1 + σ u)), respectively. The model has been non-dimensionalized using a carrying capacity to scale the prey, the attack rate as a time scale, the conversion efficiency as a predator scale, and the domain size for a spatial scale. For (2) with nonzero Dirichlet boundary conditions, its transient spatio-temporal dynamics has been explored in our recent paper [27]. If β = σ , then (2) reduces to the diffusive predator-prey model in which numerical response is assumed to be proportional to functional response, which has been extensively studied in the literature. See, for example, [26,29,30] and the references therein. We point out that transient oscillatory patterns driven by the time delay were also recently investigated in [24] for a delayed diffusive non-local blowfly equation with delay under the zero-flux boundary condition. The rest of this paper is organized as follows. Section 2 is devoted to the study of the local system without diffusion. In Section 3, we establish stability results for steady states of System (2). Hopf bifurcations and steady-state bifurcations are investigated in Section 4. In Section 5 we numerically explore possible transient spatio-temporal behaviours induced by Hopf and steady-state bifurcations. We summarize and discuss our results in Section 6.
There exists a unique solution which remains non-negative and is ultimately bounded. Moreover, the set defined in Equation (4) is positively invariant.
Note that System (3) has two boundary equilibria: E 0 = (0, 0) and E 1 = (1, 0). It is straightforward to prove the following result concerning the stability of E 0 and E 1 .

Proposition 2.2.
The equilibrium E 0 is always a saddle point with a one-dimensional stable manifold and a one-dimensional unstable manifold. The equilibrium E 1 is locally asymptotically stable if and is unstable if System (3) has a positive equilibrium, E * = (u 1 , v 1 ) if and only if Condition (6) holds. Further, this equilibrium is unique.
Proof: Recall that since E * is stable for σ < σ * , a Hopf bifurcation occurs if J * has a pair of pure imaginary eigenvalues at σ = σ * and the real parts of these eigenvalues changes from negative to positive as σ increases through σ * . Note that when σ = σ * , the Jacobian matrix of J E * has a pair of purely imaginary eigenvalues λ = ±iω(σ ), where ω = √ det(J E * ). It is straightforward to verify that dRe(λ)/dσ | σ =σ * < 0. This means that the transversality condition holds.
The direction of the bifurcation is determined from the sign of the Lyapunov number . Using the formula in [1], we obtain Thus < 0 and hence a supercritical Hopf bifurcation occurs at σ = σ * and a stable periodic solution around the positive equilibrium E * will appear if σ is less than but close to σ * [18].
. We next establish the uniform persistence for System (3) by appealing to [12,Theorem 4.2]. The following theorem shows that System (3) is uniformly persistent whenever the positive equilibrium E * exists. Proof: The proof is standard following the idea in [12].
Our next result concerns the global stability of System (3). Theorem 2.5. If σ > σ 0 , then the equilibrium E 1 is globally asymptotically stable in X \ C 1 . If σ < σ 0 holds, then the positive equilibrium E * of System (3) is globally asymptotically stable in X 0 provided that one of the following condition is valid: Proof: If σ ≥ σ 0 , then System (3) has only the trivial equilibrium E 0 and the herbivorefree equilibrium E 1 . Note that E 0 is a saddle point, which attracts solutions only with (u(0), v(0)) ∈ C 1 . Since Equation (3) is a two-dimensional system, to show that E 1 is globally asymptotically stable in X \ C 1 , we only need to exclude the existence of a cycle. Suppose that System (3) admits a cycle (a periodic solution), then by the index theory for two-dimensional ODEs [23, Theorem 6.8.2], this cycle must enclose the equilibrium E 1 only. To achieve this, the cycle must intersects C 2 , i.e. the positive u-axis. For any initial condition (u(0), v(0)) = (u(0), 0) ∈ C 2 , there is a unique solution remaining in C 2 approaching E 1 . This contradicts with the uniqueness of the solution to System (3). Thus, E 1 is globally asymptotically stable in X \ C 1 .
Recall that if σ < σ 0 , then the positive equilibrium E * of System (3) exists. Consider the function where q(u) = u/(1 + σ u) and p(u) = u/(1 + βu). Note that q(u) and p(u) are positive and increasing for u > 0. It is clear that V(u, v) ≥ 0 and the equality holds only if (u, v) = (u 1 , v 1 ). Taking the derivative of V along the solutions of System (3) and noting that q(u 1 ) = κ, we obtain Note that if either Equation (13a) or (13b) hold, then v 1 < r. According to the property of the function, g 1 (u), is {E * } and all solutions in X 0 converge to E * . Thus E * is globally asymptotically stable.
Remark 2.6. If β ∈ (0, 1], then the local stability condition of E * is the same as the global stability condition of E * . If β > 1, then the local stability of E * only requires that σ ∈ (σ * , σ 0 ), while the global stability of E * needs a stronger condition, σ ∈ (σ * , σ 0 ) ⊂ (σ * , σ 0 ). There is a gap between the local stability condition and the global stability condition (See Figure 2). Though simulations suggest that E * is globally asymptotically stable as long as it is locally asymptotically stable, we are not able to fill this gap analytically.
The above theorem only gives the existence of the limit cycle for System (3). To get the uniqueness and stability of the limit cycle, we apply [15,Theorem 4.2] to obtain the following result. (15) satisfies Q(u) ≤ 0, then System (3) has a unique limit cycle which is globally asymptotically stable with respect to the set X 0 \ {E * }.
We take parameters r = 5.5, κ = 0.08, β = 1.5 and σ = 5. Then q 1 = −6, q 2 = −13.9, q 3 = 2.71, q 4 = 0.08, q 5 = −0.04 and Q(u) < 0 for u ∈ (0, 1). Thus Theorem 2.8 applies and System (3) has a unique limit cycle which is globally asymptotically stable with respect to the set X 0 \ {E * }. This is illustrated in Figure 1. In addition, there exists a heteroclinic orbit connecting the herbivore-free equilibrium E 1 and the limit cycle. That is because E 1 is a saddle point and the limit cycle is globally asymptotically stable. Remark 2.9. We derived the following sufficient but not necessary conditions. If q 2 < 0, q 3 < 0 and q 4 < 0, where q 2 , q 3 , and q 4 are defined in Equation (14), then System (3) has a unique limit cycle which is globally asymptotically stable with respect to the set X 0 \ {E * }.
We take parameter values: r = 5.5, κ = 0.08. Then according to Theorems 2.3 and 2.5, we can sketch the stability diagram in Figure 2 in the β-σ space, where Region V is numerically sketched by a shooting method. The unstable E * is surrounded by a stable limit cycle and there exists a heteroclinic orbit connecting the equilibrium E 1 and the limit cycle. Region I: E 1 is globally asymptotically stable; Region II: E * is globally asymptotically stable; Region III: E * is locally asymptotically stable; Region IV: System (3) has at least one limit cycle; Region V: System (3) has a unique limit cycle, which is globally asymptotically stable.
As shown in Figure 2, the herbivore-free equilibrium E 1 is globally asymptotically stable if (β, σ ) is located in Region I. If (β, σ ) is in Region II, then the positive equilibrium E * is globally asymptotically stable, while if (β, σ ) is in Region III, then E * is locally asymptotically stable. If (β, σ ) is located in Region IV, then E * is unstable and there exists at least one limit cycle, while if (β, σ ) is located in Region V, then E * is unstable and there exists a unique limit cycle, which is globally asymptotically stable with respect to X 0 \ {E * }.

Stability of the spatially homogeneous steady states
We note that System (2) admits some special solutions for which both time and space derivatives vanish. These solutions are referred to as homogeneous steady states. Clearly, the equilibria E 0 = (0, 0) and E 1 = (1, 0) of System (3) are two homogeneous steady states of System (2). If σ < σ 0 , then E * = (u 1 , v 1 ) is also a positive homogeneous steady state of System (2).
Next we determine whether or not these homogeneous steady states are stable. According to Casten and Holland [3], the (local) stability of a homogeneous steady state is determined by the eigenvalues of the characteristic equation: the homogeneous steady state is stable if all eigenvalues have negative real parts and is unstable if at least one eigenvalue has positive real part. Theorem 3.1. If σ > σ 0 , then the homogeneous steady-state E 1 of Equation (2) is locally asymptotically stable. If σ < σ 0 , E 1 is unstable and E * is a positive homogeneous steady state of System (2). Moreover, E * is locally asymptotically stable provided one of the following conditions is satisfied: Proof: The linearized operator of System (2) at a general steady-state (u, v) is given by Now we consider the characteristic equation of the linear operator L: where and e n , l n are coefficients. Substituting Equation (19) derives the relation Then the eigenvalues of the linear operator L are determined by the following eigenvalue problem: (J − n 2D ) e n l n = λ e n l n , n = 0, 1, 2, . . . .
We first consider the stability of the steady-state E 1 . At the steady-state E 1 of System (2), we have for n = 0, 1, 2, . . .. Then the trace and the determinant of J n (E 1 ) are given by Clearly, if σ > σ 0 holds, then T n (E 1 ) < 0 and D n (E 1 ) > 0 for n = 0, 1, 2, . . .. This shows that all eigenvalues of J n (E 1 ) have negative real parts, and hence all eigenvalues of the linear operator L have negative real parts. Therefore, E 1 is locally asymptotically stable for System (2) provided σ > σ 0 holds. Next we consider the stability of the steady-state E * which exists if σ < σ 0 . At the steadystate E * , we have The trace and the determinant of J n (E * ) are, respectively, given by and where D(σ ) > 0 (By σ < σ 0 ). Recall that E * is stable if all eigenvalues have negative parts provided the trace is negative. It is obvious that if Equation (16) holds, then T(σ ) < 0. This implies that T n (E * ) < 0 and D n (E * ) > 0 for n = 1, 2, . . ., and hence E * is locally asymptotically stable. (2), that is, spatial diffusion does not induce instability in System (2). However, as we will show later that spatial diffusion does induce some complex transient spatio-temporal behaviours.

Hopf bifurcation
Note that System (3) undergoes a supercritical Hopf bifurcation at σ = σ * . For System (2), besides the Hopf bifurcation value σ * , in this section we show that there may be an additional finite number of Hopf bifurcation values. For convenience, we introduce the following notation: let a and b be integers. Denote N(a) = {a, a + 1, a + 2, . . .},  N(a, b) = {a, a + 1, . . . , b − 1, b} for a < b, N(a, a) = {a} and N(a, b) We regard σ as the bifurcation parameter to investigate the occurrence of Hopf bifurcation at E * when σ < σ * . If for a positive integer n ∈ N(1), there exists a critical value σ H n at which J n has a pair of simple purely imaginary eigenvalues λ = ±iω(σ H n ) satisfying the transversality condition dRe(λ)/dσ | σ =σ H n = 0, then σ = σ H n is a Hopf bifurcation value at which spatially nonhomogeneous periodic solutions bifurcate from the unstable steadystate E * . Particularly, if n = 0, then σ H 0 = σ * is a spatially homogeneous Hopf bifurcation value.
With the above preparation, we now propose steps to geometrically determine the bifurcation values.  Step 1: Determine σ n,i , i = 1,2, n ∈ N 1 , and σ n by finding the intersection points of the graph of (σ ) for σ ∈ (0, σ * ) and the horizontal lines y = n, n ∈ N(1, I 1 ); Step 2: Plot the graph of n(σ ) = 4 D(σ )/d 2 2 obtained from solving T n (σ ) = 0 and D n (σ ) = 0 for n. If an intersection point determined in Step 1 is below the graph of n(σ ), then it is a Hopf bifurcation value at which a spatially nonhomogeneous periodic solution bifurcates from the unstable steady state. Of course, if 0 < κ < κ * , then σ =σ should be excluded.

Steady-state bifurcation
In this section, we carry out a steady-state bifurcation analysis to study the existence of nonhomogeneous steady states for System (2). A steady state of System (2) is a solution to the following semi-linear elliptic system: Steady-state bifurcations occur when the system has a zero eigenvalue. In particular, if there exists a critical value σ s n and an integer n ∈ N(0) such that the linearized system has a simple zero eigenvalue, then the stability of the homogeneous steady-state changes and a nonhomogeneous steady-state appears. That is, a steady-state bifurcation occurs at the bifurcation value σ = σ s n if and Thus a steady-state bifurcation value must be an element in the set 3 defined below 3 = {σ s n ∈ (0, σ * ) : (35) and (36) hold for some n ∈ N, }.
Next we derive the conditions on the existence of steady-state bifurcations. To this end, we define a constant M 4 by where Since the non-negative function ζ(σ ) is continuous on the finite interval [0, σ * ] with ζ(σ * ) = 0 and ζ(0) > 0, is well-defined . If then the set defined by is not empty. Thus for σ ∈ 4 , we can define Note also is well-defined if Equation (41) holds. We next give a sufficient condition on the nonexistence of steady-state bifurcation. Proof: If either Equation (16a) or (16b) holds, then it follows from the proof of Theorem 3.1 that T(σ ) < 0 and D n (σ ) > 0 for n ∈ N(0) and hence it is impossible for System (2) to undergo steady-state bifurcation. Solving D n (σ ) = 0 for n, we obtain n = P + (σ ) or n = P − (σ ).
If M 4 < 4d 1 /d 2 , then P + (σ ) and P − (σ ) are not defined and thus there does not exist σ such that D n (σ ) = 0. This shows that steady-state bifurcation cannot occur in this case.
Step 2 Plot the graph of y = (σ ) for σ ∈ S σ , where (σ ) is defined in Equation (28). If the intersection point (σ s n , n) obtained from Step 1 is not on this curve, then T n (σ s n ) = 0.
Step 3 Plot the graph of y = (1 − κσ is defined in Equation (25) and If the intersection point obtained from Step 1 is not located on this curve, then D n (σ s n ) = 0.
Step 4 Determine the bifurcation value. If the intersection point (σ s n , n) obtained from Step 1 is not on the curves given by Step 2 and Step 3, then σ s n ∈ 3 is a steady-state bifurcation value.

Transient spatio-temporal patterns induced by Hopf and steady-state bifurcations
Note that a spatially homogeneous Hopf bifurcation always occurs at σ * , while the occurrence of spatially nonhomogeneous Hopf bifurcations and steady-state bifurcations depends on the specific eigen-mode cos(nx).    Figure 6, no Hopf bifurcations or steady-state bifurcations occur as σ is varied. The homogeneous Hopf bifurcation value σ * is σ * ≈ 8.96. For σ < σ * , System (2) admits a spatially homogeneous periodic solution. A numerical solution is plotted in Figure 7 for σ = 6 < σ * . As shown in Figure 7, even with the given initial oscillations, the solution quickly evolves to a spatially homogeneous solution and converges to a stable spatially homogeneous periodic solution.
To explore the dynamics under case (ii), we choose the parameter values r = 5, β = 2, κ = 0.08, d 1 = 0.008 and d 2 = 0.1. In this case, there are two Hopf bifurcation values, σ H 1 ≈ 7.83, σ * = 8.5, but no steady-state bifurcation values. If σ is near σ H 1 , then there exists a spatially nonhomogeneous periodic solution with the eigen-mode cos(x), which is unstable. This would induce a transient spatio-temporal pattern. For instance, we take σ = 7.8, which is close to σ H 1 . Figure 8 depicts the short-time behaviour (left panel) and the long-time behaviour (right panel). As can be seen in Figure 8, the solution involves transient spatio-temporal oscillations for a short-time period before converging to a stable spatially homogeneous periodic solution. To clearly demonstrate the transient spatiotemporal oscillations, we plot the time series of u at location x = 2 in Figure 9 (left panel), in which the solution undergoes a few oscillations with small amplitudes first and then maintains the oscillations with large uniform amplitude corresponding to the stable spatially homogeneous periodic solution.
For case (iii), we choose the parameter values r = 5, β = 2, κ = 0.08, d 1 = 0.008 and d 2 = 0.15. In this case, as seen in Figure 10   7.33, σ * = 8.5, and 4 steady-state bifurcation values: σ S 4,1 ≈ 2.39, σ S 3,1 ≈ 3.90, σ S 4,2 ≈ 6.87 and σ S 3,2 ≈ 7.07. If σ is near a steady-state bifurcation value, then there exists a spatially nonhomogeneous steady state with the eigen-mode cos(nx) (here, n = 3 or n = 4), which is unstable. This would induce a different transient spatio-temporal pattern. As an example, we take σ = 7.05, which is close to σ S 3,2 . Figure 11 depicts the short-time behaviour (left panel) and the long-time behaviour (right panel). As can be seen in Figures 11 and 12, the solution follows a steady state involving only temporal oscillations for a short-time period before converging to a stable spatially homogeneous periodic solution. This can be effectively illustrated in Figure 9 (right panel), where u(t, 2) is plotted against time t.
The set of parameter values r = 3.5, β = 2, κ = 0.08, d 1 = 0.008 and d 2 = 0.25 corresponds to case (iv). As seen in Figure 13    Hopf bifurcation values, σ H 1 ≈ 8.15, σ * ≈ 9.17 and 10 steady-state bifurcation values: 3.44, 5.95, 6.88, 7.71, 7.82, 8.29, 8.62, 8.82, 8.94, 8.95. Note that this set of parameter values also corresponds to case (iii). Different from the above mentioned case (iii), here, σ H 1 is located between several steady-state bifurcation values. For simulation purpose, we pick σ = 7.8, which is near the steady-state bifurcation value 7.82 corresponding to the eigenmode cos(x) and is less than the first nonhomogeneous Hopf bifurcation value σ H 1 . We  note that for some initial conditions which are small perturbation of the bifurcated spatially nonhomogeneous steady state, the corresponding solutions may evolve first along the steady state, and then follow an unstable spatially nonhomogeneous periodic solution with a small amplitude and finally converge to the stable spatially homogeneous periodic solution with a larger amplitude. This is demonstrated in Figures 14 and 15.

Summary and discussion
In this paper we have considered a diffusive plant-herbivore model subject to Neumann boundary conditions. We have derived conditions under which the spatially homogeneous steady states, E 1 and E * , are stable. Hopf and steady-state bifurcation analyses are carried out when the steady-state E * becomes unstable. Geometric methods are provided to determine and locate the bifurcation values.
Depending on the choice of the parameter values, four scenarios involving Hopf and steady-state bifurcations are possible. That is, (i) neither Hopf bifurcation values σ H n nor steady-state bifurcation values σ S n exist for positive wave numbers n; (ii) there are Hopf bifurcation values but no steady-state bifurcation values; (iii) both Hopf bifurcation values σ H n and steady-state bifurcation values σ S n exist; (iv) only steady-state bifurcation values exist and no Hopf bifurcation value exists for any positive wave number n. These four scenarios result in very interesting transient spatio-temporal behaviour. The long-time behaviour is rather simple: the plants and the herbivores are spatially uniformly distributed with or without temporal oscillations depending on if the spatially homogeneous steady-state E * is unstable or stable. However, the short-time behaviour, i.e. the transient dynamics, can be very complex. It may involve transient oscillations both in space and in time, or only in one of them. These types of transient spatio-temporal patterns are influenced by the instability of the spatially nonhomogeneous periodic solutions, or that of the spatially nonhomogeneous steady states.
Keep all parameters fixed except σ , as seen from our simulations, as we decrease the value of σ , the steady-state E * may change from being stable to being unstable, and spatially nonhomogeneous steady states or spatially homogeneous and spatially nonhomogeneous periodic solutions may appear. In this sense, σ has a stabilizing effect. The diffusion coefficients also can influence the dynamics of the plant-herbivore interactions. For example, if the diffusion coefficient of the herbivore d 2 is very small, then no steady-state bifurcation would occur.