Exploration on dynamics in a ratio-dependent predator-prey bioeconomic model with time delay and additional food supply

1.   Introduction

There are some relationships between species such as predation, competition, mutualism and parasitism, while predation is a popular representative in ecosystem. Since Lotka ^[1] and Volterra ^[2] almost simultaneously proposed the Lotka-Volterra model of predation between two populations, many scholars have studied the dynamics of different types of predator-prey models (see ^[3,4,5,6]). They usually consider the response of the predator to the prey density changing through a functional response in the modelling of predator-prey model. A variety of models are generated due to the selection of different types of functional response functions. Holling type Ⅰ, Ⅱ and Ⅲ ^[7,8,9] were first proposed by Holling based on empirical field data from 1959 to 1966. The traditional predator-prey models assumed that the functional response depends only on the prey populations ^[10,11,12]. However, several biologists questioned functional response is determined solely by prey density. They suggested that predators might interfere with each other's foraging because predators share or compete for food with each other when predatos forage. Therefore, the functional response should depend not only on the prey density but also on the predator, see for instance ^{[13,14,15,16,17,18,19,20,21,22]}. The authors considered that the growth of population is on slow time scale (days or months) and established the ratio-dependent model. By virtue of Holling type Ⅱ function, Kesh et al.^[19] proposed a ratio-dependent model and obtained the sufficient condition for permanent co-existence. Xu and Chen ^[20] and Tang and Chen ^[21] studied the local and global stability of two classes of ratio-dependent system with delay, respectively. Fan and Li ^[22] proposed a delayed discrete ratio-dependent model and obtained the sufficient conditions for the permanence.

The study on providing additional food for the predator is an important topic in ecological modeling, which leads to more complex dynamics for the predator-prey model and has received extensive attentions. In recent years, many scholars have focused on the work of providing additional food for the predator to protect species or control pest, for details to see references [23,24,25,26]. Assuming that the quality of the additional food remains constant, Srinivasu and Prasad ^[23] found that the growth of the predator could be enhanced by the provision of extra food for the predator. Basheer et al. ^[24] showed that providing extra food for the predator is effective in eliminating pests. Samaddar ^[25] and Wu ^[26] investigated the Hopf bifurcation of the model which took into account the provision of extra food for the predator and time delay.

Motivated by the view of Arditi and Ginzburg ^[27] and Srinivasu et al. ^[28], a ratio-dependent model which involves the assumption that the additional food is provided for the predator was first proposed by Kumar and Chakrabarty ^[29]. This extra food is considered to be available uniformly within the ecological domain ^[28]. The model is as follows

Where $A$ is the additional food for the predator, the number of preys and predators at time $T$ is respectively represented by $N(T)$ and $P(T)$, here $N(0) > 0, P(0) > 0$. The intrinsic growth rate of the prey is expressed by $r$ and $m'$ stands for natural mortality of the predator, the carrying capacity of the prey is described by $K$. $e_1$ and $e_2$ denote the coefficients of predator capture on prey and additional food, $n_1$ and $n_2$ represent separately the nutritional value of the prey and additional food, $h_1$ and $h_2$ are handling time per predator per unit of the ratio $\frac{N}{P}$ prey biomass and per predator per unit of the ratio $\frac{A}{P}$ additional food biomass, respectively.

Another important aspect, the real ecosystem could not be denied the fact that the processes of biological development are not normally instantaneous on account of the interaction with environment or other species. Therefore, several biologists and mathematicians began to explore and study the time-delay models inspired by the dynamics of population models with time delay. From the biological and mathematical point of view, predator-prey models with delay ^{[30,31,32,33,34,35]} tend to exhibit more complex dynamics than those without that, such as local and global stability, Hopf bifurcation, limit cycle and so on. For example, the properties of the periodic solutions of the model with delay and impulse were obtained in ^[31]. Jiang ^[34] investigated the global asymptotical and the existence of Hopf bifurcation of a diffusive model with ratio-dependent function and delay. By studying the delay differential system, Tang and Zhang ^[35] showed that the stability of the model becomes unstable state when the time reaches a certain critical point, thus Hopf bifurcation occurs. Considering the maturation time of the prey, model (1.1) turns to the following version

It is widely known that human survival and development depend on natural resources. Biological resources have the most unique development mechanism due to renewability. Over-exploitation of biological resources will not only lead to a reduction in their quantity and quality, but also may drive species to extinction. Thus, taking the path of sustainable development becomes necessary to maintain ecological balance and meets the material needs of the people. From an economic perspective, the impact of harvesting effort on ecosystem guarantees that the net economic revenue is equal to the difference between total revenue and total cost which may influence on the harvesting activities. There has been a growing literature for the modeling and analysis of bioeconomic systems. The existence of steady state was discussed in a bioeconomic model with ratio-dependent ^[36]. The dynamics of the delay and diffusion terms of a bioeconomic plankton model was investigated by Zhao et al. ^[37]. In particular, an increasing number of bioeconomic models are expressed by differential-algebraic equations ^{[38,39,40,41,42,43,44]}. A prey-predator economic model with time-delay and stage structure was proposed by Zhang et al. ^[38]. Through detailed proof, the results show that the model undergoes three different bifurcational phenomena. Chen ^[41] established a biological economic model which suggested that the model undergoes flip bifurcation and Neimark-Sacker bifurcation.

To extend the model above, a bioeconomic model with time delay incorporating two differential equations and one algebraic equation is written as

From an economic perspective, the impact of harvesting effort on ecosystem is expressed as: net economic revenue (NER) = total revenue (TR) - total cost (TC). TR $= p_1q_1E(t)P(t)$ and TC $= wE(t)$, where $p_1$ is the price per unit predator biomass, $q_1$ and $E(t)$ stand for the harvesting capacity coefficient and harvesting effort of the predator, respectively. $w$ represents the cost per unit predator harvesting efforts. $Q$ means the NER. We only consider the harvesting of the predator, the catch rate function $h(t)$ is written as $h(t) = q_1E(t)P(t)$ which meets the catch-per-unit-effort hypothesis ^[45]. All the parameters in model (1.3) are positive. The initial conditions for model (1.3) are given $N_{[-\tau, 0]}\in C_+([-\tau, 0], \mathbb{R_+})$. If $P(0)$ is provided, $E(0)$ can be expressed as $Q/(p_1q_1P(0)-w)$ thus it is content with $P(0) > w/p_1q_1$.

Significantly, a real-life application is motivated by data from Katz ^[46]. We consider that prey $N$ is represented by barnacles Balanus balanoides and predator $P$ is expressed as snails Urosalpinx cinerea. The data strongly supports the ratio-dependent form $N/P$ much better than that only depends on prey $N$. Incorporating additional food $A$ (Ostrea gigas thunberg) into our model, our target is to address our understanding on the consequences of providing additional food for the predator on the dynamics, which brings out an explicit link between practical biological control and theoretical fruits. Besides, we investigate the optimal harvesting strategy of snails Urosalpinx cinerea due to its use as a kind of Chinese medical herbs.

The rest of this manuscript is emerged as follows. Section 2 is devoted to discussing the local steady-state and the Hopf bifurcation by investigating the characteristic equation. In Section 3, the main results including the properties of the bifurcating periodic solutions and the formulas of the direction of Hopf bifurcation are investigated. The optimal harvesting problem of model (1.3) without time delay is analyzed in Section 4. Four numerical examples are displayed in Section 5 to explain our findings. The brief and powerful biological discussion and conclusions are provided in Section 6. The paper end with Appendices A, B and C to show some proof and calculations.

2.   Local stability of equilibria and Hopf bifurcation

Model (1.3) is simplified and analyzed by using the following rescaling

After nondimensionalizing the model (1.3) for convenience, we obtain

where $c = \frac{e_1}{r}$, $\alpha = \frac{n_1h_2}{n_2h_1}$, $\xi = \frac{\gamma A}{K}$, $\gamma = \frac{n_2e_2}{n_1e_1}$, $b = \frac{n_1}{rh_1}$, $m = \frac{m'}{r}$, $q = \frac{q_1}{r}$, $p = p_1Ke_1h_1r$.

Through calculation, the internal equilibrium of the model (2.1) could be written as $P_0^*(x^*, y^*, E^*)$, where

and $x^*$ is the positive solution of $Ax^4+Bx^3+Cx^2+Dx+F = 0$ where

Remark 2.1. By Descartes sign rule ^[47], the possible positive roots of equation $Ax^4+Bx^3+Cx^2+Dx+F = 0$ are as follows.

(ⅰ) If any one of the conditions $(a)$, $(b)$, $(c)$ and $(d)$ holds, the equation exists one positive root

(ⅱ) If any one of the conditions $(e)$, $(f)$, $(g)$, $(h)$, $(i)$ and $(j)$ holds, the equation exists two positive roots or no roots

(ⅲ) If any one of the conditions $(k)$, $(l)$, $(m)$ and $(n)$ holds, the equation exists three positive roots or one positive root

(ⅳ) If the condition $(p)$ holds, the equation exists four positive roots or two positive roots or no roots

According to the following theorem, we know the conditions of the Hopf bifurcation occurs in model (2.1).

Theorem 2.1. If the assumptions $A_1+A_4 > 0, A_3 > 0, A_2+A_3 > 0$ and $A_2-A_3 < 0$ are satisfied, then when $\tau < \tau_0$, equilibrium $P_0^*(x^*, y^*, E^*)$ of model (2.1) is stable; equilibrium $P_0^*(x^*, y^*, E^*)$ of Eq (2.1) is unstable when $\tau > \tau_0$; when $\tau = \tau_0$, model (2.1) undergoes a Hopf bifurcation.

The detailed proving process of the Theorem is given in Appendix A.

3.   Stability and direction of the Hopf bifurcation

We investigate the direction of Hopf bifurcation and the properties of the bifurcating periodic solutions at the positive equilibrium in this section.

Theorem 3.1. If $\mu_2 > 0\ (\mu_2 < 0)$, the model undergoes a supercritical (subcritical) Hopf bifurcation at $P^*(x_1^*, y_1^*, E_1^*)$; moreover, if $\beta_2 < 0\ (\beta_2 > 0)$ then the bifurcating periodic solutions are stable (unstable); if $T_2 > 0\ (T_2 < 0)$, the periodic of bifurcating periodic solutions increase (decrease). Here

The detailed proof of the Theorem is given in Appendix B.

4.   Optimal harvesting

In order to achieve species persistence and profit maximization, optimal harvesting problem for model (1.3) without time delay is discussed in this section. We consider instantaneous annual discount rate $\delta$, we give a continuous time stream of revenues $\mathcal{J}$ as follows

The control constrain on this interval $0 < E(t) < E^{max}$. Here, $E^{max}$ stands for the maximum harvesting effort. The optimal prey and predator densities and the corresponding optimal harvesting effort are represented by $N^{optimal}$, $P^{optimal}$ and $E^{optimal}$, respectively. We are committed to calculating the optimal harvesting effort $E^{optimal}$ which satisfies the following formula

where $U$ presents the control set written as

The calculation procedure of $N^{optimal}$, $P^{optimal}$ and $E^{optimal}$ are shown in Appendix C.

5.   Numerical simulation

In this section, we will employ several specific examples to simulate the solutions to model (1.3), and verify the analytical results of the existence of Hopf bifurcation. Our analytical results illustrate that Hopf bifurcation arises in model (1.3) when the value of bifurcation parameter $\tau$ crosses critical values $\tau = \tau_k$, where

and $\omega_{10}$ is the root of equation ${\omega_1}^4+(A_1^2-2A_2-A_4^2){\omega_1}^2+(A_2^2-A_3^2) = 0$. Here, bifurcation parameter $\tau$ is the maturation time of prey. In order to illustrate the Hopf bifurcation of (1.3). We first choose the parameter values as follows: $r = 1.01$, $K = 0.91$, $e_2 = 1.01$, $h_1 = 0.99$, $h_2 = 0.01$, $n_1 = 0.85$, $n_2 = 0.98$, $A = 0.98$, $m' = 1.99$, $q_1 = 0.91$, $p_1 = 360$, $w = 1.01$, $Q = 0.51$. The scientific significance of these parameters are shown in Table 1, under the above parameters, model (1.3) becomes the following model (5.1)

Through maple software, we get $P_1(N^*, P^*, E^*) = (0.7191, 0.2939, 0.0054)$ and $P_1^*(x_1^*, y_1^*, E_1^*) =$ $(0.7903, 0.7248, 0.0108)$ are respectively the equilibrium of models (5.1) and (A.2) with above parameters.

Substituting the value of the positive equilibrium $P_1^*(x_1^*, y_1^*, E_1^*)$ in the expression of $A_1, A_2, A_3$ and $A_4$ to calculate the values as follows: $A_1 = 0.81$, $A_2 = -0.16$, $A_3 = 0.73$ and $A_4 = 0.79$ which surely satisfy the conditions of Theorem 2.1 (that is $A_1+A_4 = 1.60 > 0$, $A_2+A_3 = 0.56 > 0$ and $A_2-A_3 = -0.89 < 0$). In addition, we further substitute $A_1, A_2, A_3$ and $A_4$ to equation (2.8) to calculate $\omega_{10} = 0.74$, then substituting $\omega_{10}$ in the expression of $\tau_k$ to yield $\tau_0 = 1.86$. A series of calculations show that $g_{20} = -12.70+3.11i$, $g_{11} = -4.43-10.67i$, $g_{02} = 11.05-3.71i$, $g_{21} = -106.05+74.57i$, $c_1(0) = -97.31-43.88i$, $Re\{\lambda'(\tau_k)\} = 1.01$, $Im\{\lambda'(\tau_k)\} = 2.35$, $\mu_2 = 96.35 > 0$, $\beta_2 = -194.62 < 0$, $T_2 = 134.80 > 0$. Hence, by Theorem 3.1, the periodic increase, and the direction is supercritical, and the periodic solutions are stable. By Matlab we will reveal how the model can arise Hopf bifurcation by choosing different time delays $\tau$ in the following figures.

When $\tau = 1.63 < \tau_0 = 1.86$, model (5.1) satisfies the condition of Theorem 2.1. So $P_1(N^*, P^*, E^*) = (0.7191, 0.2939, 0.0054)$ is asymptotically stable. Figure 1 reveals the time evolution and the phase trajectory of model (5.1) for $\tau = 1.63$. It is obvious that model (5.1) converges to the positive equilibrium.

When $\tau = 1.87 > \tau_0 = 1.86$, model (5.1) satisfies the condition of Theorem 2.1. So $P_1(N^*, P^*, E^*) = (0.7191, 0.2939, 0.0054)$ is unstable, i.e., a family of periodic solutions are bifurcated from $P_1(N^*, P^*, E^*)$. The periodic solutions and its stability are shown in Figure 2.

Combined with Figures 1 and 2, we can see that the size of maturation time $\tau$ affects the dynamic behaviors of the model, see Figure 3. This means that if time delay $\tau$ exceeds critical value, bifurcation behavior occurs at the positive equilibrium $P_1(N^*, P^*, E^*)$ of model (1.3).

We choose the parameter values as follows: $r = 0.2$, $K = 0.1$, $e_2 = 0.7$, $h_1 = 0.95$, $h_2 = 0.8$, $n_1 = 0.1$, $n_2 = 3$, $A = 0.01$, $m' = 0.9$, $q_1 = 5$, $p_1 = 1$, $w = 0.01$, $\delta = 0.01$. $(N^{optimal}, P^{optimal})$ and $E^{optimal}$ as capture rate $e_1$ increase are shown by (a) in Figure 4. It shows that with the increase of capture rate $e_1$, the optimal prey population decreases initially and increases afterwards, optimal predator population gradually decreasing, and the optimal harvesting effort gradually decreases. Besides, we choose the parameter values as follows: $r = 1$, $K = 0.45$, $e_1 = 60$, $e_2 = 0.2$, $h_1 = 1$, $h_2 = 0.02$, $n_1 = 1.5$, $n_2 = 1$, $m' = 0.01$, $q_1 = 1.44$, $p_1 = 0.02$, $w = 1.01$, $\delta = 0.1$. $(N^{optimal}, P^{optimal})$ and $E^{optimal}$ as additional food $A$ increase are shown by (b) in Figure 4. It displays that with the increase of additional food $A$, the corresponding optimal prey and predator populations gradually increase, and the optimal harvesting effort gradually decreases.

6.   Discussion and conclusions

In our paper, based on a predator-prey model with additional food proposed by Kumar et al. ^[29], we incorporate time delay and harvesting into the above model and establish a differential-algebraic bio-economic model. Through a series of analytical analysis we obtain the results about the existence of positive equilibrium, the conditions of steady-state and Hopf bifurcation, and the direction of bifurcation as well as the stability of bifurcating periodic solutions. Firstly, the local stability condition of positive equilibrium of model (2.1) with time delay is $A_2-A_3 > 0$, while that without time delay are $A_1+A_4 > 0$ and $A_2+A_3 > 0$. Theorem 2.1 implies that the presence of the mature delay can destabilize model (2.1). When $\tau < \tau_0$, equilibrium $P_0^*(x^*, y^*, E^*)$ of model (2.1) is stable. The equilibrium is unstable when $\tau > \tau_0$; when $\tau = \tau_0$, model (2.1) has a Hopf bifurcation near the equilibrium. Besides, Theorem 3.1 displays that the direction of the bifurcation is determined by $\mu_2$, the stability of bifurcating periodic solutions is determined by $\beta_2$, and the period of the bifurcating periodic solutions is determined by $T_2$. Last, optimal harvesting strategies for model (1.3) without time delay are investigated. We achieve optimal harvesting equilibrium $(N^{optimal}, P^{optimal})$ and optimal harvesting effort $E^{optimal}$.

The numerical simulation of the model is carried out through Maple and Matlab software, and we find that the results of the numerical simulation agree with the theoretical results in Theorems 2.1 and 3.1. The influence of the bifurcation parameter $\tau$ on the dynamics of prey, predator and harvesting effort is discussed. The experimental results illustrate that the positive equilibrium of the model is locally asymptotically stable when bifurcation parameter $\tau < 1.86$, see Figures 1. Conversely, equilibrium $P_1(N^*, P^*, E^*)$ is unstable and there exist bifurcating periodic solutions which are stable {when bifurcation parameter $\tau > 1.86$, see Figure 2. The bifurcation diagram of (1.3) when $\tau = 1.86$ is presented in Figure 3, and the direction of bifurcation is supercritical. This also illustrates the influence of time delay $\tau$ which is one of the causes of population size fluctuation on the Hopf bifurcation of model (1.3), } i.e. the presence of the mature delay can destabilize model (1.3). Furthermore, we display optimal harvesting strategies in Figure 4, and depict the change of optimal harvesting equilibrium and optimal harvesting effort with respect to the rate of predator capture on prey $e_1$ and extra food $A$, respectively. The phenomenon in Figure 4 (a) is expected due to the stronger impact of the coefficient of predator capture on prey $e_1$. From Figure 4 (a), we know that optimal harvesting effort approaches zero when the coefficient of predator capture on prey $e_1 = 0.21$. This means that when the coefficient of predator capture on prey $e_1$ is beyond the level 0.21, harvesting predators should not be considered in order to prevent extinction. Figure 4 (b) shows that the optimal harvesting effort gradually decreases, while optimal prey and predator populations gradually increase, since additional food A distracts the predator from attacking the prey. At the same time, in order to prevent harmful prey from growing to uncontrollable numbers, the harvesting effort on predator should be reduced.

There are many topics that require us to refine and further analyze our model. On the one hand, the discrete time delay is considered in this study of Hopf bifurcation while more rich dynamic behaviors are generated by incorporating the distributed delays. On the other hand, the influence of random factors could not be ignored in which we introduce white noise into our model which involving the multiplicative random excitation and additive random excitation to analyze the stochastic P-bifurcation and stochastic D-bifurcation. The model is displayed as follows

where $\xi(t), \eta(t)$ represent the multiplicative random excitation and additive random excitation respectively. The above mentions make our model more interesting and complicated and remain to be done in the future.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors thank the editor and referees for their careful reading and valuable comments. The work is supported by the National Natural Science Foundation of China (No.12101211) and the Program for Innovative Research Team of the Higher Education Institution of Hubei Province (No.T201812) and the Scientific Research Project of Education Department of Hubei Province (No.B2020090) and the Program of Graduate Education Innovation of Hubei Minzu University (No.MYK2023042).

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix A

In order to linearize model (2.1), we consider the following linear transformation

Thus model (2.1) could be rewritten as

The internal equilibrium of the model (A.2) becomes $P_1^*(x_1^*, y_1^*, E_1^*)$, we transform the equilibrium into zero by linear transformation $x_1 = x_1^*+x_2$, $y_1 = y_1^*+y_2$, $E_1 = E_1^*+I(x_2, y_2)$, where

Model (A.2) is expressed as the following equation according to the above transformations

By calculation, the Jacobian matrix of the model (A.3) at $(0, 0)$ is represented as follows

The formula for the determination of the characteristic equation gives that $\lambda^2+A_1\lambda+A_2+(A_3+A_4\lambda)e^{-\lambda\tau} = 0$, where

When $\tau = 0$, the characteristic equation is as follows

By Routh-Hurwitz criterion, we know that $P_1^*(x_1^*, y_1^*, E_1^*)$ is locally asymptotically stable provided that the conditions $A_1+A_4 > 0$ and $A_2+A_3 > 0$.

When $\tau > 0$, the characteristic equation is yielded to

Suppose that $i\omega_1$ is a solution of (A.5), one gets

Separating the real and imaginary parts of Eq (A.6), we obtain

Furthermore,

Let $u = \omega_1^2$. Then, Eq (A.8) becomes

If $A_1 > 0, A_3 > 0, A_2+A_3 > 0$ and $A_2-A_3 < 0$, then Eq (A.8) exists a unique positive real root $\omega_{10}$, hence (A.5) has a pair of pure imaginary roots $\pm i\omega_{10}$ at the critical value $\tau_k > 0$, which could be denoted as

If $A_2-A_3 > 0$, (A.5) has no real roots and $P_1^*(x_1^*, y_1^*, E_1^*)$ is asymptotical stability for any $\tau > 0$.

We have to prove that $\frac{d}{d\tau}Re\lambda(\tau_k) > 0$ is true for guaranteeing the existence of Hopf bifurcation. Consider the following equation

Using a computation process similar to that of Tang ^[48] and Cooke ^[49], one yields

Since $A_2+A_3 > 0$ and $A_2-A_3 < 0$, we know that $f(u)$ has a unique positive real root $u_0 = \omega_{10}^2$. In addition it can be shown that $f'(u_0) > 0$ (that is $f'(\omega_{10}^2) > 0$). Combined with the above analysis, we have the following theorem.

Appendix B

Let $u_1(t) = x_2(\tau t)$, $u_2(t) = y_2(\tau t)$, $\tau = \tau_k+\mu\ (\mu\in\mathbb{R})$, then, $\mu = 0$ is the Hopf bifurcation value. The functional differential equation is considered as follows

where $u(t) = (u_1(t), u_2(t))\in\mathbb{R}^2$ and $L_\mu:C\rightarrow \mathbb{R}^2, F:\mathbb{R}\times C\rightarrow \mathbb{R}^2$. Let us define

and

where

in which $\mu_1 = (x_1^*+y_1^*+\alpha\xi)$, $\nu_1 = (pqy_1^*-w)$ and $\phi(\theta) = (\phi_1(\theta), \phi_2(\theta))\in C$.

There exists a $2\times2$ matrix function $\eta(\theta, \mu)$ $(\theta\in[-1, 0])$ such that

At this point, bounded variation function could be chosen as

Here, $\delta(\theta)$ is the Dirac function.

For $\phi\in C^1([-1, 0], \mathbb{R}^2)$, we define

Then, Eq (B.1) is equivalent to the following equation

The bifurcating periodic solution $u(t, \mu)$ of Eq (B.1) is affected by parameter $\varepsilon$, solution $u(t, \mu(\varepsilon))$ of Eq (B.1) exists amplitude $o(\varepsilon)$, nonzero exponential $\beta(\varepsilon)\ (\beta(0) = 0)$ and period $T(\varepsilon)$. Under the above assumptions, $\mu$, $T$, and $\beta$ have the following expansion form

Next, we compute the coefficients $\mu_2$, $T_2$, $\beta_2$. For $\psi\in C^1([0, 1], (\mathbb{R}^2)^*)$, assign the conjugate operator $\mathcal{A}^*$ of $\mathcal{A}$ is as follows

the adjoint bilinear form is given by

Obviously, $\pm i\omega_{10}\tau_k$ are both eigenvalues of $\mathcal{A}(0)$ and $\mathcal{A}^*$. Now we calculate the eigenvector of $\mathcal{A}(0)$ corresponding to $i\omega_{10}\tau_k$. Let $q(\theta) = (1, \Delta)^Te^{i\omega_{10}\tau_k\theta}$, then, $\mathcal{A}(0)q(\theta) = i\omega_{10}\tau_kq(\theta)$. Combining Eq (B.6) and the definition of $\mathcal{A}(0)$, it is easy to get

we further obtain

By a similar argument as that in the above part, the eigenvector of $\mathcal{A}^*$ corresponding to $-i\omega_{10}\tau_k$ could be achieved. Let $q^*(s) = \mathcal{D}(1, \Delta^*)e^{i\omega_{10}\tau_ks}$, we also have

From Eq (B.12) one gets

Hence,

Apparently, $q^*(s)$ $q(\theta)$ meets $\langle q^*(s), q(\theta)\rangle = 1$ and $\langle q^*(s), \bar{q}(\theta)\rangle = 0$. When $\mu = 0$, let $\mu_t$ be the solution of Eq (B.9), we set

Therefore,

here, the local coordinates of center manifold $C_0$ corresponding to $q$ and $q^*$ are expressed as $z$ and $\bar{z}$. Note that $W$ is real if $u_t$ is real, thus only real solutions are considered. For solution $u_t \in C_0$, since $\mu = 0$,

that is

where

Comparing Eqs (B.18) and (B.20), one has

where $q(\theta) = (1, \Delta)^Te^{i\omega_{10}\tau_k\theta}$, furthermore

Combining Eqs (B.3) and (B.4), we obtain

Comparing with the corresponding coefficients of the Eq (B.23), it yields to

In order to deduce $g_{21}$, now we have to calculate $W_{20}(\theta)$ and $W_{10}(\theta)$. It follows from Eqs (B.9) and (B.18) that

where

Consider on the center manifold $C_0$ close enough to the origin, we have

Substituting Eq (B.30) into (B.28) and comparing the coefficients on $z^2$ and $z\bar{z}$, it easy to see

For $\theta\in[-1, 0)$, from Eq (B.28) we can get

Comparing the coefficients of Eq (B.29) gives that

We have from Eqs (B.31) and (B.33) that

Noticing that $q(\theta) = q(0)e^{i\omega_{10}\tau_k\theta}$, then

where $M_1 = (M_1^{(1)}, M_1^{(2)})\in\mathbb{R}^2$ is a constant vector. By the same method as that in the above part we can calculate

where $M_2 = (M_2^{(1)}, M_2^{(2)})\in\mathbb{R}^2$ is also a constant vector.

Calculating constant vectors $M_1$ and $M_2$ is our main work in the following discussion. It follows from Eqs (B.28) and (B.29) that

where

and

According to the definition of $\mathcal{A}$ and Eq (B.31) we have

and

where $\eta(\theta) = \eta(0, \theta)$. Substituting Eqs (B.35), (B.36) and (B.37) into (B.31), we notice that

it is easy to see

Hence, we further get

and

From Eqs (B.35) and (B.36), we can determine the values of $W_{20}(\theta)$ and $W_{11}(\theta)$. $g_{21}$ be represented by the parameters of Eq (B.1). Therefore, we further calculate the values of $c_1(0)$, $\mu_2$, $\beta_2$ and $T_2$, as follows

Appendix C

Before finding the maximum value of $\mathcal{J}(E)$, the concrete Hamiltonian function is given

where $\lambda_i = \lambda_i(t), i = 1, 2,$ are adjoint variables. The condition that the Hamiltonian function $H$ satisfy is presented by

The adjoint variables $\lambda_1$ and $\lambda_2$ satisfy the following adjoint equations according to Pontryagin's Maximum Principle.

In order to find the solutions $\lambda_1(t)$, $\lambda_2(t)$ of Eq (C.6), we delete $\lambda_2$ in Eq (C.6) to obtain a second order differential equation with respect to $\lambda_1$

where

Then, we obtain $\lambda_1(t) = O_1e^{-\delta t}$, where $O_1 = \frac{C_1}{\delta^2-I_1\delta+I_2}$. Set up the equation in a similar way as above $\frac{d^2\lambda_2}{dt^2}+I_1\frac{d\lambda_2}{dt}+I_2\lambda_1 = C_2e^{-\delta t}$, we get $\lambda_2(t) = O_2e^{-\delta t}$, where $O_2 = \frac{C_2}{\delta^2-I_1\delta+I_2}$, $C_2 = -Ep_1q_1\bigg(\delta+\frac{rN}{K}-\frac{e_1^2h_1\frac{N}{P}}{(1+e_1h_1\frac{N}{P}+e_2h_2\frac{A}{P})^2}\bigg).$ Substituting $\lambda_2(t)$ into (C.5) yields

Therefore, the values of $(N^{optimal}, P^{optimal})$ and $E^{optimal}$ can be obtained by Eq (C.9) and solving the biological equilibrium.

The optimal harvesting effort at any time is determined by the following optimal harvesting solution

where $E^{min}$ is the minimum harvesting effort.