A diffusive plankton system with time delay and habitat complexity effects under Neumann boundary condition

In this paper, we establish a delayed semilinear plankton system with habitat complexity effect and Neumann boundary condition. Firstly, by using the eigenvalue method and geometric criterion, the stability of the equilibria and some conditions for determining the existence of Hopf bifurcation are studied. Through analyzing the stability of positive equilibrium, we found that at the positive equilibrium the system may switch finitely many times from stable to unstable, then from unstable to stable, finally becoming unstable, i.e., the time delay induces a “stability switch” phenomenon. Secondly, the properties of Hopf bifurcation are derived by applying the normal form method and center manifold theory, including the bifurcation direction and the stability of bifurcating periodic solutions. Finally, some numerical simulations are given to illustrate the theoretical results, and a biological explanation is given.


Introduction
The plankton model is an important subject in marine biological systems. Chattopadhy et al. [1] showed that the delay of toxin release has a great impact on algal blooms. Wang [2] studied the toxic phytoplankton-zooplankton model and analyzed the effects of time delay and harvesting on the system. Sharma et al. [3] studied the plankton model with multiple delays. Roy et al. [4] established a nontoxic phytoplankton-zooplankton model and toxin-releasing phytoplankton-zooplankton model, respectively, and proved that nontoxic phytoplankton is beneficial to the growth of zooplankton, while toxin-producing phytoplankton is harmful to the growth of zooplankton. Furthermore, many scholars [5][6][7] have demonstrated that the toxins produced by phytoplankton can be used as a biological control quantity. Chattopadhayay et al. [8] pointed out that when the toxin rate exceeds the critical value, Hopf bifurcation occurs at the positive equilibrium. However, as long as the toxin rate is controlled close to the critical value, the system is stable at the positive equilibrium, so the harmful algal blooms are effectively controlled.
The theoretical study demonstrates [9][10][11][12][13] that the toxin released by phytoplankton may be a very strong regulation factor for the feeding rate of zooplankton. Habitat refers to the spatial scope of the environment where organisms appear, generally to the place where organisms live or the eco-geographical environment in which organisms live. The habitat complexity is the inhomogeneity of morphological characteristics within the structure itself and the heterogeneity of object arrangement in space. Research shows that the majority of population habitats are complex due to heterogeneity [14][15][16], for example, marine habitat become very complex in coral reefs, mangroves, sea grass beds, and salt marshes [17]. In lakes, the heterogeneity of habitat usually represents the vegetation depth and gradient diversity in coastal areas [18]. In addition, a large number of experimental studies show that habitat complexity reduces the encounter rate between predator and prey, thus reducing the predation rate [19][20][21][22][23][24]. Habitat complexity not only reduces the interaction between phytoplankton and zooplankton, but also reduces the available space of interacting species. Therefore, it is necessary to introduce habitat complexity into the plankton system.
Based on experiments and mathematical modeling, many scholars have established different mathematical models to describe the population dynamics [25][26][27][28]. Yang et al. considered the Holling type II plankton model with diffusion term in [29], and proposed the following model: in which = [0, lπ] (l > 0); P(x, t) and Z(x, t) represent the phytoplankton and zooplankton population densities at time t and distance x, respectively; d 1 and d 2 are diffusion terms; r is the intrinsic growth rate of phytoplankton; K indicates the maximum capacity of phytoplankton environment; c and e are the maximum capture rate and conversion rate of zooplankton; μ is the natural mortality of zooplankton population; ρ indicates the toxin intensity; f is the proportion of phytoplankton that can be caught by zooplankton, therefore, the phytoplankton with ratio 1f can aggregate to form a rough sphere, its surface area can be expressed as a function of ρP 2/3 . In this paper, we introduce the habitat complexity effect into system (1.1). Comparing with the processing time h, the habitat complexity is more likely to affect the attack coefficient β, therefore, we use β(1m) to replace β, where m (0 < m < 1) is a one-dimensional parameter used to measure the intensity of β. Assume that habitat complexity is homogeneous throughout the habitat, then the total amount of phytoplankton which is caught, V (P), can be expressed as where T s is the available search time and T is the total time. By calculation, we have Therefore, for system (1.1), the functional response function with habitat complexity effect is modified as Based on model (1.1), we introduce production delay, habitat complexity effect, and diffusion term to establish a toxic plankton model with Holling type II functional response function: where = [0, lπ] (l > 0). The rest of this paper is organized into sections. In Sect. 2, by analyzing the roots of the characteristic equation, we discuss the stability of diffusion system without delay at the equilibria (including boundary equilibria and positive equilibrium). In Sect. 3, we study the existence of Hopf bifurcation for the delayed diffusion system and the bifurcation direction, and the stability of periodic solutions is discussed by employing the center manifold and normal form theory. In Sects. 4 and 5, a biological explanation is given and some numerical simulations are carried out.

Stability of positive equilibrium point
We assume that system (1.2) has only one positive equilibrium, denoted as E * = (P * , Z * ). When τ = 0, we move E * = (P * , Z * ) to (0, 0). Making a transformationP = P -P * ,Z = Z -Z * , and omitting the bar, (1.2) becomes Defining the real Sobolev space , then system (2.1) can be written as an abstract functional differential equatioṅ where For system (2.1), the linearized equation at (m, 0, 0) iṡ where We use μ n = n 2 l 2 (n = 0, 1, 2, . . . ) to represent the nth eigenvalue of -ϕ xx = μϕ, ϕ x | x=0,lπ = 0. Define the linear operator It is easy to obtain that the eigenvalue of L(m) can be given by the eigenvalue of L n (m), and the eigenequation of L n (m) is The characteristic roots of (2.2) are Theorem 2.2 Assume (H 0 ) and (H 1 ) hold and K 2 < P * < K . Then the following conclusions are true: , then T n (m) > 0, D n (m) > 0, the roots of Eq. (2.1) have negative real parts, and system (2.1) is locally asymptotically stable at E * = (P * , Z * ); 1) has at least one root with positive real part, and system (2.1) is unstable at E * = (P * , Z * ).

Theorem 2.3
Assume (H 0 ) and (H 1 ) hold and 0 < P * < K 2 . Then the following conclusions are true: 1) has at least one root with positive real part, and system (2.1) is unstable at E * = (P * , Z * ).

Stability and bifurcation analysis of the delayed system
In nature, the change of population size is not only related to the current state, but also depends on the previous state. Considering the influence of the past state on the population size, we take the time delay as a bifurcation parameter to study the delay effect on the dynamic properties of the system, including the stability of positive equilibrium, the existence and direction of Hopf bifurcation, and the stability of bifurcating periodic solutions.
Thus the characteristic equation is equivalent to f n (λ, τ ) = λ 2 + A n λ + B n + C n e -λτ = 0, We make the following hypotheses: Proof We seek the critical value τ such that Eq. (3.4) has a pair of pure imaginary roots.
In summary, the conclusions are true, and the roots of Eq. (3.7) are Equation (3.6) has at least one positive root z ± n , ω ± n = z ± n .

Numerical simulations
In this section, we give some numerical simulations to verify the correctness of conclusions.
Therefore, the bifurcating periodic solutions are asymptotically stable, and the period increases.