Dirichlet problem for a diffusive logistic population model with two delays

In this paper, we investigate a diffusive logistic equation with non-zero Dirichlet boundary condition and two delays. We first exclude the existence of positive heterogeneous steady states, which implies the uniqueness of constant positive steady state. Then, we analyze the local stability and local Hopf bifurcation at the positive steady state. We show that multiple delays can induce multiple stability switches. Furthermore, we prove global stability of the positive steady state under certain conditions and obtain global Hopf bifurcation results. Our theoretical results are illustrated with numerical simulations.


1.
Introduction. In 1838, Verhulst proposed the following logistic equation to describe the self limiting growth of a biological population [17], where u is the population density at time t, r is the intrinsic growth rate, and K is the carrying capacity of the environment. Obviously, every nontrivial solution of (1) converges monotonically to K as t → ∞. To catch the observed sustained oscillations in the growth of many species in the nature, Hutchison [8] incorporated the effect of delays into (1). By assuming egg formation to occur τ units of time before hatching, the following delayed logistic model was proposed, The dynamics and properties of solutions to (2) have been studied very extensively, but there are still some open problems. It is well-known that if rτ < π 2 then the equilibrium u = K is locally asymptotically stable. However, so far, one can only prove the global stability of the positive equilibrium under the condition rτ ∈ (0, 1.5706]. It remains a conjecture that the local stability of the equilibrium implies its global stability; see Bánhelyi et al. [1]. When rτ > π 2 , the solutions of (2) oscillate about the carrying capacity K.
In (2), the regulatory effect only depends on the population at a fixed earlier time t − τ , rather than at the present time t. As pointed out by Kitching [9], the life cycle of the Australian blowfly Lucila cuprina has multiple time delay features. Based on this observation, (2) has been generalized to the following model with multiple delays, Some results on the linear stability of the positive equilibrium, oscillations of solutions, and the occurrence of limit cycles have been obtained; see, for example, [2,4,7,10,15,21] and references therein. In particular, Ruan [15] considered and proved that, (i) if a 1 ≥ a 2 (that is, the instantaneous dependence dominates) then the steady state u * = 1/(a 1 + a 2 ) is asymptotically stable for all delay τ ≥ 0; and (ii) if a 1 < a 2 then u * is asymptotically stable when τ ∈ [0, τ 0 ) and unstable . In the case where a 1 < a 2 , a Hopf bifurcation occurs at u * when τ passes through τ 0 . Recently, Yan and Shi [21] studied and observed multiple stability switches. Inherent in the above-mentioned study on logistic equations is the assumption that spatial properties are uniform and homogeneous, and that disturbances in carrying capacity are transmitted instantaneously throughout the population. In other words, spatial inhomogeneities and migrations are ignored. However, a fundamental goal of theoretical ecology is to understand how the interactions of individual organisms with each other and with the environment determine the distribution of populations and the structure of communities. In recent decades, the role of spatial effects in maintaining biodiversity has received a great deal of attention. Most theoretical approaches are based on or start with single species models. Among the simpler ways of introducing the effects of migration is through diffusion.
With heterogeneous environment in consideration, an extension of (3) is the following diffusive delayed logistic equation, where u(x, t) represents the density of a species at location x and time t, r(x) is the intrinsic growth rate at location x, and d is the random diffusion rate. When Ω is unbounded (usually Ω = R n ), traveling wavefronts [6,20,23] and asymptotic spreading [12] are investigated. When Ω is bounded, various boundary conditions can be imposed, which include Dirichlet boundary condition, Neumann boundary condition, and Robin boundary condition. In this case, the study focuses on the existence of steady states and their stability (see, for example, [3,4,5,13,16,22] and references therein). In many studies of Hopf bifurcations, the diffusion rate is often chosen as the bifurcation parameter. To the best of our knowledge, only Chen and Wei [4] and Su et al. [16] chose time delay as the bifurcation parameter. Su et al. [16] proposed the following diffusive model and analyzed the stability of the unique positive steady state solution. They also investigated the occurrence of Hopf bifurcation from this positive steady state solution, the direction of the Hopf bifurcation, the stability of the bifurcating periodic orbits, and global continuation of the Hopf bifurcation branches. Chen and Wei [4] studied a different diffusive model where λ, d, a i and τ i are all positive constants for i = 1, 2, . . ., N , Ω is a bounded domain in R n with a smooth boundary ∂Ω. They proved that if τ i is a "dominated delay", then it may affect the stability of the positive equilibrium. Moreover, Hopf bifurcation can occur when λ is near a critical value λ * . On the other hand, if τ i is not a "dominated delay", then it does not have any effect on the stability/instability of the positive equilibrium when λ is close to λ * . Here λ * is the principal eigenvalue of the following eigenvalue problem In both (5) and (6), the homogeneous Dirichlet boundary condition is chosen. Though instantaneous density dependence is contained in (5), there is only one delay. Equation (6) involves multiple delays, but it does not include the instantaneous density dependence. As we know, the absence of instantaneous density dependence in a growth model could make the prediction of population inaccurate. Moreover, as mentioned earlier, the growth of a population usually is affected by multiple delay factors. Therefore, it is desirable to consider both multiple delays and instantaneous density dependence. In this paper, we propose to study the following diffusive delayed logistic equation with nonhomogeneous Dirichlet boundary condition and two delays, where a, b, c > 0, and u * = 1/(a + b + c) is the carrying capacity. The remaining of this paper is organized as follows. In Section 2, we establish some preliminary results such as global existence, uniqueness, positiveness and boundedness of the solutions for our model system. In Section 3, we prove the nonexistence of the heterogeneous positive steady states using phase portrait method. In Section 4, we investigate the stability of the positive equilibrium, obtain multiple stability switches, and conduct local bifurcation analysis of the proposed model. In Section 5, we further develop unbounded global Hopf bifurcation theory via analyzing global continuity of periodic solutions. In Section 6, we conclude our paper with a brief discussion and numerical simulations.
2. Preliminary results. We start by considering the following diffusive logistic population model with two delays, where Ω is a connected bounded open domain in R N (N ≥ 1) with smooth boundary, and By constructing a pair of lower solution and upper solution, we establish the existence, uniqueness and boundedness of the solutions to (7).
Throughout the rest of this paper, we assume that N = 1 and the space Ω is normalized as Ω = (0, π).
3. Nonexistence of heterogeneous steady states. A steady state u(x) of (7) satisfies the boundary value problem, Let v(x) = u (x). Then equation (8) can be rewritten as It is easily observed that system (9) is a Hamilton system, where the Hamilton function is Then H(u, v) remains a constant along the solution curve of system (9). Obviously, system (9) has two fixed points (0, 0) and (u * , 0). It is easily verified that (0, 0) is a center and (u * , 0) is a saddle. Then we use the Hamiltonian function (10) to draw the phase portrait of (9) in Fig.1. Note that a steady state solution of (7) is equivalent to a trajectory of (9) that starts on the line u = u * at x = 0 and comes back on the same line at x = π. Thus we arrive at the following result regarding nonexistence of positive heterogeneous steady states of (7). Proposition 2. The solution u(x) ≡ u * is the unique positive steady state solution of (7).
Proof. Obviously, the solution u(x) ≡ u * is the unique constant solution of (7), which is also a positive steady state. To establish the uniqueness of the positive steady state, it suffices to exclude the existence of positive heterogeneous steady states of (7). Assume to the contrary that there exists a trajectory S = (u(x), v(x)) which starts on the line u = u * at x = 0 and comes back on the same line at Denote X to be the homoclinic orbit starting and ending at the point (u * , 0). By the Hamilton function (10), X consists of all points (u(x), v(x)) satisfying u(x) ≤ u * and . It is easily seen that v(0)v(π) < 0, otherwise, S and X would intersect with each other, which contradicts with the uniqueness of solution to (7). Then both v(0) and v(π) are positive or negative. Specifically, if both v(0) and v(π) are positive, then S should always stay in the first quadrant, which implies that u(x) is increasing and thus S cannot move back onto the line u = u * , a contradiction. Similarly, if both v(0) and v(π) are negative, then u(x) is decreasing, which leads to a contradiction. Therefore, we have excluded the existence of heterogeneous positive steady states of (7). The proof is complete. 4. Stability and Hopf bifurcation analysis. In this section, we study the stability of the constant steady state and the existence of Hopf bifurcation. We first transform the constant steady state to the origin via a translation U (x, t) := u(x, t) − u * . Then (7) can be rewritten as be the Hilbert space with usual inner product < ·, · > and where p = au * , q = bu * , r = cu * . The corresponding characteristic equation is Next, by using the delay τ as a bifurcation parameter, we analyze the stability switches at u * and the existence of Hopf bifurcation. Note that when τ = 0, the eigenvalues are which implies that u * of (7) is locally asymptotically stable if τ = 0. Thus a stability change at u * can only happen when one or more eigenvalues cross the imaginary axis to the right. It follows from dn 2 + p + q + r = dn 2 + 1 > 0 that 0 can not be an eigenvalue of (12). We then need to look for a pair of purely imaginary eigenvalues λ = ±iω with ω > 0 for some τ > 0. Substituting λ = iω into (12) gives By separating the real and imaginary parts, we have It then follows that cos ωτ = −q(dn 2 +p−r) Squaring and adding the above equations lead to (14) where z = ω 2 . Thus iω (with ω > 0) is a root of (12) if and only if the above equation has a positive root for some positive integer n. Note that the discriminant is Fix n ∈ N. We consider the following cases.
then (14) has a double positive root.
then (14) has two positive roots.
The equation (14) has no positive root if and only if one of the following three assumptions is satisfied. (H1) dn 2 + p ≥ 3r and q = dn 2 + p + r.
For simplicity, we assume that (H) r < 4d + p, q < 4d + p + r and q 2 < max{2(4d + p) 2 − 2r 2 , 8r(4d + p − r)}. It then follows from Lemma 4.1 that (14) with n ≥ 2 has no positive root. We are left to consider the case when n = 1. For illustration, we introduce the following four regions to classify the root distribution of (14); see Fig. 2.
The equation (14) has no positive root in A, exactly one positive root in B, and two positive roots in C 1 ∪ C 2 . Moreover, it has exactly one positive root in Now, we have the following result on the stability of u * .   (7) belong to the region A, then u * is locally asymptotically stable for any τ ≥ 0. If further dp 2 + p − q − r > 0, then u * is globally asymptotically stable.
Proof. If the parameters of model (7) belong to the region A, based on the arguments after (12), u * is locally asymptotically stable at τ = 0, and it does not change stability for all τ > 0 since 0 cannot be an eigenvalue of (12) and the equation (14) has no positive root. Therefore, u * is locally asymptotically stable for any τ ≥ 0.
To prove the global stability of u * , we need construct a Lyapunov functional V : C → R as follows.
The time derivative of V along a positive solution of (11) is given by It then follows from the boundedness of u(x, t) in Proposition 1, that is, U (x, t) + u * < 1/a, that Note that p = au * , q = bu * , r = cu * . Then dp 2 + p − q − r > 0 implies that da 2 u * + a − b − c > 0. Therefore, we have dV /dt ≤ 0. It is easily seen that the largest compact invariant subset of {U (x, ·) ∈ C : dV /dt = 0} is the singleton {U (x, t) ≡ 0}. Thus we have proved the global stability of u * for (7) by the LaSalle invariance principle [11].
It then follows from (13) that This completes the proof.
If the parameters of model (7) belong to the region B, then (14) with n = 1 has exactly one positive root z ≥ r 2 − (d + p) 2 + q 2 /2. Note that We have ω 2 + (d + p) 2 − r 2 > 0. It follows from the second equation of (13) that sin ωτ > 0; namely, ωτ belongs to either the first or the second quadrant. Solving the first equation of (13) yields It is readily seen that ω 2 + (p + d − r)(p + d − 3r) > 0, which, together with Lemma 4.3, implies that Thus the transversality condition holds and a Hopf bifurcation occurs at τ = τ j for each j = 0, 1, 2, · · · . When τ goes beyond a bifurcation point τ j , a pair of eigenvalues cross the imaginary axis from left to right. Recall that all eigenvalues have negative real parts at τ = 0. Thus the real parts of all eigenvalues remain negative when τ ∈ [0, τ 0 ). For j ∈ N and τ ∈ (τ j−1 , τ j ), there are exactly 2j eigenvalues which have positive real parts. Especially, u * is always unstable when τ > τ 0 .

5.
Global Hopf bifurcation analysis. In this section, we provide global Hopf bifurcation analysis on our model system by applying the global Hopf bifurcation theorem [18,19]. For simplicity, we only consider the case where parameters of model (7) belong to the region B ∪ L; the other cases can be investigated in a similar manner. Let z(t) = u(·, τ t) − u * and rewrite (7) as where is the nonlinear reaction function. Denote by {T (t)} t≥0 the semigroup generated by A T with zero Dirichlet boundary condition on (0, π). We then have It follows from z(t + ω) = z(t) that Iterating the above equation gives Note that the principal eigenvalue of A T is negative. We obtain T (t) → 0 as t → ∞. Especially, T (t + nω)z(0) → 0 as n → ∞. Consequently, On the other hand, if z(t) is a periodic solution of (25), it also satisfies (24). By [18,Chapter 6.5], the integral operator on the right-hand side of (25) is differentiable, completely continuous, and G-equivariant. It follows from Proposition 2 that u * is the unique positive steady state solution of (7). Note that 0 can not be an eigenvalue of (12) for any τ ≥ 0. Hence the assumption (H1) in [18, Chapter 6.5] is satisfied. On account of Theorem 4.4, the characteristic equation (12) has exactly one pair of purely imaginary eigenvalues ±iω when τ = τ j . Thus the assumption (H2) in [18, Chapter 6.5] holds. Following [18, Chapter 6.5], we introduce the local steady state manifold X) is a real isometric Banach representation of the group S 1 = {z ∈ C : |z| = 1} and E S 1 = {x ∈ E : gx = x for all g ∈ S 1 }. It is obvious that , ±iω is an eigenvalue of (12) if and only if τ = τ j . Thus the assumption (H3) in [18,Chapter 6.5] is also satisfied. Consequently, (u * , τ j , 2π/(ωτ j )) is an isolated singular point in M . Next, we define a closed subset Γ of X × R 2 + by Γ = Cl{(z, τ, T ) ∈ X × R 2 + : z is a nontrivial T -periodic soltuion of (23)}. By Theorem 4.4, the connected component of (u * , τ j , T j ) in Γ is nonempty and we denote it by C j (u * , τ j , T j ). Define It is readily seen that M = M * × (0, ∞). We then have the following result.
Proof. From Proposition 1, we have u(x, t) ≥ 0 and lim sup t→∞ u(x, t) ≤ 1/a for all x ∈ [0, π]. We claim u(x, t) ≤ 1/a for all (x, t) ∈Ω × R + . Otherwise, if u(x 1 , t 1 ) > 1/a for some (x 1 , t 1 ) ∈Ω × R + , then where T is the period of u(x, t). This contradicts the fact that lim sup t→∞ u(x 1 , t) ≤ 1/a. Thus 1/a is a uniform upper bound of u(x, t). Since 0 is a uniform lower bound, the proof is completed. Proof. We prove by contradiction. Assume that v(x, t) is a nontrivial periodic solution of (7) with period τ . Then we have First, we note that the solution z(t; z 0 ) to the ordinary differential equation with the positive initial condition z 0 > 0 converges to u * . We choose z 0 = max{u * , max [0,π] u 0 (x, 0)}. It is readily seen that z(t; z 0 ) and 0 formulate a pair of upper and lower solutions to (26). Hence system (26) has a unique solution v(x, t) such that 0 ≤ v(x, t) ≤ z(t; z 0 ). By strong maximum principle, we further have v(x, t) > 0 for any t > 0 and x ∈ [0, π]. Fix a t 1 > 0 and denote It is obvious that z(t; M 1 ) and z(t; m 1 ) formulate a pair of upper and lower solutions to the system satisfied by w(x, t) = v(x, t + t 1 ). Note that the system for w(x, t) is similar to (26) with a different initial condition w(x, 0) = v(x, t 1 ). Thus we have By squeeze theorem, we obtain This contradicts to the assumption that v(x, t) is a periodic solution of system (26) with period τ . Thus model (7) has no nontrivial periodic solution of period τ and the proof is completed.
We are now ready to analyze the structure of C j and prove global existence of periodic solutions of (7).
Theorem 5.4. Assume that (H) holds and the parameters of model (7) belong to the region B ∪ L. Let C j = C j (u * , τ j , T j ) be the connected component of (u * , τ j , T j ) in Γ. Then we have the following results.
(i) If j = l, then C j ∩ C l = ∅.
(ii) For each j ∈ N, the global Hopf branch C j is unbounded with unbounded τcomponent, bounded period component in (1/(j+1), 1/j), and bounded solution component in (0, 1/a]. (iii) For any τ > τ j with j ∈ N, there exist at least j periodic solutions for model (7).
Proof. Near the bifurcation point τ j with j ∈ N, it follows from the local Hopf bifurcation theorem that ωτ ∈ (2jπ, 2jπ + 2π). Consequently, the period T = 2π/(ωτ ) ∈ (1/(j + 1), 1/j). By Lemma 5.3, model (7) has no nontrivial periodic solution of period τ . Therefore, the scaled model (23) does not have nontrivial periodic solutions of period 1/j for any j ∈ N. Since the Hopf bifurcation branches are continuous, the periods on C j are bounded by 1/j and 1/(j + 1). Near the bifurcation point τ 0 , we have ωτ ∈ (0, π) and T > 1. Again, by continuity argument, the periods on C 0 are always greater than 1. Therefore, any two global Hopf branches C j and C l with j = l do not intersect. This proves (i).
For each j ∈ N, the periods on C j are bounded in the interval (1/(j + 1), 1/j). Recall from Lemma 5.2 that the periodic solutions on C j are bounded in the interval (0, 1/a). It then follows from Theorem 5.1 and the unboundedness of C j that the τ -component should be unbounded. This proves (ii).
Finally, (iii) is a simple consequence of (i) and (ii).
6. Conclusion and numerical simulations. In this paper, taken spatial heterogeneity into consideration, we proposed and studied a diffusive logistic population model with two delays (one is twice of the other) and instantaneous density dependence. To guarantee the existence of nontrivial homogeneous steady states, we imposed a nonhomogeneous Dirichlet boundary condition, namely, the population densities equal the carrying capacity on the boundary.
After showing that the carrying capacity is the only steady state, we studied its local stability. We found that delay can destablize the steady state and hence lead to Hopf bifurcation. Depending on parameter ranges, there may exist multiple sequences of Hopf bifurcation points and stability switches may occur. This is not observed in [4,16], where only one sequence of critical values is obtained. The reason is that instantaneous density dependence is not incorporated in the model in [4], while the model in [16] does not consider multiple delays. Our study shows that the combination of instantaneous density dependence and multiple delays may induce more complicated dynamics. We further considered the global continuation of these local bifurcation periodic solutions and proved the existence of multiple periodic solutions for some ranges of the time delay.
To exclude the existence of heterogenous steady states, we transformed the second order ordinary differential equation satisfied by the steady state into a Hamiltonian system and employed the technique of phase portrait. The success crucially depends on the special boundary condition. If one of the boundary conditions at the two end points is different from the carrying capacity, heterogeneous steady states may exist. We expect the model dynamics would be more complicated when this happens.
In the global Hopf bifurcation analysis, we proved that the delays on all but the first Hopf branch are unbounded. We conjecture that the delays on the first Hopf branch are also unbounded; or equivalently, the periods on this branch are bounded.
In our model, there are only two delayed terms and one delay is twice of the other, which means that there is always a "dominated delay". Therefore, the time delay always has effect on the stability of the only steady state, which agrees with the findings in Chen and Wei [4]. It will be interesting and challenging to extend the results in [4] to the more general case including instantaneous density dependence and multiple delays (not necessarily multiples of one).