Hopf bifurcation in an age-structured prey-predator model with Holling Ⅲ response function

1.   Introduction

Ecosystem consists of varieties of living organisms and nonliving substances where the living organisms are not independent and have continuous interactions, such as mutualism, competition and predation ^[1,2,3]. It is obvious that the prey-predator interaction is of the greatest attention among other interactions in an ecosystem because of its universal existence and immense importance ^[4,5]. Moreover, due to the complexity of ecosystems, mathematical models have been considered as an alternative method of improving our knowledge of the ecological processes ^[6]. Such an approach goes back to the work of Lotka ^[7] and Volterra ^[8], who first proposed a prey-predator model in their pioneer work, which stimulated considerable theoretical researches on mathematical models of population dynamics.

Additionally, a growing evidence suggests that the functional response plays a crucial role in the interaction between the predator and prey. In recent years, various functional response functions, which describe the number of preys consumed by predator per unit time, have been used considerably in modelling population dynamics ^[9,10], for example, Holling ^[11] introduced several prey-predator models with three functional responses (named type Ⅰ, type Ⅱ and type Ⅲ) depending on different predation pressure. Specifically, the following mathematical model with Holling type Ⅲ functional response was analyzed by ^[12]:

where $V(t)$ and $U(t)$ denote the biomass of prey and predator at time $t$, respectively. Here, $r$, $c$, $k$, $d$ are positive constants that denote prey logistic growth rate, coefficient for the conversion that predator intake to per capital prey, carrying capacity of prey and mortality rate of the predator, respectively. The term $\frac{mV^{2}(t)}{b+V^{2}(t)}$ is Holling type Ⅲ functional response.

Recently, Holling type Ⅲ prey-predator models have been widely explored by incorporating some factors. Some studies have shown that a large refuge of prey can stabilize the positive equilibrium and protect prey from its extinction ^[13,14]. Moreover, a number of works have argued that diffusion-driven spatiotemporal patterns appear in diffusive prey-predator models with Holling type Ⅲ functional response ^[15,16,17]. In addition, many researchers are also interested in delay-induced dynamics in Holling type Ⅲ prey-predator models, and these studies suggested that delay can destabilize the positive equilibrium, and then the periodic oscillation occurs ^[18,19,20]. Apart from above researches, Holling type Ⅲ prey-predator models have been devoted by researchers to investigating the dynamics induced by other factors such as Allee effects ^[21,22] and competition ^[23,24].

Although all of above studies provide a better understanding of the dynamics of Holling type Ⅲ prey-predator models, they generally consider the fertility rate of a species as a constant for most eco-models. However, in natural ecosystem, a general mechanism for species is that a new reproducing generation needs time to mature to gain fertility, which is so called "age-dependent fertility". Actually, there exist many instances in real ecosystems which support that the fertility is age dependent, such as human, temora longicornis, hoverfly in nature ^[35]. The age-dependent fertility has been widely explored in prey-predator models, where age-dependent fertility is described by age-structure and there is a growing result suggested that age structure can influence the population dynamics ^[25,26,27]. In order to investigate the effects on model stability of several specific biological, delay causing mechanisms related to fertility, Cushing et al. ^[28,29,30] presented several prey-predator models and hypothesized age structure of predator can significantly influence its reproductive ability. These studies focused on the effects of mature period, of gestation periods and of age-related differentials ^[30]. In our investigation, the biological motivation is to understand further how mature period affect the model stability.

The study of the dynamics induced by age structure in prey-predator models is mainly focused on the boundedness of solutions, the existence and stability of steady states ^[31,32,33]. However, due to some factors such as noise and environmental factors, a constant population density may not exist in nature ^[34,35,36]. In recent years, Liu et al. investigated several prey-predator models with age structure ^{[37,38,39,40]}, and found that the age-dependent equilibrium may loses its stability via Hopf bifurcation and a non-trivial periodic oscillation phenomenon appears, which implied that the oscillation behaviors exists in age-structure models. Here, we are interested in investigating how age structure affect the changes of the population density. Motivated by the idea developed in Liu et al. ^[39], we present the following prey-predator model with age structure:

Here, $u(t, a)$ is distribution function of the predators over age $a$ at time $t$, and $\int_{0}^{+\infty}u(t, a){\rm d}a$ denotes the predator population density at time $t$. The function $\beta(a)$ is age-specific fertility function.

In natural ecosystems, predator individuals cannot reproduce in the immature stage, that is, $a\in (0, \tau)$, so the corresponding fecundity is zero. Obviously, when predator individuals are mature, the reproduction may occur. Thus, we have the following assumption:

Assumption 1.1. Assume that

where $\tau > 0$, $\beta^{*} > 0$.

Let $R = \int_0^{+\infty}\beta(a) e^{-da} {\rm d}a$, where $e^{-da}$ denote the survival probability. According to the results of ref. ^[25], the predator population is stable if $R\geq 1$, but the predator population tends to be extinct if $R < 1$. And the results of refs. ^[37,38] show that the results with the condition $R = 1$ are similar to ones with $R = M$, where $M > 1$ is a constant. Thus, we assume $R = 1$ in this paper, which implies $\beta^{*} = de^{d\tau}$.

In this paper, we propose a prey-predator model and investigate how the mature period effects the interaction dynamics between preys and predators. The rest of paper is organized as follows: In section 2, we express the model as a non-densely defined abstract Cauchy problem and obtain the positive equilibrium. The characteristic equation is obtained based on integrated semigroup theory, by which we derive the existence of Hopf bifurcation in Section 3. Some numerical simulations are given in Section 4. Finally, the paper ends with the conclusions in Section 5.

2.   Preliminaries

2.1.   Rescaling time and age

In this subsection, by the time-scaling $\bar{t} = \frac{t}{\tau}$ and the age-scaling $\bar{a} = \frac{a}{\tau}$, the change of variables $\bar{V}(\bar{t}, \bar{a}) = V(\tau\bar{t}, \tau\bar{a})$ and $\bar{u}(\bar{t}, \bar{a}) = u(\tau\bar{t}, \tau\bar{a})$, we can normalize $\tau$ in model (2). Then we have the following model:

and we denote the new function $\beta(a)$ as

Denote $V(t) = \int_{0}^{+\infty}v(t, a){\rm d}a$ in model (2), we can rewrite the first equation of model (3) as:

where

and $r = \theta-\mu$ represents the intrinsic growth rate of the prey, where $\theta$ is the birth rate of the prey.

Set $\chi(t, a) = \left(\begin{array}{c} v(t, a)\\ u(t, a)\\ \end{array} \right)$, model (2) can be rewritten as follows:

where $Q = \left(\begin{matrix} \mu & 0 \\ 0 & d\\ \end{matrix} \right)$, and $F(\chi(t, \cdot)) = \left(\begin{matrix} N(u, v)\\ \frac{cmV^{2}(t)\int_{0}^{+\infty}\beta(a)u(t, a){\rm d}a}{b+V^{2}(t)}\\ \end{matrix} \right).$

Considering the following Banach space:

where $\left\|\left(\begin{matrix} \alpha \\ \psi \end{matrix} \right) \right\| = \parallel \alpha \parallel _{\mathbb{R}^{2}}+\parallel \psi \parallel _{L^{1}((0, +\infty), \mathbb{R}^{2})}.$

And we define the linear operator $A_{\tau}: D(A_{\tau})\rightarrow X$ as

with $D(A_{\tau}) = (\{\textbf{0}_{\mathbb{R}^{2}}\}) \times W^{1, 1}((0, +\infty), \mathbb{R}^{2}) \subset X.$ The operator $H:\overline{D(A_\tau)} \rightarrow X$ is defined by

Due to the fact that $X_0: = \overline{D(A_\tau)} = \{\textbf{0}_{\mathbb{R}^2}\} \times L^1((0, +\infty), \mathbb{R}^2)$, the linear operator $A_\tau$ is non-densely defined.

Let

model (4) can be rewritten as the following non-densely defined abstract Cauchy problem:

model (5) has been studied including global existence and uniqueness of solution by ^[41,42].

2.2.   Existence of non-negative steady state

Define $\overline{x}(t) = \left(\begin{matrix} \textbf{0}_{\mathbb{R}^2}\\ \overline{\chi}(a) \end{matrix} \right) \in X_0$ is a steady state of model $(5)$, then we have

It follows that

From the first equation of model (6), we get

and

According to the second equation of model (6), we have

Obviously, model (5) always has a boundary equilibrium:

and the unique positive equilibrium:

with

if and only if $cm > 1$, $k-\sqrt{\frac{b}{cm-1}} > 0$.

In the remainder of our paper we assume that $cm > 1$, $k-\sqrt{\frac{b}{cm-1}} > 0$.

3.   Existence of Hopf bifurcation

3.1.   linearized equation

Here we will obtain the linearized equation of model (5) around the positive equilibrium $\bar{x}_{\tau}$. Taking the change of variable $x(t) = y(t)+\overline{x}_\tau$ into model (5), we get:

which is equivalent to

where

with

Next we define

satisfying $M(0) = 0$, and $DM(0) = 0$.

Therefore the model $(7)$ can be rewritten as:

3.2.   characteristic equation

In this subsection, we will get the characteristic equation of model (5) around the positive equilibrium $\bar{x}_\tau$. Denote

From the results of ^[43], we derive the following lemma.

Lemma 3.1. For $\lambda\in\Omega\cap\rho(A_\tau)$, we can obtain:

with $\left(\begin{matrix} \eta\\ \varphi \end{matrix} \right)\in X,$ and $\left(\begin{matrix} \textbf{0}_{\mathbb{R}^2}\\ \psi \end{matrix} \right)\in D(A_\tau).$

Furthermore, $A_\tau$ is a Hille-Yosida operator and satisfies

Denote $A_0$ is the part of $A_\tau$ in $\overline{DA(\tau)}$, for $\left(\begin{matrix}\textbf{0}_{\mathbb{R}^2}\\\psi\end{matrix}\right)\in D(A_0)$, we obtain:

where $\overline{A_0}(\psi) = -\psi^{'}(a)-\tau Q\psi, D(\overline{A_0}(\psi)) = \{\psi\in W^{1, 1}((0, +\infty), \mathbb{R}^2):\psi(0) = 0\}.$

Note that $\tau DH(\overline{x}_\tau):D(A_\tau)\subset X\rightarrow X$, is a compact bounded linear operator, we get

then we obtain

and from ^[42], we have

Consequently, we derive the following proposition.

Lemma 3.2. The linear operator $F_\tau$ is a Hille-Yosida operator, and its part $(F_\tau)_0$ in $\overline{D(F_\tau)}$ satisfies $\omega_{0.ess}((F_\tau)_0) < 0$.

Set $\lambda\in\Omega$, $(\lambda I-A_\tau)$ and $(\lambda I-F_\tau)$ are invertible and

then, we have $(I-\tau DH(\bar{x}))(\lambda I-A_{\tau}^{-1})$ is also invertible.

Let

which is equivalent to

We obtain

that is

Taking the formula of $DB(\overline{\chi}_\tau)$ into consideration, we get

where

and

From the proof of Lemma(3.5) in ^[44], we derive the following lemma.

Lemma 3.3. Under Assumption 1.1, the following results hold

● $\sigma(F_\tau)\cap\Omega = \sigma_p(F_\tau)\cap\Omega = \{\lambda\in\Omega:det(\Delta(\lambda)) = 0\}$,

● if $\lambda\in\delta(F_\tau)\cap\Omega$, we can obtain the formula for resolvent

where $\varphi(a) = (\Delta(\lambda))^{-1}(\chi+K(\lambda, \psi))e^{-\int_{0}^{a}(\lambda I+\tau Q){\rm d}l}+\int_{0}^{a}\psi(s)e^{-\int_{s}^{a}(\lambda I+\tau Q){\rm d}l}{\rm d}s$. On the other hand, we get

From (9) and (13), the characteristic at the positive equilibrium is:

where

Let $\lambda = \tau\xi$, we get

and

3.3.   the existence of Hopf bifurcation

In this section, we consider the parameter $\tau$ as a bifurcation parameter and investigate the existence of Hopf bifurcation by applying the Hopf bifurcation theory in ^[43]. We get

Obviously $\xi = 0$ is not an eigenvalue of Eq 15 when $p_{0}+q_{0}\neq 0$, and the real parts of eigenvalues of Eq 15 are negative for $\tau = 0$ when $p_1+q_1 > 0.$

Suppose that Eq 15 has a pure imaginary root $\xi = i\omega(\omega > 0)$, thus we get

Separating real and imaginary parts of the above equation, we have

It follows that

Set $\theta = \omega^2$, then Eq 17 turns into

Consequently, we can clearly see that Eq 17 has only one positive root $\omega_0 = \sqrt{\theta_0}$ when $p_0^2-q_0^2 < 0$.

From Eq 16, one can get that $g(\xi) = 0$ at $\tau = \tau_0^{k}, k = 0, 1, 2\cdots$ has a pair of purely imaginary roots $\pm i \omega_0$, where

and

Assumption 3.1. Assume that $cm > 1$, $k-\sqrt{\frac{b}{cm-1}} > 0$, and $p_1+q_1 > 0$, $p_0^2-q_0^2 < 0$.

Lemma 3.4. Under Assumptions 1.1 and 3.1, we have that Eq 17 has a simple root $\omega = \omega_0$, which given by

The proof is similar to ^[39], and hence we omit it.

Lemma 3.5. Under Assumptions 1.1 and 3.1, Eq 17 has a simple root $\xi(\tau) = \alpha(\tau)+i \omega(\tau)$ satisfying $\alpha(\tau_{0}) = 0$ and $\omega(\tau_{0}) > 0$, from ^[39] we have

Thus we conclude the following results.

Theorem 3.1. Under Assumptions 1.1 and 3.1, a Hopf bifurcation occurs around $\bar{x}_{\tau}$ for $\tau$ near $\tau = \tau_{0}$.

4.   Numerical simulations

In this section, we give two examples to illustrate the existence of Hopf bifurcation. The following parameter set are used

The age-specific fertility function becomes

By a simple computation, we can get that model (1) exists a positive equilibrium $E_{0}(1.36, 0.89)$ with the given parameter set. Based on software Python 2.7.15, we carry out some numerical simulations for model (1) using of the finite difference methods from ^[45], and the corresponding time series and phase diagram are shown in Figure 1, which indicates that the positive equilibrium of model (1) is stable. Then taking age structure into model (1), that is, the model (3), we can readily get $cm-1 = 1.4$, $k-\sqrt{\frac{b}{cm-1}}\approx 18.54$, such that there exists a positive steady state $E^{*}(0.54\tau e^{-0.5\tau a}, 1.46)$. Additionally, we obtain the conditions $p_1+q_1\approx 0.23$, $p_0^2-q_0^2\approx -1.01$ that guarantee the existence of critical value $\tau_0\approx 3.4$ from Eq 19.

Example 4.1. We give out an example corresponding to stability of the age-dependent equilibrium by choosing $\tau = 3$ in Figure 2. Obviously, the age-dependent equilibrium is locally asymptotically stable when $\tau = 3$ is below $\tau_{0}$, which is consistent with the theoretical results.

Example 4.2. In constrast, Figure 3 shows the oscillatory behavior induced by mature period. The results imply that when $\tau = 8$ is beyond $\tau_{0}$, the age-dependent equilibrium loses its stability via Hopf bifurcation and the corresponding periodic solution emerges. Under the same parameter set, Figure 1 shows that the positive equilibrium is stable in the absence of age structure. Obviously, the results shown in Figure 3 are induced by age structure.

In addition, we can clearly see that the existence of positive steady state depends on $a$ and $\tau$ from its expression. Figure 4 (a) shows that, under a fixed mature period value, the amplitude and period of periodic solutions have a decreasing tendency with the increase of age. In contrast, Figure 4 (b) shows that the mature period can enhance the prey-predator oscillation when the age is identical. Obviously, the mature period and age of species have a significant impact on the dynamic behavior of the prey-predator models.

5.   Conclusions

In this paper, we proposed a prey-predator model with age structure, and studied the effect of mature period on prey-predator dynamics. Specifically, the condition of the existence of non-negative steady state was derived, by which we found that the steady state exists whenever $cm > 1$ and $k-\sqrt{\frac{b}{cm-1}} > 0$. According to general results of the non-densely defined abstract Cauchy problem ^[41,42], we obtained the global existence, positivity and uniqueness of solution of model (3). By integrated semigroup theory, the characteristic equation is derived, by which we studied the effect of the mature period $\tau$ on stability of the steady state of model (3). Particularly, the theoretical results shown that when bifurcation parameter $\tau$ passes through the critical value, the age-dependent equilibrium loses its stability via Hopf bifurcation and the corresponding periodic solutions emerge, which indicates that the oscillatory behavior exists in changes of density of prey-predator populations. Actually, in real ecosystems, population oscillations is likely to be common due to many external factors ^[36], and the oscillatory behavior induced by mature period indeed exist in an ecological context.

Additionally, the numerical simulations showed the variation in the amplitude of prey-predator populations corresponding to the mature period and age. The numerical results suggested that the positive equilibrium is stable for model (3) in the absence of age structure. However, when age structure is incorporated, there exists a critical value $\tau_{0}$ such that the positive steady state loses its stability when mature period $\tau$ exceeding $\tau_{0}$, which means prey-predator population produces periodic oscillatory behaviors. And the numerical analysis shown that the population density tends to zero with the increase of age, which consist with the natural law of ecosystem. Furthermore, we observed that the mature period $\tau$ can affect the characteristics of periodic solutions such as amplitude and period. To be specific, when mature period $\tau$ exceeding the critical value $\tau_{0}$, the intensity of prey-predator population biomass oscillations is positively correlated with the mature period, while there is a negative correlation between age and the intensity of prey-predator biomass oscillation. The numerical analysis clearly shows that the mature period is the principle factor giving raise to oscillation, and the amplitude and period of periodic solution are greatly depending on $\tau$ and $a$.

In fact, there is a growing evidence that a constant population density may not exist due to the existence of some factors such as noise, physical factors and so on. In addition, the experimental results and field observation ^[34,35] have confirmed the existence of oscillations in reality. In present paper, we found that prey-predator population exhibits fluctuation under some conditions and mature period is considered as the factor causing population oscillation, which may be more reasonable compared with prey-predator model without age structure. Our results not only will help in further investigating the dynamics of population models but also provide an in-depth understanding of the oscillating coexistence of populations.

Acknowledgments

This work was supported by the National Key Research and Development Program of China (Grant No. 2018YFE0103700), the National Natural Science Foundation of China (Grants No. 61901303 and No. 61871293), Wenzhou basic argicultural science and technology project (N0. S20180006).

Conflict of interest

The author declares that there are no conflicts of interest.