Harvesting a population model with Allee effect in a periodically varying environment

1.   Introduction

The population growth behavior has been described through numerous models. For instance, the logistic model ^[1,2] displays the population evolution from the initial density to the carrying capacity limit, no matter how small the initial population is. However, some populations only evolve when the initial population is large enough, otherwise the population is extinct. Some reasons for this behavior have been described by Courchamp et al. ^[3] in several categories, such as genetic diseases, demographic stochasticity, etc. Such a phenomenon is called a strong Allee effect ^[4,5]; for further details, see ^[6,7,8] and the references therein. There are numerous examples in the real world that adapt the feature of the Allee effect, which have been addressed in ^[9,10,11,12]. Mathematically, the simplest model that displays this feature ^[13] is given by

where $N(T)$ is the population density at time $T\geq0$, whereas $P, K, N_0$, and $S$ are positive constants of per capita growth rate, carrying capacity, the initial population and Allee threshold, $0 < S < K$, respectively. Model (1.1) clearly has the property that when the initial population is greater than the Allee threshold, $0 < S < N_0$, the solution travels to the carrying capacity, describing the population at the equilibrium. In contrast, when $0 < N_0 < S$, the solution decreases to zero, showing the extinct situation of the population. Since the solution of such model can only be provided in an implicit form, Idlango et al. ^[14] obtained an analytical approximate expression for the population evolution and showed a very good comparison with the numerical data.

Note that adding a constant effort harvest term that is proportional to the population density, $HN(T)$, into (1.1) gives the logistic population model subjects to the harvesting and the strong Allee effect ^[15] as follows

describing the situation of subcritical and supercritical harvesting under the restriction of the Allee effect.

It must be mentioned that the qualitative properties of the population models are regularly influenced by the seasonal and even violent changes in the environment, such as climate change, food sources, hunting, etc ^[16]. Thus, it is seen that within the progress of the mathematical biology, the time-varying coefficients are often incorporated with the models. For example, the dynamical models in a view of a periodic environmental change have been investigated by Rosenblat ^[17] and Legović et al. ^[18]. More recently, Alharbi ^[19,20], Alsharidi et al. ^[21], and Idlango et al. ^[22,23] considered single species models subjects to variable coefficients. Accordingly, the model (1.2) can be written in terms of variable coefficients as follows

Clearly, when the time-varying parameters are adopted, as in (1.3), the population models turn to nonautonomous systems and the exact solution is far of achievement, so the numerical approximate method must be resorted to; see, for example, the analysis of those in ^[17,18] among others ^[24,25]. However, the numerical method often responds to particular parameter values, which leads to providing a limited information about the population evolution, as in ^[25]. As a result, it is essential to consider the analytical approximate method in order to obtain a general representation that is capable to describe the influence of the time variation in the model coefficients and its impacts on the population evolution for various general scenarios.

In many cases, it is possible to consider slow fluctuations in the environmental factors, such as climate changes, adequate food, etc. This may lead to proposing a slowly varying population model. Ludwig et al. ^[26] used such observation in the analysis of the spruce budworm, but leaving their calculation unsolved analytically. Thus, Shepherd et al. ^[27], Grozdanovski et al. ^[28], and Alharbi ^[19,20] considered a multi-scaled perturbation method to investigate population models in a slowly varying environment for various cases. Accordingly, the analytical approximate expressions have been obtained and validated numerically. The comparison between the two methods shows a very good agreement. However, in their analysis, the associated coefficients functions must vary in a slow time scale compared to the intrinsic time scale of the overall population, otherwise the considered hypothesis is not valid in obtaining the analytical expression of the population evolution. In addition, utilizing the multi-scaled perturbation method often provides asymmetric rate of convergence between the leading and higher order differential equations, which must be aligned by eliminating the disturbed terms, arising in the higher order equations.

Therefore, in this paper, we aim to investigate the single species population model (1.3) that exhibits the Allee effect and a slowly periodic fluctuation of the harvesting rate. The assumption of the harvesting periodicity may arise in response to relatively slow fluctuations of the environment. This allows us to utilize the regular perturbation approach based on the small variation to obtain the analytical approximate solutions for the population evolution in the situation of subcritical and the supercritical harvesting with respect to the Allee threshold condition. Such analytical expressions will successfully be able to show the harvesting variation on the same time scale of the overall population. Moreover, the rate of convergence for the leading and higher order population expressions is symmetric without any elimination of the disturbed terms. In addition, due to dealing with the nonlinear system, a numerical computation will be obtained using a fourth-order Runge-Kutta technique in order to validate the analytical approximation. The comparison between the two methods is expected to show that when the time progresses, the analytical expressions will display a very good convergence with the numerical data, while such accuracy falls off when the time approaches the initial density. The finding results are appropriate for a wide range of parameter values, which provides an advantage of avoiding extensive recalculations while displaying the population evolution.

2.   Preliminaries

Since we are only interested in the influence of the periodic (or seasonal) fluctuation of the harvesting, which might be caused by a disease, a fishing, or a hunting, on the population behavior, the coefficients $P(T), K(T)$, and $S(T)$ in (1.3) are assumed fixed. So, we arrive at

In what follows, we render the system (2.1) dimensionless. By expressing the nondimensional time, population, and harvesting effort scale by $t = PT$, $N(T) = Kn(t)$, and $H(T) = H_0\eta(t)$, respectively, (2.1) can be written in a dimensionless form as

where $\eta(t)$ fluctuates periodically and is given by

and $\sigma$, $\alpha$, and $\mu$ are given by

respectively, $\epsilon$ is a dimensionless, small, and positive parameter, and $H_0$ is the baseline of the harvesting function.

As the model (2.2) stands, two situations that depend on the harvesting have been identified. One is the subcritical harvesting, where the population evolves to equilibrium, or is extinguished due to the Allee effect. The other is supercritical harvesting, where the population is extinct for any initial density $\mu$.

3.   Perturbation approach

It is important to note that, for the general variation, assuming at least one of the model coefficients varies with time, the system (2.2) can only be solved using numerical techniques. However, applying (2.3) in (2.2) and assuming $\epsilon$ is sufficiently small, this permits obtaining the analytical approximate expressions for the population evolution model using the perturbation approach. Accordingly, we propose the perturbation expansion based on $\epsilon$ as follows

and substituting (3.1) into (2.2) obtains, on equating like powers of $\epsilon$, a leading order term

a higher order term

and the initial conditions

where $\psi$ and $\delta$ are defined by

respectively. Therefore, the leading order term (3.2) is an autonomous differential equation, whereas the higher order term (3.3) describes a slowly varying population. Since $\alpha$ measures the ratio of the baseline of the harvesting effort, $H_0$, to the growth rate, we conclude from the second of (3.5) and (3.2) that when $\alpha < \sigma\psi^2-1$, $\delta > 0$, so the harvesting is subcritical. In contrast, when $\alpha > \sigma\psi^2-1$, $\delta < 0$, so the harvesting becomes supercritical. Consequently, we consider these two cases.

4.   Subcritical harvesting

In the following, we consider $\delta > 0$, so the harvesting is subcritical. Note that sometimes the initial density, $\mu$, is driven below the critical threshold, so it will, in this section, be expected that the population either survives to a positive limiting state or is extinguished, depending on the $\mu$ value.

With the application of the initial condition, given in the first of (3.4), the differential equation (3.2) may be solved implicitly in the form

where

Straightforward, applying the stability analysis to the leading order term (3.2) displays the behavior of the solution (4.1), which is governed by the initial density, $\mu$. In order to discuss that, we first introduce the critical points as follows

$(1)$ $cp_1 = 0$, which is stable.

$(2)$ $cp_2 = \frac{\sigma\psi-\sqrt{\sigma\delta}}{\sigma}$, which is unstable and acts as breakpoints ^[29].

$(3)$ $cp_3 = \frac{\sigma\psi+\sqrt{\sigma\delta}}{\sigma}$, which is stable.

Thus, when the initial population obeys $cp_2 < \mu < cp_3$, $n_0(t)$ increases to the surviving limiting state, $cp_3$, as $t\rightarrow \infty$, whereas it decreases to the limit when the $\mu$ value is beyond $cp_3$. In contrast, when the $\mu$ value lies below $cp_2$, the population tends to extinguish, i.e., $n(t)\rightarrow0$ as $t\rightarrow \infty$; see Figure 1.

Note that $n_0(t)$ appears implicitly in (4.1), so our approach is apparently not able to obtain the expression of the higher order term, $n_1(t)$. However, it can be observed that as $t\rightarrow \infty$, the term $e^{-\phi t}\rightarrow 0$. Such an observation shows a possibility of obtaining the asymptotic expression of $n_0(t)$ by defining the solution of (4.1) as a power series in terms of a small $e^{-\phi t}$. Once this achieved, the linear differential Eq (3.3) can, in principle, be solved utilizing an integrating factor together with the application of the initial condition, given in the second of (3.4). Accordingly, we express $n_0(t)$ as a Poincar$\acute{\text{e}}$ expansion in $e^{-\phi t}$

where the subscript $s$ denotes association with the survival case. Substituting (4.2) in (4.1), expanding in power of $e^{-\phi t}$, and equating like powers of $e^{-\phi t}$ yields the coefficients $f_0, f_1, \cdots$. Thus, with the choice $f_0 = \frac{\sigma\psi+\sqrt{\sigma\delta}}{\sigma}$, the expansion for $n_{0s}(t)$ when $n_{0s}(t)\rightarrow cp_3$ as $t\rightarrow \infty$ becomes

where $\chi$ is given by

determining the leading order term when $\mu > cp_2$.

Similar to the above but more precise calculation, when $\mu < cp_2$, $n_{0}(t)\rightarrow 0$ as $t\rightarrow \infty$. Thus, the expansion of $n_{0}(t)$ is given by

where $\omega$ is given by

and the subscript $d$ denotes association with the extinction case. Note that in the expansion (4.3), the convergence rate to $cp_3$ (survival) depends on $\phi$, while in the expansion (4.5), the convergence rate to $cp_1$ (extinction) depends on $(1+\alpha)$.

Therefore, the solution of the leading order term $n_0(t)$ has split into two parts: $n_{0s}(t)$, given in (4.3), is to describe the surviving population, and $n_{0d}(t)$, given in (4.5), is to represent the extinct population due to driven $\mu$ below the critical threshold. Such division will be maintained while solving the higher order terms.

Substituting the leading terms (4.3) into (3.3) yields

Considering the leading terms indicated in (4.7) and solving, we obtain, on applying the initial condition, given in the second of (3.4),

determining the higher order expression of the surviving population in terms of the periodic function $\gamma(t)$.

In a like manner, substituting the leading terms (4.5) into (3.3), we arrive at

As earlier, preserving the leading terms mentioned in (4.9) and solving, we obtain, on applying the initial condition, given in the second of (3.4),

determining the higher order expression of the extinct population.

To recapitulate, considering the subdivisions of the leading and higher order terms, the perturbation expansion (3.1) for the surviving population at the positive limiting state becomes

whereas the expansion for the extinct population due to the Allee effect becomes

where $n_{0s}(t), n_{1s}(t), n_{0d}(t)$ and $n_{1d}(t)$ are given by (4.3), (4.5), (4.8), and (4.10), respectively.

5.   Supercritical harvesting

Considering $\alpha > \sigma\psi^2-1$, $\delta$ is always negative, the harvesting is supercritical. Accordingly, we redefine $\delta$, given in the second of (3.5), in such a way that $\delta = -\delta_e$, where $\delta_e$ is given by

and the subscript $e$ denotes association with supercritical harvesting.

Thus, considering (5.1) together with the first of (3.4), the solution of (3.2) yields

where

Again, $n_{0e}(t)$ only appears implicitly. So, applying analogous reasonings of that in Section 4 but when $n(t)\rightarrow 0$ as $t \rightarrow \infty$, we define the asymptotic expansion of $n_{0e}(t)$ as a power series in terms of a small $e^{-(1+\alpha)t}$ as follows

where $\tau$ is given by

determining the leading order term of the supercritical harvesting.

In a like manner, substituting (5.4) into (3.3), we arrive at

Retaining the leading terms mentioned in (5.6) and solving, we obtain, on considering the initial condition given in the second of (3.4),

determining the higher order term of the supercritical harvesting.

Therefore, the perturbation expansion (3.1) regarding the supercritical harvesting becomes

where $n_{0e}(t)$ and $n_{1e}(t)$ are given by (5.4) and (5.7), respectively.

6.   Results and discussion

The analytical approximate expressions for the logistic population model that exhibits the Allee effect together with the slowly periodic harvesting term are readily obtained and the expansions (4.11), (4.12), and (5.8) represent the solutions. These expansions have the advantage of not requiring specific parameters values to be provided, which leads to avoiding the extensive recalculations while displaying the population behavior. It is important to note that our expansions are limited in the first two terms due to displaying a very good convergence with the numerical solution and avoiding the swell in the expression. In addition, since the asymptotic expansions of the leading order terms, given in (4.3), (4.5) and (5.4), are constructed in terms of small exponential terms, we expect that when $t$ is close to the initial density, the perturbation expansions (4.11), (4.12), and (5.8) diverge from the numerical solutions, whereas these expansions provide remarkable agreement as $t\rightarrow \infty$.

For $\alpha < \sigma\psi^2-1$, the harvesting is subcritical and the population survives at the equilibrium state when $\mu > cp_2$ and the expansion (4.11) displays the evolution. On the other hand, when the initial density is below the Allee threshold (i.e., $\mu < cp_2$), the population is extinguished and the expansion (4.12) exhibits the population behavior.

When $\alpha$ is increased beyond the limit $\sigma\psi^2-1$, the harvesting is turned to supercritical, which leads to declining the population to zero in a finite time, as shown in the expansion (5.8).

Note that $\epsilon = 0$ corresponds to the exact solution (no slowly varying term), and the leading order terms (4.1) and (5.2) describe the population evolution in the situations of subcritical and supercritical harvesting, respectively, including the Allee threshold condition. However, when $\epsilon > 0$, (4.11), (4.12), and (5.8) describe the slowly varying population behavior.

As mentioned earlier, the implicit representation of the leading order terms of the subcritical and supercritical harvesting make the explicit expressions of the higher order terms unattainable. However, for a particular case, setting the periodic function $\gamma(t) = \sin t$, the solution of (4.8), (4.10), and (5.7) yields

determining slowly varying expressions for subcritical harvesting (survival and extinct case due to the Allee threshold) and supercritical harvesting, respectively

From (6.1)–(6.3), note that for a large class of time, we arrive at

and, hence, the perturbation expansions (4.11), (4.12), and (5.8) become as $t\rightarrow \infty$

describing the population behavior at the limiting states.

It must be pointed out that the leading and higher order terms of our perturbation expansions have a symmetric convergence rate, whereas in the analysis of Alharbi ^[19,20] and Idlango et al. ^[22,23], they needed to be aligned by eliminating the disturbed terms, arising in the higher order terms.

Figures 3 and 4 were generated using the expansions (4.11), (4.12), and (5.8), respectively, together with the numerical data obtained using a fourth-order Runge-Kutta technique. It is clearly seen that the convergence between the two solutions is very good for $t\rightarrow \infty$, but such convergence falls off at $t$ close to zero as expected from the asymptotic expressions above. In addition, Table 1 is obtained by calculating the absolute difference between the approximate expressions and the numerical data as $t\rightarrow \infty$. This shows a remarkable discrepancy between the analytical approximate and numerical solution when $t$ is small, but when $t$ progresses, the two solutions match each other. Graphically, Figure 2(a), (b) confirms the result of Table 1 by showing a monotonically decreasing discrepancy curve for subcritical and supercritical harvesting as the time increases.

Figure 3(a), (b) regards the population evolution in the case of subcritical harvesting. They show that when $\mu > cp_2$, the population survives to the equilibrium state, $cp_3$, from below when $cp_2 < \mu < cp_3$ (Figure 3(a)) and from above when $\mu > cp_3$ (Figure 3(b)).

In contrast, Figure 4(a) displays that when $\mu$ locates below the critical Allee threshold, $\mu < cp_2$, the population is instantly extinct in finite time. Figure 4(b) relates to supercritical harvesting. It can be seen that for any initial population, $\mu$, the population dies out as $t\rightarrow \infty$.

Figure 5(a), (b), which was generated using the expression of the equilibrium given in (6.4), shows the influence of the associated parameters $\sigma$ and $\alpha$ on the population behavior. Figure 5(a) shows that fixing the value of $\alpha$ while increasing $\sigma$ leads to increasing the variation of the population at the equilibrium. In contrast, in Figure 5(b), when fixing the value of $\sigma$, the population density reduces as $\alpha$ increases.

7.   Conclusions

In this paper, we investigate the single species population model that exhibits the Allee effect and the slowly varying periodic harvesting. The slow periodicity of the harvesting term may arise in response to relatively slow fluctuations of the environment. This assumption allows us to utilize the perturbation approach based on the small variation to obtain the analytical approximate solutions for the population evolution in the situation of subcritical and the supercritical harvesting with respect to the Allee threshold condition. In addition, due to the absence of the exact solution of such a system, we obtain the numerical computation using a fourth-order Runge-Kutta technique to validate the analytical results. The comparison between the two methods shows a remarkable agreement as $t\rightarrow \infty$, while such an agreement falls off when the time is close to the initial density. In addition, we show that setting $\epsilon = 0$ extinguishes the slowly varying term and the exact solution will be available in an implicit form. The finding results are appropriate for a wide range of parameter values, which leads to avoiding the extensive recalculations while displaying the population behavior.

Use of AI tools declaration

The author has not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.